WIN-fft结果的解释.doc

FFT变换结果的解释
作者：dhwik
FFT是每个做信号处理都无法回避的问题，大家都在用，在使用的过程中，免不了会产生一些疑惑，变换后横纵坐标的标称值表示的是什么含义，它与原始信号对应频率的幅度是什么关系，为什么在画幅度-频率图时，对变换后的取模结果直流部分要除以采样数据长度L，交流部分要除以L/2，而频率轴部分要除以FFT变换点数N，为啥要做这么奇怪的处理， 要问答这些问题，还得从傅里叶级数说起。
在做解说之前，先交代一下：FFT是DFT的快速算法（废话），做DFT变换的序列x[n]隐含着周期性，是对x[n]做了周期延拓的，其本质与周期信号的傅里叶级数是一样的。对一L点实序列x[n]，n=1,2,3….L，做N点FFT产生一个N点复数X[k]，k=1,2,3…..N，这里的讨论都是基于实序列的。
下面先说FFT结果对应频点幅度与原始信号的相应频点幅度关系。
傅里叶前辈曾告诉我们：任何一个周期信号可以用一系列具有谐波关系的复指数信号之和来表示（详参文献[1]），它的周期就是基波周期，这就是傅里叶级数的数学理论，其数学表达式如下
（1）
其中是周期信号的基波周期，其谐波频率为，同时假定它是带限信号，否则没法由采样数据重构原始信号了，我们对这个表达式做个变形，重写如下
（2）
由于傅式谱图关于纵轴对称，也即=（实际是共轭对称，即=），对上式的后两项利用逆欧拉公式合并，得到
（3）
对上式中由傅里叶级数积分确定，即
（4）
对上式积分用n（为采样周期，它应满足采样定理）代替t，代替dt，L代替（L为用采样一个周期采样得到的采样点数），积分换成求和，得到上式的一个逼近，其精度与有关
（5）
用采样值x（n）代替采样时刻值，化简得到
（6）
对于上式去掉是不是感觉很熟悉了呢，它就是我们常说的DFT，但还不完全是，注意这里k没有限制。对上式换种写法
=X（k） （7）
将其带入（3）得
（8）
（8）式可看成是x(t)幅度频率的时域表示，（6）是频域表示。（8）中的是其时域对应频点的幅度，用（k非零）表示，得到频域幅度与其时域幅度关系：=，变形得，对于X（0），由（6）和（8）得，到此FFT结果对应频点幅度与原始信号的相应频点幅度关系得到解释。
再看频率对应关系，我们对（5）式做个变形，重写如下

式中= ，（4）式的近似是基于抽样数据的，它的频谱X（jw）是离散且周期性的，即有无限多的频率成分，因此要求限制我们的计算在一个频带宽度=内，通常选择这个区间是 ，这与FFT是一致的，那么（6）式变化成
，k=0,1,2….L-1 (9)
去掉后就是我们真正的DFT了。通过DFT计算出来的结果中第一个点X（0）表示直流分量，第L个点X（L）表示抽样频率，中间被L-1个点平分，每点的频率依次增加，如第k个频点频率即是，其中又被称为物理分辨率，它表示DFT相邻两点间隔所代表的频率宽度，例如L=10，Ts=0.1表示DFT相邻两点间隔1Hz，这也说明如果两个频点的距离小于1Hz，在DFT的结果中是显示不出来的，要想显示出来必须通过增加采样数据L的长度，使小于1Hz，也即增加L可以提高频率分辨率。而FFT为了加快计算速度，采用了基2算法，要求计算点数N必须是2的某次方，对于采样数据长度不足N或大于N分别采用末尾补零和截尾的办法，当采样数据长度不足N，又引出了一个计算分辨率F=，它表示FFT相邻两点间隔频率宽度，补零相当于增加了频域采样点数，有关计算分辨率和物理分辨率的讨论，请参文献[2]。而FFT中横轴的刻度所表示的频率正是基于计算分辨率的，第n点的频率代表着实际频率:
==，记作公式（10），
至此频率对应关系也说完了，还有一点是FFT求出来的谱图关于1/2对称，这是由于抽样造成的频域的周期化的结果，可以把这部分看成相对于纵轴与0<Fk<1/2Fs对称的负频率，具体参看文献[1]。到此所有解说结束。
说了那么多都不如一个具体的例子体会的深切。下面以一个简单的例子说明上述过程。
数学表达式如下
x(t)=1+2*cos(2*pi*50*t)+3*cos(2*pi*100*t)
其中直流幅度为1，2个交流频率：50赫兹和100赫兹的幅度分别为2和3。采样频率Fs=256Hz , 采样点数为256点。256点FFT变换的结果如下图.

先看横轴：计算分辨率F==Fs/N=256/256=1Hz,按照公式（10）=得到第50点所代表的实际频率为50Hz，同理第100点表示100Hz ，结果与测试原信号频点一致。
再看纵轴：从图中看到，在第1点、51点、101点、157点、207点附近有较大值，因后两点是负频率，我们只取出前三点来细看
y =

1.0e+002 *
第1点：2.5600
第51点：2.5600 + 0.0000i
第101点：3.8400 + 0.0000i
计算他们的模值得到，第1点、51点、101点、分别是256、256、384，按照上面的推导结果可以计算出直流幅度A0=256/L=256/256=1, 50Hz幅值A1=256/L*2=2, 同理100Hz幅值为3，可见FFT分析出来的结果是正确的。
为了看到补零对分析结果的影响，对上面的例子做512点FFT变换结果如下图

可以看到纵轴幅度跟上一幅图比没有变化，这也说明了为啥求幅度要用L不用N的原因，也即零值对幅度无贡献。而横轴频率点位置加倍，是因为计算分辨率F==Fs/N=256/512=0.5Hz减半，相邻两点的间隔所代表的频率减半，利用公式（10）让可算出100点和200点对应的实际频率分别是50Hz和100Hz。同时看出补一倍的零使频域的采样点数也增加了一倍。至于图中出现毛边是由于数据截断造成频谱泄露引起的，详情参考文献[1，2]，最后一根突起可看成X[-1]频率。
参考文献
【1】 James H.McClellan等著，周利清等译，信号处理引论 ，电子工业出版。
【2】 王艳芬等著，数字信号处理原理及实现。
附录 Clear;
Fs=256; %采样频率
N=256; %采样点数
t=[0:1/Fs:(N-1)/Fs]; %采样时刻
x=1+2*cos(2*pi*50*t)+3*cos(2*pi*100*t); %原始信号
y = fft(x,256) %做FFT变换
plot(x); %显示原始信号
title('原始信号');
figure;
yy = abs(y); %取模
plot(yy); %显示FFT幅值结果
title('FFT原始幅值');