自己动手实现fft.m函数

背景

自己动手使用matlab实现fft.m函数

知识回顾

1、四种形式的Fourier变换

首先，我们需要系统回忆一下4种形式的Fourier变换

1)、时域连续，周期信号

$$x(t)=\sum _{k=-\infty }^{\infty }X(k{\Omega }_0)e^{jk{\Omega }_{0}t}$$
$$X(k{\Omega }_0)=\frac{1}{T}{\int }_{T}x(t)e^{-jk{\Omega }_{0}t}dt$$
其中。${\Omega }_0=\frac{2\pi }{T}$，T为周期。

2)、时域连续，非周期信号

$$X(j\Omega )={\int }_{-\infty }^{\infty }x(t)e^{-j\Omega t}dt$$
$$x(t)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }X(j\Omega )e^{j\Omega t}d\Omega$$

3)、时域离散，周期信号

$$X(k)=\sum _{n=0}^{N-1}x(n)e^{-j\frac{2\pi }{N}nk}$$
$$x(n)=\frac{1}{N}\sum _{k=0}^{N-1}X(k)e^{j\frac{2\pi }{N}nk}$$
其中，N为离散信号的周期

4)、时域离散，非周期信号

$$X(e^{jw})=\sum _{n=-\infty }^{\infty }x(n)e^{-jwn}$$
$$x(n)=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }X(e^{jw})e^{jwn}dw$$
以上4种情况的傅里叶变换，最为常用的是第3种和4种，也称之为离散傅里叶级数DFS,离散时间傅里叶变换DTFT。第1种情况和第2种情况分别被成为傅里叶级数FS和傅里叶变换FT。我们今天的思路也是由DFS介绍到DFT再最后介绍FFT。

由上述4种情况可以总结出的是，如果信号在时域是周期的，那么在频域则必定是离散的。如果在时域是非周期的，那么在频域必定是连续的。

2、离散傅里叶变换DFT

我们知道，在计算机种要实现信号的频谱分析及其他方面的处理工作时，对信号的要求应该由2点，一是时域和频域都应该是离散的，二是应该有限长。在上述4种情况中，只有DFS在时域和频域都是离散的。但是，两者却都是无限长的。我们仔细观察这个变换对
$$X(k)=\sum _{n=0}^{N-1}x(n)e^{-j\frac{2\pi }{N}nk}$$
$$x(n)=\frac{1}{N}\sum _{k=0}^{N-1}X(k)e^{j\frac{2\pi }{N}nk}$$
其中，N为离散信号的周期。
仔细观察我们可以知道，$e^{j\frac{2\pi }{N}nk}$和$e^{-j\frac{2\pi }{N}nk}$相对n和k都是以N为周期的（因为$e^{jw}$是以$2\pi$为周期的）。利用这个性质，我们将DFS重写，得到离散傅里叶变换DFT。
$$X(k)=\sum _{n=0}^{N-1}x(n)e^{-j\frac{2\pi }{N}nk}=\sum _{n=0}^{N-1}x(n)W_N^{nk}$$
其中，$k=0,1,...N-1$
$$x(n)=\frac{1}{N}\sum _{k=0}^{N-1}X(k)e^{j\frac{2\pi }{N}nk}=\frac{1}{N}\sum _{n=0}^{N-1}x(n)W_N^{-nk}$$
$$W_N=e^{-j\frac{2\pi }{N}}$$
其中，$n=0,1,...N-1$。上面，便得到了DFT。DFT对应的是在时域、频域都是有限长、且又都是离散的一类变换。显然，DFT并不是一个新的Fourier变换形式，它实际上来自于DFS，只不过仅在DFS的时域和频域各取一个周期而已。由这一个周期作延拓，就可以得到整个的DFS。

在实际工作中常常遇到的是非周期序列，它们可能是无限长，也可能是有限长。对这样的序列做Fourier变换，理论上应该是要做DTFT，但是做完DTFT之后，得到的频域是周期且连续的$X(e^{jw})$，$X(e^{jw})$无法直接在计算机上进行数字运算。那么我们该如何处理呢？

如果$x(n)$是一个有限长序列，则我们令其长度为N。如果$x(n)$是一个无限长序列，我们可以用矩形窗将其截为N点，然后把这N点视为一个周期序列$\tilde{x}(n)$的一个周期。也就是说，$\tilde{x}(n)$是由N点$x(n)$做周期延拓形成的。这样，便可对$\tilde{x}(n)$做DFT变换。

我们可以发现，求出一点的X(k)，需要N次复数乘法。计算N点X（k），则需要$N^2$次乘法。运算量很大，而快速傅里叶变换FFT算法可以将N点DFT的计算量降低到$\frac{N}{2}lo{g}_2^N$次，从而使得DFT成为信号处理中最重要也是最方便的运算。

