# A Sparse Algorithm for Computing the DFT Using Its Real Eigenvectors
## 4 Computing Orthonormal Eigenvectors of the DFT
### Matveev Method  
The first step in the Matveev method is to compute the four projection matrices of the DFT corresponding to each unique eigenvalue. If \(P_{e}\) is a projection matrix corresponding to an eigenvalue \(e\), \(F\) is the Fourier transform matrix and \(I\) is the identity matrix,  
\[ \begin{split} P_{1}&=\frac{1}{4}\Big{(}F^{3}+F^{2} +F+I\Big{)}\\ P_{-1}&=\frac{1}{4}\Big{(}-F^{3}+F^{2}-F+I\Big{)} \\ P_{j}&=\frac{1}{4}\Big{(}jF^{3}-F^{2}-jF+I\Big{)} \\ P_{-j}&=\frac{1}{4}\Big{(}-jF^{3}-F^{2}+jF+I\Big{)} \end{split}  \tag{18}\]  
The rank of each of these matrices is tied to their corresponding eigenvalue multiplicity. Thus, the matrix \(P_{e}\) will be rank \(m_{e}\). Therefore, to compute the orthogonal eigenvectors only the first \(m_{1}\) and \(m_{-1}\) columns of \(P_{1}\) and \(P_{-1}\) are considered; the first column of \(P_{j}\) and \(P_{-j}\) are ignored and only the next \(m_{j}\) and \(m_{-j}\) columns are considered. The first columnsof \(P_{j}\) and \(P_{-j}\) are ignored because they are zeros. The kept columns form independent eigenvectors of \(F\) corresponding to each respective eigenvalue. Let the elements of the matrices in (18) be represented as \(P_{e}(m,n)\). For example, if the indexing starts from 0, the element in the second row and first column of \(P_{-1}\) will be \(P_{-1}(1,0)\). Therefore, the last element will be \(P_{-1}(N-1,N-1)\). Let the columns of these matrices be represented as \(p_{e,k}\), where \(p_{e,k}\) is the \(k^{\text{th}}\) column vector of \(P_{e}\). For example, the third column vector of \(P_{j}\) will be \(p_{j,2}\), again assuming the indexing starts from 0. Note that these projection matrices are symmetric. Therefore, the \(k^{\text{th}}\) column will also be the \(k^{\text{th}}\) row. The orthogonal eigenvectors are given by the determinants of the following matrices. For eigenvalues \(e=\{1,-1\}\),  
\[\beta_{e,0}=p_{e,0}\]  
\[\beta_{e,1}=\left| \begin{array}{cccc}P_{e}(0,0)&p_{e,0}\\ P_{e}(1,0)&p_{e,1}\end{array} \right|\]  
\[\beta_{e,k}=\left| \begin{array}{cccc}P_{e}(0,0)&P_{e}(0,1)&\cdots&P_{e}(0,k- 1)&p_{e,0}\\ P_{e}(1,0)&P_{e}(1,1)&&P_{e}(1,k-1)&p_{e,1}\\ \vdots&&\ddots&\vdots&\vdots\\ P_{e}(k,0)&P_{e}(k,1)&\cdots&P_{e}(k,k-1)&p_{e,k}\end{array} \right|\]  
where \(\beta_{e,k}\) form a set of orthogonal eigenvectors corresponding to eigenvalues \(e=\{1,-1\}\) and \(1\leq k\leq m_{e}-1\). \(\beta_{e,1}\) gives the eigenvector when \(k=1\). Because these matrices contain scalars \(P_{e}(m,n)\) and vectors \(p_{e,k}\), the determinant can be computed by expanding along the last column and multiplying these vectors by their corresponding cofactors. Similarly, for eigenvalues \(e=\{j,-j\}\),  
\[\beta_{e,0}=p_{e,1}\]  
\[\beta_{e,1}=\left| \begin{array}{cccc}P_{e}(1,1)&p_{e,1}\\ P_{e}(2,1)&p_{e,2}\end{array} \right|\]  
\[\beta_{e,k}=\left| \begin{array}{cccc}P_{e}(1,1)&P_{e}(1,2)&\cdots&P_{e}(1,k )&p_{e,1}\\ P_{e}(2,1)&P_{e}(2,2)&&P_{e}(2,k)&p_{e,2}\\ \vdots&&\ddots&\vdots&\vdots\\ P_{e}(k+1,0)&P_{e}(k+1,2)&\cdots&P_{e}(k+1,k)&p_{e,k+1}\end{array} \right|\]  
where \(\beta_{e,k}\) is the general equation. Here also \(1\leq k\leq m_{e}-1\). Note that if \(e=\{j,-j\}\), for the projection matrix \(P_{e}\) the elements of the first row \(P_{e}(0,k)\), first column \(P_{e}(k,0)\) and the first column vector \(p_{e,0}\) are not considered since they are zeros. The determinant for these matrices can be computed using the method explained in the previous paragraph for the real eigenvalues. In both cases, the resulting vectors \(\beta_{e,k}\) are orthogonal but not unit length. After normalization, the columns of the eigenvector matrix matrix \(V\) are  
\[v_{e,k}=\frac{\beta_{e,k}}{\left\|\beta_{e,k}\right\|}\]  
Although (18) contains powers of the Fourier transform matrix, it is not necessary to numerically calculate those powers. Using the properties of the DFT, \(F^{3}=F^{\text{H}}\), where \(F^{\text{H}}\) represents the complex conjugate transpose of \(F\). \(F^{2}\) will be a real matrix such that \(F^{2}=\left[