Consider the discrete Fourier transform of order $N$ defined by the following equations :
$$\begin{cases} \displaystyle y_{k}=\sum_{n=0}^{N-1} Y_{n} \omega_{N}^{n k},& \quad k=0,1,2, \ldots, N-1, \\\\ \displaystyle Y_{n}=\frac{1}{N} \sum_{k=0}^{N-1} y_{k} \omega_{N}^{-n k},& \quad n=0,1,2, \ldots, N-1, \end{cases} \tag{1}$$
where $\displaystyle \omega:=e^{i\frac{2\pi}{N}}$.
I was wondering if it's possible to write the DFT in matrix form $Y=S_{n}y$. What I noticed is that the matrix used to describe the DFT $(y_{k})\xrightarrow{\mathscr{F}_{N}}{(Y_{n})}$ is given by : $$ \Omega_{N}=\left(\omega_{N}^{n k}\right)=\left[\begin{array}{ccccc} 1 & 1 & 1 & \ldots & 1 \\ 1 & \omega_{N} & \omega_{N}^{2} & \ldots & \omega_{N}^{N-1} \\ 1 & \omega_{N}^{2} & \omega_{N}^{4} & \ldots & \omega_{N}^{2(N-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & \omega_{N}^{N-1} & \omega_{N}^{2(N-1)} & \ldots & \omega_{N}^{(N-1)^{2}} \end{array}\right] $$ However, I am not able to obtain the form $Y=S_{n}y$.
Aha, my favourite subject. There is a whole paper on the matrix identities of FFT:
Rose, Donald J., Matrix identities of the fast Fourier transform, Linear Algebra Appl. 29, 423-443 (1980). ZBL0463.15016.
Anyway, if you don't mind, I'll be abusing some notations here, particularly the $\otimes$ symbol to denote the Kronecker product, and also $\mathrm{diag}$ to mean both diagonal matrix and block diagonal matrix. The following discusses the formulation of radix-2 FFT.
Specify $P_2^N$ to be the permutation matrix that shuffles all odd rows to the top half, and even rows to bottom half (if you use 0-indexing, change odd to even and vice versa). This permutation matrix has the property $$ (P_2^N)^{\mathsf T}P_2^N = I_N $$
Applying this to $\Omega_N$ we get the following factorisation
$$ \begin{aligned} P_2^N\Omega_N &=\begin{bmatrix} \Omega_{N/2} &\Omega_{N/2}\\ \Omega_{N/2}W_N & -\Omega_{N/2}W_N \end{bmatrix}\\ &= \begin{bmatrix} \Omega_N & 0\\ 0 & \Omega_N \end{bmatrix}\begin{bmatrix} I_{N/2} & I_{N/2}\\ I_{N/2}W_N & -I_{N/2}W_N \end{bmatrix}\\ &= \begin{bmatrix} \Omega_N & 0\\ 0 & \Omega_N \end{bmatrix}\begin{bmatrix} I_{N/2} & 0\\ 0 & W_N \end{bmatrix}\begin{bmatrix} I_{N/2} & I_{N/2}\\ I_{N/2} & -I_{N/2} \end{bmatrix}\\ &=(I_2\otimes \Omega_{N/2})\mathrm{diag}(I_{N/2},W_N)(\Omega_2\otimes I_{N/2})\\ \therefore \Omega_N&=(P_2^N)^{\mathsf T}(I_2\otimes \Omega_{N/2})\mathrm{diag}(I_{N/2},W_N)(\Omega_2\otimes I_{N/2}) \end{aligned} $$ with $$ \begin{aligned} W_N&=\mathrm{diag}(1,\omega_N^1,\omega_N^2,\dots,\omega_N^{N-1}) \\ \Omega_2&=\begin{bmatrix}1 & 1 \\ 1 & -1\end{bmatrix} \end{aligned} $$ The definition can be applied recursively to $F_{N/2}$. This is the so-called decimation in frequency (DIF) radix-2 FFT.
The main ingredient to allowing the effective savings of multiplication count is the fact that
$$ w_N^{n+N/2} = -w_N^{n} $$
which is the reason why the lower right block matrix in the first line of the equation above is simple the negative of the lower left block matrix. This fact is applicable to even-numbered radices.
Meanwhile, the transpose of $\Omega_N$ yields itself $$ (\Omega^N)^{\mathsf T} = \Omega^N $$
so if you shuffle the columns instead, you get the decimation in time (DIT) formulation of radix-2 FFT:
$$ \Omega_N = (\Omega_2\otimes I_{N/2})\mathrm{diag}(I_{N/2},W_N)(I_2\otimes \Omega_{N/2})P_2^N $$
In either DIT or DIF case, if the same decimation method is used for all recursions, the permutation matrix all appear on the same side.
$$ (I_{N/2}\otimes P_2^2)\dots(I_2\otimes P_2^{N/2})P_2^N $$
Altogether, this forms the bit-reversal permutation.
We now have a clue on how to formulate decimations of other higher radices. For example, for radix-4, you need a different permutation $P_4^N$ that sorts $y_n$ into four row blocks $y_{4n}, y_{4n+1}, y_{4n+2}, y_{4n+3}$. The formula is then
$$ \begin{aligned} \Omega_N &=(P_4^N)^{\mathsf T}(I_4\otimes \Omega_{N/4})\mathrm{diag}(I_{N/4},W_N,W_N^2,W_N^3)(\Omega_4\otimes I_{N/4}) &(\textrm{DIF})\\ &=(\Omega_4\otimes I_{N/4})\mathrm{diag}(I_{N/4},W_N,W_N^2,W_N^3)(I_4\otimes \Omega_{N/4})P_4^N &(\textrm{DIT}) \end{aligned} $$
where $$ W_N^k = \mathrm{diag}(1,\omega_N^k,\omega_N^{2k},\dots,\omega_N^{N-k}) $$
Check out also the following paper for another interesting approach:
Zhou, Yicong; Cao, Weijia; Liu, Licheng; Agaian, Sos; Chen, C. L. Philip, Fast Fourier transform using matrix decomposition, Inf. Sci. 291, 172-183 (2015). ZBL1355.65188.
(On the other hand, as $\Omega_N$ is already orthogonal, I doubt you would get any interesting results through QR or SVD methods.)