Discrete Fourier Transform离散傅里叶变换算法

DFT

公式：

A k = ∑ n = 0 N − 1 e − i 2 π N k n a n , k ∈ { 0 , 1 , . . . N − 1 } A_k = \sum_{n=0}^{N-1}e^{-i\frac{2\pi}{N}kn}a_n , k\in \{0,1,... N-1\} Ak​=n=0∑N−1​e−iN2π​knan​,k∈{0,1,...N−1}

一般我们通常写成
A k = ∑ n = 0 N − 1 W N k n a n , k ∈ { 0 , 1 , . . . N − 1 } A_k = \sum_{n=0}^{N-1}W_{N}^{kn}a_n , k\in \{0,1,... N-1\} Ak​=n=0∑N−1​WNkn​an​,k∈{0,1,...N−1}

其中

$$W_N=e^{-i\frac{2\pi }{N}}$$

$W_N^k$称为旋转因子, twiddle factor，之所以称之为旋转因子，是因为它在单位圆上根据k值不同而做不同角度的旋转。 (更多阅读请参考欧拉公式和傅里叶变换原理)

下面是N=2,N=4,N=8的旋转，

从图中可以看到，旋转因子具有如下性质

$$\begin{gathered} W_N^0=1 \\ {W}_N^{\frac{N}{2}}=-1 \\ {W}_N^k=W_N^{k+N} \end{gathered}$$

---

iDFT

公式：
a n = 1 N ∑ k = 0 N − 1 W N − k n A k , k ∈ { 0 , 1 , . . . N − 1 } a_n = \frac{1}{N}\sum_{k=0}^{N-1}W_N^{-kn}A_k, k\in \{0,1,... N-1\} an​=N1​k=0∑N−1​WN−kn​Ak​,k∈{0,1,...N−1}

---

具体示例

2个sample的DFT

旋转因子

$$W_2=e^{-i2\pi /2}=e^{-i\pi }=cos\pi -isin\pi =-1$$

$A_k$:
$$A_k=\sum _{n=0}^1(-1{)}^{kn}{a}_n=(-1{)}^{0k}{a}_0+(-1{)}^{1k}{a}_1$$

所以
$$\begin{gathered} A_0=a_0+a_1 \\ {A}_1=a_0-a_1 \end{gathered}$$

4个sample的DFT

旋转因子

$$W_4=e^{-i2\pi /4}=e^{-i\pi /2}=-i$$

$A_k$:
$$A_k=\sum _{n=0}^3(-1{)}^{kn}{a}_n=(-i{)}^{0k}{a}_0+(-i{)}^{1k}{a}_1+(-i{)}^{2k}{a}_2+(-i{)}^{3k}{a}_3$$

所以
$$\begin{gathered} A_0=a_0+a_1+a_2+a_3 \\ {A}_1=a_0-i{a}_1-a_2+i{a}_3 \\ {A}_2=a_0-a_1+a_2-a_3 \\ {A}_3=a_0+i{a}_1-a_2-i{a}_3 \end{gathered}$$

可以把它抽取出来参数，作为一个矩阵相乘，
$$\left\{\begin{array}{c}{A}_0\\ A_1\\ A_2\\ A_3\end{array}\right\}=\left\{\begin{array}{cccc}1& 1& 1& 1\\ 1& -i& -1& i\\ 1& -1& 1& -1\\ 1& i& -1& -i\end{array}\right\}\left\{\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\end{array}\right\}$$

为了快速计算$A$，我们可以将计算稍作改变：
$$\begin{gathered} A_0=(a_0+a_2)+(a_1+a_3) \\ {A}_1=(a_0-a_2)-i(a_1-i{a}_3) \\ {A}_2=(a_0+a_2)-(a_1+a_3) \\ {A}_3=(a_0-a_2)+i(a_1-a_3) \end{gathered}$$
这样很多运算结果可以公用。

如果使用图线的方式来表示乘与加，有如下表示方法，

