Hanning窗在插值FFT算法中应用的研究

Hanning窗在插值FFT算法中应用的研究丁文浩摘要：用FFT进行谐波分析时，频谱的泄漏和栅栏效应使得结果（在信号非同步采样时）几乎没有意义。施加旁瓣小的窗可以减小泄漏。在泄漏的影响可以忽略的前提下，可以通过插值克服栅栏效应。本文1）推导了Hanning窗频谱的精确表达式，指出周期性Hanning窗（在主瓣下）具有严格的线性相位；2）在此基础上，得到近似表达式，指出即使在截断长度N不太大时（比如16N），近似表达式也非常逼真；3）在关心的频率范围内，对近似表达式进行Taylor级数展开，得到近似多项式表达式，大大减小了运算量，同时避免了0/0的情况；4）给出了频率、幅值、相位的插值公式。仿真结果显示，这一算法基本消除了非同步采样导致的方法误差。关键词：Hanning窗；插值FFT；谐波分析；多项式StudyonHanningWindowEmployedinInterpolatedFFTalgorithmDINGWen-haoAbstract:WhenusingFFTforharmonicanalysis,spectrumleakageandpicket-fenceeffectmaketheresultsalmostmeaninglessduringnon-synchronoussampling.Byemployingwindowswithsmallside-lobes,theleakagecanbelargelyreduced.Withthepreconditionthattheeffectofleakageisnegligible,picket-fenceeffectcanbeovercomebyinterpolatingthediscretespectrum.Thispaper1)derivestheexactexpressionofthespectrumofHanningwindow,whichrevealsthatHanningwindowhasstrictlinearphaseunderthemain-lobe;2)further,anapproximateexpressionisobtained,whichissopreciseevenwhenNissmall(e.g.,N=16);3)withinthefrequencyrangeconcerned,viatheTaylorseriesexpansion,getsanapproximatepolynomialexpression,whichfeaturesmuchlesscalculationandno0/0condition;4)providesinterpolationformulasforfrequency,amplitude,andphase.Simulationresultshowsthatthealgorithmessentiallyeliminatesthemethoderrorcausedbynon-synchronoussampling.Keywords:Hanningwindow;interpolatedFFT;harmonicanalysis;polynomial0.引言用FFT进行谐波分析时，频谱的泄漏和栅栏效应使得结果（在信号非同步采样时）几乎没有意义。施加旁瓣小的窗可以减小泄漏。在泄漏的影响可以忽略的前提下，可以通过插值克服栅栏效应[1~4]。Hanning窗的第一旁瓣峰值为-32dB，然后以31/或-18dB/OCT的速度下降[5]。因此，可以通过增加截断长度内信号的周期数来显著地减小邻近频率（包括负频率）的旁瓣的叠加影响。当此影响可以忽略时，信号频谱上对应于各个频率成分的主瓣形状将近似于窗函数频谱的主瓣形状。如果此形状可以用某个函数来描述，则可以利用主峰两侧的谱线（通过插值）计算出主峰的高度，从而克服栅栏效应。文章将给出完整的插值算法，并通过仿真验证其精度。1.公式推导1.1.Hanning窗频谱的精确表达式周期性Hanning窗的定义如下：1,,1,0,2cos5.05.0)(NnnNnw（1）其连续频谱为22()0.5()0.25()(),[,)WDDDNN（2）其中，2/)1()2/sin()2/sin()(NjeND（3）是矩形窗的连续频谱，称为Dirichlet核。由于(/2)/22sin(/2)()sin(/2)jjNNNDeeNN（4）因此，(/2)(/2)/2/2/222()0.5sin(/2)0.5sin(/2)sin(/2)sin(/2)cos(/2)0.5sin(/2)sin(/2)sin(/2)1sin(/)jjjNNjNjNeeeWNeNNNeN/2(),[,)jNrWe（5）其中，22sin(/2)cos(/2)()0.5sin(/2)sin(/2)1sin(/)rNWN（6）这就是Hanning窗频谱的精确表达式。注意到，在主瓣下，()0rW，因此，)()(rWW（7）2/)(NW（8）可见，Hanning窗频谱（在主瓣下）具有严格的线性相位。而同样是余弦窗的Hamming窗却不具有线性相位，这使得其频谱表达式难以简化，从而得不到实用的插值公式。余弦窗的一般表达式为[1]02()(1)cos,0,1,,1KkKkkwnaknnNN（9）系数满足约束条件[5]01Kkka（10）当01a时，0()wn为矩形窗；当010.5aa时，1()wn为Hanning窗；当010.54,0.46aa时，1()wn为Hamming窗；当0120.42,0.50,0.08aaa时，2()wn为Blackman窗。容易证明，除矩形窗外，相位线性的充分条件是0(1)0Kkkka（11）可见，Blackman窗也具有线性相位。不过，其插值公式已相当复杂[1]。1.2.Hanning窗频谱的近似表达式为简单起见，进行变量替换N/2（12）得到()(),[/2,/2)jrWWeNN（13）其中，22sin()cos(/)()0.5sin(/)sin(/)1sin(/)rNWNNN（14）此处，我们关心的频率范围是11。当N比较大时，利用sin(/)/NN（15.a）cos(/)1N（15.b）sin(/)/NN（15.c）得到2sin()1()()21rrNWW（16）这就是Hanning窗频谱的近似表达式。令人惊奇的是，即使在N不太大时（比如16N），近似表达式（甚至在全频域）都非常逼真，如图1。这跟近似各项误差的互相抵消作用有关。不过很遗憾，本文作者未能给出严格的数学证明。图116N时，近似表达式与精确表达式的比较Fig.1Comparisonbetweenapproximateandexactexpression(whenN=16)图2是近似表达式的误差曲线。从图2（b）可以得到一个重要结论是：即使在N不太大时，插值算法也能达到相当高的精度。比如，当信号的频率接近Nyquist频率时，这时如果按照消除频谱泄漏影响的原则选取截断长度内信号的周期数，则截断长度N不会很大，此时，同样可以用插值算法来克服栅栏效应。仿真部分将用一个检测电网载波信号的例子来验证这种情况下的应用。（a）全频域（b）插值算法关心的频域(a)wholefrequencyrange(b)frequencyrangeconcernedbythealgorithm图2近似表达式的误差曲线Fig.2Errorcurvesoftheapproximateexpression1.3.Hanning窗频谱（主峰附近）的近似多项式表达式在0附近，对式（16）进行Taylor级数展开（以8阶为例）：222202468ˆ()()rrWWaaaaa（17）其中，021203!41425!61647!81869!10.6449340668482260.1668083584351270.02394346568695690.00220438213069792aaaaaaaaa（18）这就是Hanning窗频谱（主峰附近）的近似多项式表达式。注意到，式（17）的运算量只有5次乘法和4次加法。从后面插值公式的推导可知，应用近似多项式表达式时，关心的频率范围可进一步缩小为0.50.5。在此频域内，误差仍能保持在-95dB以下，如图3。图3近似多项式表达式的误差曲线Fig.3Errorcurvesofapproximatepolynomialexpression1.4.插值公式考虑单一频率信号(2)()xxjFtxxtAe（19）其N点采样序列(2/)()()(/)e,0,1,,1xsxxxjFFnsxjnxxnxnFAeAnN（20）其中，sF为采样频率；2/xxsFF（21）为归一化角频率（或称为频域角频率）。()xn加窗后的连续频谱为()(),[,)xjwxxXAWe（22）或者，按照式（12）进行变量替换，为()(),[/2,/2)xjwxxXAWeNN（23）其中，Nxx/2（24）取round()xxk，即最接近x的整数，有()()()()xxxxjwxxxxjxxjjxrxXkAWkeAWeAWee（25）其中，round()xxxxxk（26）显然，5.0x，因此()()wxxrxXkAW（27）()wxxxXk（28）于是，解得幅值和相位的插值公式：11()()ˆ()()xwxwxrxrxAXkXkWW（29）()xwxxXk（30）注意到，()wXk是加窗后信号的离散频谱，可以由FFT得到；而()wxXk是其中幅值最大的那条谱线。频率的插值公式可以由式（26）、（24）、（21）推得：()sxxxFFkN（31）剩下的关键问题是如何得到x。这里将根据幅值次最大与最大两条谱线的比值，通过窗函数频谱的近似表达式来估算。分四舍、五入两种情况：ⅰ）当1xxxkk时，定义(1)()wxwxXkXk，则(1)12()rxxxrxWW（32）解得211x（33）ⅱ）当1xxxkk时，定义(1)()wxwxXkXk，则(1)12()rxxxrxWW（34）解得211x（35）2.仿真验证例1．电网谐波分析借用文献[2]中的输入信号911()cos(2)mmmxtAFt（36）其中，基波频率149.5FHz；其它参数的设定值见表1。取2500sFHz，1024N。这样，截断长度内包含略大于20个（基波）周期的信号。仿真结果见表1。表1例1的仿真结果Tab.1SimulationresultofExample1谐波次数m幅值mA相位m设定值（V）测量值（V）相对误差（%）设定值（°）测量值（°）绝对误差（°）1310310.0005e-555.0010.00121.551.54418-0.3761010.2970.29739.309.300400.0042020.0090.00941.241.240180.015120119.965-0.035515.515.4998-0.001150150.0010.0016