基于MATLAB的FFT离散频谱分析

周盼　甘丽群　刘华超

摘要：通过仿真形象揭示了MATLAB上FFT过程中频谱泄漏和栅栏效应现象，重点从理论上分析了其产生原因。针对以上两情况导致的频率和幅值上的误差，介绍了常用的校正方法的原理、推导过程以及运用的局限性。指出了幅值校正过程中所有恢复系数全为[2N]的常见错误。为MATLAB上FFT的理解和应用提供一定的帮助作用。

关键词：MATLAB；FFT；频谱泄漏；栅栏效应；频率幅值校正

中图分类号：TP391 文献标识码：A 文章编号：1009-3044（2016）36-0277-04

Analysis of Discrete Spectrum Based on FFT using MATLAB

ZHOU Pan1， GAN Li-qun2 ， LIU Hua-chao2

（1. School of Electronic Engineering ， Xi'an Shiyou University， Xi'an 710065，China； 2. School of electrical and Information Engineering，Chongqing University of Science and Technology，Chongqing 401331，China）

Abstracts： The phenomenon of spectrum leakage and packet fence effect in the process of fast Fourier transform was revealed through simulation， and the cause of which was analyzed emphatically in theory. Aiming at errors of frequency and amplitude caused by FFT in MATLAB， the principle， derivation and the application limitations of the common correction methods was introduced in detail. The frequent fault which regarded the all recovery coefficients for amplitude correction as [2N] was also pointed out. The paper can provide some help for the understanding and application of FFT.

Key words： MATLAB； FFT； spectrum leakage； fence effect； frequency amplitude correction

信號处理广泛运用于语音图像处理、通信、生物医学等领域，是一门极其重要的学科。随着快速傅里叶变换的诞生，更加倾向于从频域分析采集到信号的信息。由于时域中截断导致的频谱泄漏以及频域的离散化产生的栅栏效应，若采用FFT对信号进行检测，不可避免地存在频率、幅值和相位的误差。然而很多工程实际应用如电力谐波检测、振动信号分析中，需要较准确地测量出信号的频率和幅值。因此深入了解频谱泄露和栅栏效应的产生原因，以及如何对FFT后频率和幅值进行校正具有很大的科研和工程意义。MATLAB作为一款应用方便、功能强大的仿真软件，一直受到广大科研工作者的运用。本文围绕MATLAB上FFT中频谱泄漏和栅栏效应进行了较为全面的介绍，对与其相关的常见问题做了深入分析。

1 频谱泄漏

设周期信号为[x（n）]，其表达式为：

[x（nΔT）=A0cos（2πf0nΔT+φ0）] （1）

其中[A0]、[f0]、[φ0]、[ΔT]分别为幅值、频率、相位、采样频率，且[ΔT=1fs]。对信号频率进行归一化处理令：

[f0=λ0fsN]，[Δf=fsN]为频率分辨率，则信号可表示为：

[x（n）=A0cos（2πλ0Nn+φ0）， n=0，1，2，...，N-1] （2）

周期信号为无限长，计算机无法对其进行傅里叶变换，因此需对采集到的信号作截断处理。对于MATLAB上的FFT对信号的截断相当于加矩形窗。截断后的信号表示为[xw（n）=w（n）x（n）]，[w（n）]为矩形窗。由卷积定理可知[x（n）]的傅里叶变换的幅值为：

[Xw（k）=W（k）*X（k）=A02ejφ0W（k-λ0）+A02e-jφ0W（k+λ0）]（3）