所以4个sample的DFT可以画成如下，

这便是DFT中蝴蝶算法的基本原理。
下面是8个sample的运算图

---

自己实现DFT

根据上述算法，在matlab上做如下DFT实现

定义DFT函数

function [f, Xk] = mydft(xn, fs, N)
%UNTITLED2 此处显示有关此函数的摘要
%   此处显示详细说明
%   xn:  the discrete samples
%   fs:  frequency
%   N:   the amount of samples
    n = 0:1:N-1;
    k = n;
    WN= exp(-1i*2*pi/N);
    nk = n' * k;
    Xk = xn(1:N) * WN.^nk;
    f = 0:fs/N:fs/N*(N-1);
end

Demo代码如下，

可以看到和Matlab内置函数fft结果是一样的

复数域的三维观测

频率在时间轴上形成了“波”形的二维形状，但是它仅仅展示出了振幅（能量）信息。
傅里叶变化增加了复数域轴，使它展示出了相位信息，每一个点也成为了三维空间上的点。

下面我们在三维视角下观测，操作步骤如下，

绘制波形

N=200;

fs=20;
i=0:1:N-1;
xn = sin(2*pi*i*fs/N);

绘制振幅谱

Xk = fft(xn,N);
mag = abs(Xk);

plot(1:1:N/2, mag(1:N/2));

绘制real部分

scatter(1:1:N, real(Xk));

绘制img部分

scatter(1:1:N, imag(Xk));

投影图

scatter(real(Xk), imag(Xk));

绘制3D点

scatter3(1:1:N, real(Xk), imag(Xk));

有一些有意思的结论，

1. real部分是对称的，img部分是共轭兑成的
2. 如果是完整的sin和cos波，在频率点处，可以看到跳点，原来的“直线”像断在这里一样
3. 如果是完整的sin和cos波，投影图是趋向于X型的

---

选取一段钢琴音，观测它的3D图

real部分图

imagine部分图

投影图基本呈现X形状

---

附录

8sample DFT 样例

旋转因子

$$W_8^1=e^{-i2\pi /8}=e^{-i\pi /4}=\frac{\sqrt{2}}{2}-i\frac{\sqrt{2}}{2}$$
$$W_8^2=-i$$
$$W_8^3=-\frac{\sqrt{2}}{2}-i\frac{\sqrt{2}}{2}$$
$$W_8^4=-1$$
$$W_8^5=-\frac{\sqrt{2}}{2}+i\frac{\sqrt{2}}{2}$$
$$W_8^6=i$$
$$W_8^7=\frac{\sqrt{2}}{2}+i\frac{\sqrt{2}}{2}$$
$$W_8^8=1$$

旋转因子矩阵

$$\left\{\begin{array}{c}{W}_8^0{W}_8^0{W}_8^0{W}_8^0{W}_8^0{W}_8^0{W}_8^0{W}_8^0\\ \\ W_8^0{W}_8^1{W}_8^2{W}_8^3{W}_8^4{W}_8^5{W}_8^6{W}_8^7\\ \\ W_8^0{W}_8^2{W}_8^4{W}_8^6{W}_8^0{W}_8^2{W}_8^4{W}_8^6\\ \\ W_8^0{W}_8^3{W}_8^6{W}_8^1{W}_8^4{W}_8^7{W}_8^2{W}_8^5\\ \\ W_8^0{W}_8^4{W}_8^0{W}_8^4{W}_8^0{W}_8^4{W}_8^0{W}_8^4\\ \\ W_8^0{W}_8^5{W}_8^2{W}_8^7{W}_8^4{W}_8^1{W}_8^6{W}_8^3\\ \\ W_8^0{W}_8^6{W}_8^4{W}_8^2{W}_8^0{W}_8^6{W}_8^4{W}_8^2\\ \\ W_8^0{W}_8^7{W}_8^6{W}_8^5{W}_8^4{W}_8^3{W}_8^2{W}_8^1\\ \end{array}\right\}$$