一窥数字傅里叶变换矩阵

前言

DFT技术虽然很神奇，但总觉得是个黑洞，摸不着头脑，今天研究一下DFT matrix，能更加清楚的了解一些内核的东西。

DFT vs IDFT matrics

这是一个N*N方阵，矩阵的基本元素由$W=e^{-j2\pi /N}$来构成，每个位置元素$a_{ij}=W^{(i-1)\ast (j-1)}=e^{-j2\pi (i-1)(j-1)/N}$。即矩阵元素的索引是从自然数1开始的，但指数运算需要从0（据说现在小学已经把0化为自然数了。。。）开始。简单的列一下这个矩阵：$$\begin{gathered} \mathbf{W}=\left[\begin{array}{ccccc}{e}^{-j0}& e^{-j0}& e^{-j0}& \cdots & e^{-j0}\\ e^{-j0}& e^{-j2\pi /N}& e^{-j4\pi /N}& \cdots & e^{-j2(N-1)\pi /N}\\ e^{-j0}& e^{-j4\pi /N}& e^{-j8\pi /N}& \cdots & e^{-j4(N-1)\pi /N}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ e^{-j0}& e^{-j2(N-1)\pi /N}& e^{-j4(N-1)\pi /N}& \cdots & e^{-j2(N-1)(N-1)\pi /N}\end{array}\right] \\ =\left[\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& e^{-j2\pi /N}& e^{-j4\pi /N}& \cdots & e^{-j2(N-1)\pi /N}\\ 1& e^{-j4\pi /N}& e^{-j8\pi /N}& \cdots & e^{-j4(N-1)\pi /N}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& e^{-j2(N-1)\pi /N}& e^{-j4(N-1)\pi /N}& \cdots & e^{-j2(N-1)(N-1)\pi /N}\end{array}\right] \end{gathered}$$
这个矩阵是主对角对称矩阵。即${\mathbf{W}}^T=\mathbf{W}$

酉矩阵特性

上述矩阵满足$${\mathbf{W}}^{†}\mathbf{W}=NI$$
即厄米特共轭与原矩阵的乘积是N倍的单位阵。具有酉矩阵的特点（差了N倍，不知道失算uniary，还是paraunitary）。上面的式子还暗示了${\mathbf{W}}^{†}/N={\mathbf{W}}^{-1}$，另外一个推论是${\mathbf{W}}^{†}={\mathbf{W}}^{\ast }$

矩阵DFT运算

假设一组序列为$X=[x(0),x(1),...x(N-1){]}^T$，这组序列的DFT变换$$X(k)=\sum _{n=0}^{N-1}x(n)e^{-j2\pi kn/N}=\sum _{n=0}^{N-1}x(n)W^{kn}=\mathbf{W}X$$这里的n和k都是0开始的了。$$x(n)=\frac{1}{N}\sum _{n=0}^{N-1}X(k)e^{j2\pi kn/N}=\frac{1}{N}\sum _{n=0}^{N-1}X(k)W^{-kn}$$

矩阵的值长什么样

我们利用python代码来一窥DFT矩阵的样子

    pi = math.pi
    M = 16
    W = np.exp((2*pi/M)*1j)
    nature = np.arange(M)
    DFT_Matrix = np.ones([M,M],np.complex)
    for row in range(M):
        DFT_Matrix[row]=W**(-nature*(row))

当M=2时，得到如下的值，最后一个值$e^{-j2\pi /2}=-1$

 [1.+0.j 1.+0.j]
 [ 1.+0.0000000e+00j -1.-1.2246468e-16j]

M=4时，

[1.+0.j 1.+0.j 1.+0.j 1.+0.j]
[ 1.0000000e+00+0.0000000e+00j  6.1232340e-17-1.0000000e+00j\n -1.0000000e+00-1.2246468e-16j -1.8369702e-16+1.0000000e+00j]
[ 1.+0.0000000e+00j -1.-1.2246468e-16j  1.+2.4492936e-16j\n -1.-3.6739404e-16j]
[ 1.0000000e+00+0.0000000e+00j -1.8369702e-16+1.0000000e+00j\n -1.0000000e+00-3.6739404e-16j  5.5109106e-16-1.0000000e+00j]

操练

根据【2】提供的图形，对256个点进行DFT变换，将时域频域画在一张图里，就差不多这个样子了。图中每个算法采样点换成256，然后增加几行代码：

    Y = np.matmul(DFT_Matrix, y)
    y_idx = np.linspace(0, 2 * pi, 256)
    plt.plot(y_idx, np.absolute(Y))

纯属娱乐，但这个工具对于不要求运算速度的话，还是很直观的。

总结

这个小算法能够比较清晰的将DFT算子展现在你面前，帮助理解DFT的玄妙之处

参考

1.Multirate Systems and Filter Banks, P.Vaidyanathan
2.Python信号分析 | 信号的表示（二）【三角、复指数、矩形脉冲、阶跃】