3、FFT

最后终于到了FFT变换了。
我们知道，一次复数相乘，是需要4次实数相乘和2次实数相加的。因此，对于FFT来说，计算量确实非常大。但是，在计算DFT的过程中，实际上存在着大量的重复运算。
首先观察$W_N=e^{-j\frac{2\pi }{N}}$因子，具有如下特点，

$W^0=1$,$W^{\frac{N}{2}}=-1$
$W_N^{N+r}=W_N^r$,$W^{\frac{N}{2}+r}=-W^r$
例如，对于4点DFT，本来是需要计算4x4=16次复数乘法的。我们可以写成如下矩阵形式
$$\left[\begin{array}{c}X(0)\\ X(1)\\ X(2)\\ X(3)\end{array}\right]=\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& W^1& -1& -W^1\\ 1& -1& 1& -1\\ 1& -W^1& -1& W^1\end{array}\right]\left[\begin{array}{c}x(0)\\ x(1)\\ x(2)\\ x(3)\end{array}\right]$$
我们利用线性代数的知识，将乘法矩阵的第2列和第3列进行交换，得到
$$\left[\begin{array}{c}X(0)\\ X(1)\\ X(2)\\ X(3)\end{array}\right]=\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& -1& W^1& -W^1\\ 1& 1& -1& -1\\ 1& -1& -W^1& W^1\end{array}\right]\left[\begin{array}{c}x(0)\\ x(2)\\ x(1)\\ x(3)\end{array}\right]$$
，再将上式展开，可以得到
$$\begin{gathered} X(0)=[x(0)+x(2)]+[x(1)+x(3)] \\ X(1)=[x(0)-x(2)]+[x(1)-x(3)]W^1 \\ X(2)=[x(0)+x(2)]-[x(1)+x(3)] \\ X(3)=[x(0)-x(2)]-[x(1)-x(3)]W^1 \end{gathered}$$
经过上述的一个矩阵交换运算，成功将16次复数乘法化简为了1次复数乘法。还记得前面说的吗，降低至$\frac{2}{N}lo{g}_2^N=4$次，4次实数乘法刚好对应一次复数乘法。

上式即是FFT算法的一个雏形。问题的关键在于如何利用W因子的周期性及对称性，导出一个高效的快速算法。

总地来说，FFT算法的算法由2个发展方向，一是针对N等于2的整数次幂的算法，如基2算法，基4算法，实因子算法和分裂基算法。另一个是N不等于2的整数次幂的算法，如素因子算法等。这里介绍最经典的基2算法。

令$N=2^M$，M为正整数。我们可以将$x(n)$按奇偶分为两组，即令$n=2r$和$n=2r+1$，$r=0,1,...\frac{N}{2}-1$。
因此，我们得到
$$\begin{gathered} X(k)=\sum _{r=0}^{\frac{N}{2}-1}x(2r)W_N^{2rk}+\sum _{r=0}^{\frac{N}{2}-1}x(2r+1)W_N^{(2r+1)k} \\ =\sum _{r=0}^{\frac{N}{2}-1}x(2r)W_{N/2}^{rk}+W_N^k\sum _{r=0}^{\frac{N}{2}-1}x(2r+1)W_{N/2}^{rk} \end{gathered}$$
式中，
$$W_{N/2}=e^{-j\frac{2\pi }{N/2}}=e^{-j4\pi /N}$$
这样，经过层层分解，共经过M-1层分解，可将N全部分解完毕。这里还是不详细展开讲了。

4、fft代码

4.1 dft与idft代码实现

首先，我们回顾一下dft的计算公式
$$X(k)=\sum _{n=0}^{N-1}x(n)e^{-j\frac{2\pi }{N}nk}=\sum _{n=0}^{N-1}x(n)W_N^{nk}$$
其中，$k=0,1,...N-1$

function [Xk]=dft(xn, N);
% 计算离散傅里叶变换，输入为xn，N为输入信号的长度
n=[0:1:N-1];
k=[0:1:N-1];
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^nk;
Xk=xn*WNnk;

我们再回顾一下idft公式的实现
$$x(n)=\frac{1}{N}\sum _{k=0}^{N-1}X(k)e^{j\frac{2\pi }{N}nk}=\frac{1}{N}\sum _{n=0}^{N-1}x(n)W_N^{-nk}$$

function [xn]=idft(Xk,N);
%计算反变换
n=[0:1:N-1];
k=[0:1:N-1];
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^(-nk);
xn=(Xk*WNnk)/N;

4.2 fft代码实现