其中[W（k-λ0）]、[W（k+λ0）]分别为傅里叶变换后的正负频率成分。当[5<λ0

[Xw（k）?A02W（k-λ0）] （4）

图1为FFT变换过程中的频谱泄漏示意图。（a）为原始周期信号的时域图，理论上其频谱图上只有在其真实频率[λ0]处才有能量，为一个脉冲信号，如（b）所示；截断函数的时域、频域如（c）（d）所示；截断后的信号如（e）所示，信号的长度由无限长变成了时间T以内。由卷积定理可知，时域的相乘对应频域的卷积，则截断后信号的频谱（f）为（a）（d）的卷积。而任何函数[y（n）]与脉冲函数的卷积，就是脉冲函数处当做坐标原点对[y（n）]进行重构[2]。图（f）即为图（d）在图（a）处的重构，值得注意的是，重构过程中需要把图（d）中所有点的幅值乘以脉冲函数的幅值。通过以上分析可知，截断导致信号的频率泄漏到整个频段，且频谱图为矩形窗的形状。由于不同的窗函数具有不同的主瓣宽度和旁瓣衰减率，实际应用中可根据需要对信号加不同的窗函数，以抑制FFT变换过程中的频谱泄露[3-5]。

2 栅栏效应

对于MATLAB中FFT，式（4）中[k]只能取整数，[k]表示的频率为[kΔf]。信号的实际频率为[f=λfsN]，其中：[λ=l+δ]，[l]表示整数，[δ]表示小数。当信号的实际频率为[Δf]的整数倍，即整周期采样时，[δ]为0。[λ=5]的FFT后的仿真频谱如图2所示，谱线最高点对应的频率为5，表现出了真实频率。然而实际工程应用很难达到整周期采样，即[δ]不为0。[λ=4.8]的FFT后的仿真频谱如图3所示，谱线最高点对应的频率为5，未能正确表示出信号的真实频率，此时的FFT变换就带来了误差。两图相比较可看出，信号整周期采样时，除了最高谱线有幅值外，其他谱线幅值均为0；非整周期采样时，所有谱线幅值均不为0。此现象可从理论上解释如下。矩形窗的离散傅里叶变换[6]为：

[W（k）=sin（kπ）sin（kπ/N）] （5）

当[k≤0.5]时，[k?N]，[kπ/N→0]，[sin（kπ/N）?kπ/N]，则式（5）可化为：

[W（k）=Nsin（kπ）kπ] （6）

将（6）代入（4）中，可得FFT变换后信号的幅值为：

[Xw（k）?A0N2sin（（k-λ0）π）（k-λ0）π] （7）

因整周期采样时[λ0]为整数，而[k]也为整数，则[k-λ0]全为整数。当[k=λ0]时，利用洛必达法则可得：

[sin（k-λ0）π（k-λ0）π?]1 （8）

此时

[Xw（λ0）?A0N2] （9）

即为信号的真实幅值。当[k≠λ0]时，[sin（k-λ0）π≡0]

此时

[Xw（k）≡0] （10）

非整周期采样时[λ0]为非整数，而[k]为整数，则[k-λ0]全为非整数。此时，对于所有[k]，[Xw（k）]全为非0，且[k]不可能等于[λ0]，频谱图也表现不出信号的真实频率和幅值。栅栏效应主要是由于频域的离散化造成的，可采取补零的措施减小频域抽样的间距，达到FFT变换后表现的最高谱线更加接近真实谱线的效果。

3 FFT中频谱泄露和栅栏效应产生误差的校正

工程实际中，由于信号的频率未知，很难达到整周期采样，利用FFT将会导致最大为[0.5Δf]频率误差，进而产生幅值误差。而在很多工程领域需要较精确地检测出信号的频率和幅值，因此对FFT变换后的结果进行校正显得非常重要。

3.1 频率误差校正

频率的误差来源于时域的加窗和频域的离散化，早期就从原理上发现了校正方法[7-8]。为了分析的简单性，现介绍矩形窗下的校正方法。设最高幅值谱线频率为[l]，当次高谱线为[l+1]时，由于实际的最高谱线位于FFT变换后的最高谱线与次高谱线之间，此时真实频率可表示为[λ0=l+δ]，将其代入（7）可得：

[Xw（k）?A0N2sin（（k-l-δ）π）（k-l-δ）π] （11）

现定义次高谱线与最高谱线幅值比为：

[α=Xw（l+1）Xw（l）] （12）

将（11）代入（12）得：

[α=sin（（l+1-l-δ）π）（l+1-l-δ）π.（l-l-δ）πsin（（l-l-δ）π）=δ1-δ] （13）

因此：

[δ=α1+α] （14）

当次高谱线为[l-1]时，此时真实频率可表示为[λ0=l-δ]，同理可求出：

[δ=-α1+α] （15）

最终，频率的校正量可表示为：

[λ0=l±α1+α] （16）

当次高谱线位于最高谱线右侧时取‘—，反之，取‘+。

3.2 幅值误差校正

若信号为整周期采样，由（9）可知，将表现出的最大值乘以[2N]的系数即可得到信号的真实幅值。对[x（n）=cos（2π*64nN+φ0）]，[N=256]的信号作FFT，其频谱如图4所示，可看出信号的幅值为128，乘以恢复系数可得真实幅值为128*2/256=1。当信号为非整周期采样时，不能像有关文献提出的直接乘以[2N]的系数恢复真实幅值。对[x（n）=cos（2π*64.5*nN+φ0）]，[N=256]的信号作FFT，如图5所示。若也按[2N]恢复真实幅值，则为81.48*2/256=0.637，误差达到0.37。现对FFT变换后如何恢复真实幅值做一下理论解释。

整周期采样时，频域抽样时肯定会抽到真实频率[λ0]，将[k=λ0]代入（7）可得出FFT变换后表现出的最大值为[A0N2]，因此乘以[2N]的恢復系数，即可得到真实幅值[A0]。非整周期采样时，频域抽样根本抽不到真实频率[λ0]，表现出来的最大幅值为最靠近[λ0]的谱线，设最大幅值对应的频率为[l]，幅值为[Y]。利用（7），那么此时的真实幅值可校正为

[YR=Y2N（l-λ0）πsin（（l-λ0）π）] （17）

式中的[λ0]为信号的真实频率，可通过3.1节的频率校正先得到[λ0]的近似值。当信号为最大非整周期采样时，如图所示，[l-λ0=0.5]，此时[（l-λ0）πsin（（l-λ0）π）=1.57]，用（17）可校正出真实幅值为81.48*2/256*1.57=0.99，就基本上接近真实幅值了。

4 总结

离散傅里叶变换过程中频谱泄漏和栅栏效应的充分理解，可以帮助更好地学习信号处理中的其他知识。目前频率和幅值的校正理论都是建立在负频率干扰被忽略前提下的，因此对于运用到密集频谱的校正中，将会具有很大的误差。采集到信号中难免掺杂有噪声，可把频谱校正理论的抗干扰性当做重点研究内容。

参考文献：

[1] Luo J， Xie Z， Xie M. Interpolated DFT algorithms with zero padding for classic windows[J].Mechanical Systems & Signal Processing，2016，s70–71：1011-1025.

[2] 丁康，谢明，杨志坚 - 科学出版社[M] ， 2008，13-14.

[3] 祁才君，王小海. 基于插值FFT算法的间谐波参数估计[J]. 电工技术学报，2003，18（1）：92-95.

[4] 黄纯，江亚群. 谐波分析的加窗插值改进算法[J]. 中国电机工程学报，2005，25（15）：26-32.

[5] 钱昊，赵荣祥. 基于插值FFT算法的间谐波分析[J]. 中国电机工程学报，2005，25（21）：87-91.

[6] 谢明，丁康. 频谱分析的校正方法[J]. 振动工程学报，1994（2）：172-179.

[7] Jain V K， Collins W L， Davis D C. High-Accuracy Analog Measurements via Interpolated FFT[J]. Instrumentation & Measurement IEEE Transactions on，1979， 28（2）：113-122.

[8] 潘文，钱俞寿，周鹗. 基于加窗插值FFT的电力谐波测量理论—（Ⅰ）窗函数研究[J]. 电工技术学报，1994（1）：50-54.