快速傅里叶变换FFT

1. 离散傅里叶变换

时域离散线性时不变系统理论和离散傅里叶变换是数字信号处理的理论基础，数字滤波和数字谱分析是数字信号处理的核心。快速傅里叶变换并不是一种新的变换，而是离散傅里叶变换（Discrete Fourier Transform，DFT）的一种高效算法。

如果信号在频域是离散的，则该信号在时域就表现为周期性的时间函数；相反，在时域上是离散的，则该信号在频域必然表现为周期性的频率函数。不难设想，如果时域信号不仅是离散的，而且是周期的，那么由于它时域离散，其频谱必是周期的，响应的频谱必是离散的。换句话说，一个离散周期时间序列，它一定具有既是周期又是离散的频谱。我们还可以得出一个结论：一个域的离散必然造成另一个域的周期延拓，这种离散变换，本质上都是周期的。

一个连续信号经过理想采样后的表达式为
$$x_a(t)=\sum _{-\infty }^{\infty }{x}_a(nT)\delta (t-nT)$$
其频谱函数时是上式的傅里叶变换，容易得出其傅里叶变换为
$$X_a(j\Omega )=\sum _{-\infty }^{\infty }{x}_a(nT)e^{-j\Omega nT}$$
式中，$\Omega$为模拟角频率，单位为弧度/秒，它与数字角频率$\omega$之间的关系为$\omega =\Omega T$。对于数字信号来说，处理的信号其实是一个数字序列。因此可以用$x(n)$代替$x_a(nT)$，同时用$X(e^{j\omega }))$代替$X_a(j\omega /T)$，则可以得到时域离散信号的频谱表达式，即
$$X(e^{j\omega })=\sum _{-\infty }^{\infty }x(n)e^{-j\omega n}$$
显然，$X(e^{j\omega })$是以周期为$2\pi$的周期函数。

对于一个长度为N的有现场序列，在频域表现为周期性的连续谱$X(e^{j\omega })$。如果我们将有限长序列以周期N进行周期性拓展，则在频域必将表现为周期性的离散谱$X(e^{j{\omega }_s})$，且单个周期的频谱形状与有限长序列相同。因此，我们可以将$X(e^{j{\omega }_s})$看做在频域对$X(e^{j\omega })$等间隔采样的结果。根据采样理论，要想采样后能够不失真地恢复出原信号，采样速率必须满足一定的条件。假设时域信号的周期长度为NT，则在频域的一个周期内，采样点数$N_0$必须大于等于N。

用离散角频率变量$k{\omega }_s$代替$X(e^{j{\omega }_s})$中连续变量${\omega }_s$，且取$N_0=N$，则有限长序列的频谱表达式为
$$X(e^{j{\omega }_s})=\sum _{n=0}^{N-1}x(n)e^{-j(2\pi /N)kn}$$
令以N为周期的周期函数$W_N^{kn}=e^{-j(2\pi /N)kn},\tilde{X}=X(e^{jk{\omega }_s})$，$\tilde{x}(n)$为序列$x(n)$以N为周期的延拓，可以写为
$$\tilde{X}(k)=\sum _{n=0}^{N-1}\tilde{x}(n)W_N^{kn}$$
将式两边同乘以${\sum }_{n=0}^{N-1}{W}_N^{-kn}$，可以得到
$$x(n)=(1/N)\sum _{n=0}^{N-1}\tilde{X}(K)W_N^{-kn}$$
由于上两式中，n和k有这样的关系$0⩽k⩽N-1,0⩽n⩽N-1$，也就是说只涉及一个周期内的N个样本。因此，也可以用有限长序列$x(n)$和$x(K)$，即各取一个周期来表示这些关系式。我们定义有限长序列$x(n)$和$x(K)$之间的关系为离散傅里叶变换（DFT）。
$$X(k)=\tilde{X}(k)R_N(k)=\sum _{n=0}^{N-1}x(n)W_N^{kn},0⩽k⩽N-1$$
$$x(n)=\tilde{x}(n)R_N(n)=(1/N)\sum _{n=0}^{N-1}X(k)W_N^{-kn},0⩽n⩽N-1$$

2. DFT存在的问题

DFT是分析信号频谱的有利工具，在应用DFT分析连续信号的频谱时，其中涉及序列补零、混叠失真、频谱泄漏和栅栏效应等问题。