首先介绍一下matlab自带的fft函数，用法如下：
X=fft(x,N)。其中x为输入序列，N为FFT变换的长度。当缺少输入参数N时，则以序列x的长度作为FFT变换的长度。当出现x的长度大于N时，则对序列x进行截断N进行FFT变换。当出现x的长度小于N时，则对序列x进行补零操作，达到N长度后再进行FFT变换。
这里需要注意的一点是，在正式进行DFT变换之前，需要对原输入信号先进行变址运算。
代码如下

function y=myditfft(x)

% 本程序对输入序列 x 实现DIT-FFT基2算法，点数取大于等于x长度的2的幂次，因此本函数不需要特地输入N
% x为给定时间序列
% y为x的离散傅立叶变换
% 
m=nextpow2(x);N=2^m;            % 求x的长度对应的2的最低幂次m
if length(x)<N 
    x=[x,zeros(1,N-length(x))]; % 若x的长度不是2的幂，补零到2的整数幂
end  
nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1;   % 求1:2^m数列的倒序
y=x(nxd);                             % 将x倒序排列作为y的初始值
for mm=1:m             % 将DFT作m次基2分解,从左到右，对每次分解作DFT运算
    Nmr=2^mm;u=1;           % 旋转因子u初始化为WN^0=1
    WN=exp(-i*2*pi/Nmr);    % 本次分解的基本DFT因子WN=exp(-i*2*pi/Nmr)
      for j=1:Nmr/2         % 本次跨越间隔内的各次蝶形运算
        for k=j:Nmr:N       % 本次蝶形运算的跨越间隔为Nmr=2^mm
          kp=k+Nmr/2;        % 确定蝶形运算的对应单元下标
          t=y(kp)*u;        % 蝶形运算的乘积项
          y(kp)=y(k)-t;     % 蝶形运算
          y(k)=y(k)+t;      % 蝶形运算
        end
      u=u*WN;               % 修改旋转因子,多乘一个基本DFT因子WN
      end
end

function xtc_fft(x,NFFT)
%% 进行快速傅里叶变换
i=NFFT;
M=nextpow2(NFFT);
BitReverse(x,NFFT);
le=2;le2=1;
for l=1:M
    uR=1;uI=0;
    sR=cos(PI);sI=-sin(PI);
    for j=1:le2
        for i=j-1:2:NFFT-1
            ip=i+le2;
            tR=real(x(ip))*uR-imag(x(ip))*uI;
            tI=real(x(ip))*uI+imag(x(ip))*uR;
            real(x(ip))=real(x(i))-tR;
            imag(x(ip))=imag(x(i))-tI;
            real(x(i))=real(x(i))+tR;
            imag(x(i))=imag(x(i))+tI;
        end
        tR=uR;
        uR=tR*sR-uI*sI;
        uI=tR*sI+uI*sR;
    end
    le=le*2;
    le2=le2*2;
end

end

function BitReverse(x,NFFT)
%% 变址运算，把自然顺序变成倒位序.NFFT的长度一般是输入x的长度，输入x是一个复数
j=1;
k=1;
nextValue=NFFT/2;
for i=1:NFFT
    if i<j
        temp=x(j);
        x(j)=x(i);
        x(i)=temp;
    end
    k=nextValue;
    while k<=j
        j=j-k;
        k=k/2;
    end
    j=j+k;
end

end

5、关于补零问题的一些解释

5.1 分辨率

分辨率包括时间分辨率和频率分辨率。频率分辨率是通过一个频域的窗函数来观察频谱时所看到的频率的宽度，时间分辨率则是指通过一个时域的窗函数来观察频谱时所看到的频率的宽度。显然，这样的窗函数越窄，相应的分辨率就越好。

具体来说，频率分辨率指的是所用的算法能将信号中两个靠的很近的谱峰保持分开的能力。显然，频率分辨率的这一定义主要用来评价各种算法在谱分析方面的性能。我们分别针对连续时间信号的傅里叶变换，离散时间信号的傅里叶变换DTFT和离散傅里叶变换DFT来说明频率分辨率的定义。

如果信号$x_T(t)$的长度为T秒，通过傅里叶变换后得到$X_T(j\Omega )$的频率分辨率为
$$\Delta f=1/T(Hz)$$
这点比较好理解。因为长度为T秒的$x_T(t)$可以看作一个无穷长的信号$x(t)$和一宽度为T的矩形窗相乘产生的结果。该矩形窗的频谱为熟知的sinc函数。它的主瓣宽度反比于T，因此$X_T(j\Omega )$能够分辨的最小频率间隔不会小于$1/T$。

未完待续

参考资料

《数字信号处理 理论、算法与实现》胡广书编著
《数字信号处理教程 matlab释义与实现》陈怀琛著