2.1 栅栏效应和序列补零

用DFT计算频谱，只能给出频谱的${\omega }_k=2\pi k/N或{\Omega }_k=2\pi k/NT$，即频谱的采样值，而不可能得到连续的频谱函数。就好像通过一个栅栏看信号频谱，只能在离散点上看到信号频谱，这种现象被称为栅栏效应。

在DFT变换计算过程中，如果序列长为N个点，则只要计算N点DFT。这意味着对序列$x(n)$的傅里叶变换在$(0,2\pi )$区间只计算N个点的值，其频率采样间隔为$2\pi /N$。如果序列长度较小，频率采样间隔${\omega }_s=2\pi /N$可能太大，以至于不能直观地说明信号的频谱特性。然而，有一种非常简单的方法能解决这一问题，这种方法能对序列的傅里叶变换以足够小的间隔进行采样。数字频率间隔$\Delta {\omega }_k=2\pi /L$，这里的L是DFT变换的点数，显然，要提高数字频率间隔，只需增加L即可。当序列长度N较小时，可以采用在数字序列后面增加L-N个零值的办法，对L点序列进行DFT变换，满足所需的频谱采样间隔。这样做可以保证在原来频谱形状不变的情况下，使谱线加密，即使频域采样点数增加，从而使原来看不到的频谱分量可以看到。

注意，补零可以改变频谱密度，但窗函数的宽度不能改变。换句话说，必须按照数据记录的有效长度选择窗函数，而不能按补零后的长度来选择窗函数。

2.2 频谱泄漏和混叠失真

对信号进行DFT计算，首先必须使其变成时宽有限的信号，方法是将序列$x(n)$与时宽有限窗函数$\omega (n)$相乘。例如，选用矩形窗函数来截断信号，在频域中则相当于信号的频谱与窗函数频谱的周期卷积。卷积将造成频谱失真，且这种失真主要表现在原频谱的展宽，这个现象即称为频谱泄漏。因为泄漏将导致频谱扩展，从而使信号最高频谱可能超过采样频率的一半，造成混叠失真。

在进行DFT运算时，时域截短是必要的，因而频谱泄漏是不可避免的。为尽量减少泄漏的影响，可以用适当形状的窗函数，如海明窗、汉宁窗等。需要注意的是，在进行DFT之前，预加窗函数可以改善频谱泄漏情况，但必须对数据进行重叠处理以补偿窗函数边缘对数据的衰减，通常采用汉明窗进行50%重叠的处理。

2.3 频率分辨率与DFT参数的选择

在对信号进行DFT变换，分析信号的频谱特征时，通常采用频率分辨率来表征在频率轴上所能得到的最小频率间隔。对于长度为N的DFT变换，其频率分辨率$\Delta f=f_s/N$，其中$f_s$为时域信号的采样频率。需要注意的是，这里的数据长度N必须是数据的有效长度。如果在$x(n)$中有两个频率分别是$f_1、f_2$的信号，则在对$x(n)$用矩形窗进行截断时，要分辨这两个频率，必须满足
$$2f_s/N<\left|f_1-f_2\right|$$
DFT变换时的补零没有增加序列的有效长度，所以并不能提高分辨率；但补零可以使数据N为2的整数次幂，以便于使用接下来要介绍的快速傅里叶变换算法。补零对原$X(k)$起到插值作用，一方面克服栅栏效应，平滑谱的外观；另一方面，由于数据截短时引起的频谱泄漏，有可能在频谱中出现一些难以确认的谱峰，补零后有可能消除这种现象。

3. FFT算法

通常 $x(n)、X(k)$和$W_N^{kn}$都是复数，因此每计算一个$X(k)$值，必要要进行N次复数乘法和N-1次复数加法。而$X(k)$共有N个值，所以要完成全部DFT的运算要进行$N^2$次复数乘法和N(N-1)次复数加法。乘法运算比加法运算复杂，且运算时间更长，所耗资源更多，因此可以用乘法运算来衡量一个算法的运算量。由于复数乘法最终还得通多实数乘法运算完成，而每个复数乘法需要4个实数乘法运算，因此完成全部DFT运算需要进行$4N^2$次实数乘法运算。

因此，直接计算DFT，乘法次数与$N^2$成正比。随着N的增大，运算次数迅速增加。例如，当N=8时，需要64次复数乘法，而当N=1024时，则要1048576次复数乘法，即100多万次复数乘法。如果信号处理要求实时进行，对计算速度的要求实在是太高了。

仔细观察DFT和IDFT的运算，会发现系数$W_N^{kn}$具有对称性和周期性，即
$$(W_N^{kn}{)}^{\ast }=W_N^{-kn}$$
$$W_N^{n(N+k)}=W_N^{k(N+n)}=W_N^{kn}$$
$$W_N^{-nk}=W_N^{b(N-k)}=W_N^{k(N-n)}$$
$$W_N^{N/2}=-1,则W_N^{k+N/2}=-W_N^k$$
利用$W_N^{kn}$的周期性，在DFT运算中有些项目可以合并，从而减少运算量。又由于DFT的运算量与$N^2$成正比，因此N越小越有利，我们可以利用对称性和周期性将大点数的DFT分解成很多小点数的DFT。FFT算法正是基于这样的基本思路。为了不断地进行分解，FFT算法要求DFT的运算点数为$N=2^M$，M为正整数。这种N为2的整数幂的FFT，称为基-2FFT。
FFT算法可分为两大类：按时间抽取（Decimation In Time，DIT）和按频率抽取（Decimation In Frequency，DIF）。为提高运算速度，将DFT的计算逐次分解成较小点数的DFT。如果算法是逐次分解时间序列$x(n)$得到的，称为按时间抽取的FFT算法；如果算法是逐次分解频域序列$X(k)$得到的，称为按频域抽取的FFT算法。
FFT算法运算量约为$(N/2)lo{g}_{2}N$次复数乘法运算，当N较大时，其运算速度比DFT大大提高。

4. FFT算法的MATLAB仿真

分别取128点的$x(n)$、将128点$x(n)$以补零的方式加长到512点、取512点$x(n)$计算FFT。

源码

f1=2; f2=2.05;%单正弦信号的频率
fs=10;        %采样频率
%对128点时域序列进行FFT分析
N=128;        %FFT分析的点数
n=0:N-1;
xn1=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs);%产生128点时域信号序列
XK1=fft(xn1);                           %进行傅立叶变换，并进行归一化处理
MXK1=abs(XK1(1:N/2));
% MXK1=20*log10(abs(XK1(1:N/2)));
% MXK1=MXK1-max(MXK1);
%对补零后的512点时域序列进行FFT分析
M=512;
xn2=[xn1 zeros(1,M-N)];%在时域信号序列后补零
XK2=fft(xn2);          %进行傅立叶变换，并进行归一化处理
MXK2=abs(XK2(1:M/2));
%MXK2=20*log10(abs(XK2(1:M/2)));
%MXK2=MXK2-max(MXK2);
%对512点时域序列进行FFT分析
n=0:M-1;
xn3=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs);%产生128点时域信号序列
XK3=fft(xn3);                           %进行傅立叶变换，并进行归一化处理
MXK3=abs(XK3(1:M/2));
% MXK3=20*log10(abs(XK3(1:M/2)));
% MXK3=MXK3-max(MXK3);
%绘图
subplot(321);
x1=0:N-1; 
plot(x1,xn1);xlabel('n','fontsize',8);title('128点x(n)','fontsize',8);
subplot(322);
k1=(0:N/2-1)*fs/N;
plot(k1,MXK1);xlabel('f(Hz)','fontsize',8);title('128点xn的FFT变换','fontsize',8);
subplot(323);
x2=0:M-1; 
plot(x2,xn2);xlabel('n','fontsize',8);title('512点补零x(n)','fontsize',8);
subplot(324);
k2=(0:M/2-1)*fs/M;
plot(k2,MXK2);xlabel('f(Hz)','fontsize',8);title('512点补零xn的FFT变换','fontsize',8);
subplot(325);
plot(x2,xn3);xlabel('n','fontsize',8);title('512点x(n)','fontsize',8);
subplot(326);
plot(k2,MXK3);xlabel('f(Hz)','fontsize',8);title('512点xn的FFT变换','fontsize',8);

仿真结果

显然，补零对分辨率没有影响，只是对频谱图起到平滑作用。并且128点$x(n)$不满足频率分辨率要求，没法区分出两种频率成分，而512点中可明显区分出来。