加窗FFT在频谱分析中的应用

1、加窗FFT在频谱分析中的应用1离散傅里叶变换(DFT)1.1DFS板周期为N序列对应基频N，各次谐波序列ekn=瞠卜0J-N -1 n 展开为 n = IDFS Xk = -1 Xk e，5 Nk =0由正交性祝e，艾蚌n=0(.2兀) eJNrnI J得：Xk = DFS n=Nne-j寸knn=0*1N k = r0 others=N 切 8 k -1 - mNm=s1.2DFT考虑有限长序列xn,长度为N,定义周期序列 n=云x n + rN = x(n)Nr=-s和xn = x(n) R n通过取 n和Xk 的一个周期，定义DFT如下：X k = DFT xn)= 2 xnW kn0

2、k N -1N n=0 xn = IDFTXk= 2 XkW-kn 0 n N -1 NNk=0对于长度为N的有限长序列，由DFT可重构DTFT，2-1X (k帅一k=0N sin 其中中Co)=一e-j 2 2 -N sin 21.3FFT提高DFT运算效率的途径:1、利用W心的性质：周期性、可约性、共轭对称性N2、把长序列化解成短序列进行计算，通过迭代运算来减少DFT的运算量时间抽取FFT算法简介如下：X(0)X(1)X(3)oX(|)oX(5)0X(6)冲1、对时域n序号进行奇偶分解：将二进制表述的数位码按高低位码进行倒置得到，以n=8为例：正序：000,001,010,011,100,

3、101,110,111 艮口 0,1,2,3,4,5,6,7倒序：000,100,010,110,001,101,011,111 艮口 0,4,2,6,1,5,3,72、第一次分解：做二个N/2点的DFT第二次分解：做四个N/4点的DFT做N/2个两点DFT使用迭代算法3、本地即可实现，不需要额外申请空间。2利用DFT做频谱分析过程:米样 (2)加窗 (3)DFT (或需补零) (4)插值得到连续频谱2.1米样与混叠x (t),傅立叶变换为 X (jQ),xn = x (nT),DTFT为 X(ej)X(ej) = T * X (j t + j )，Q和3之间有关系为：w = QT 或Q =

4、/T。r=-s防混叠滤波2.2分辨率与频谱重构实际分析的信号只能是有限长信号，必须进行加窗处理，从而导致分辨率有 限，如同从百合窗中观察频谱，只能看到特定频点的幅频信息，此时需要对频谱 进行重构，从而得到完整频谱，根据奈奎斯特采样率，这种重构是可以无失真的 实现的。xn = A cos n + A cos n x n =(A cos n + A cos n)R n001110011 NX侦)=A。点k=-81116 (2冗k) + 6 (+-2冗k)+ A兀另1k=-s6 (一-2兀k) +1+ 6 (+-2冗k)v ( ) A L () w () A L () w (*Xe j= We jo

5、+o+ We 浪-w+1We 迥+3+ W e /o-o12 RR2 R 1 R 1两个余弦函数和加矩形窗后的频谱图2.3泄露与栅栏现象在对应X。为零的点，x (k证零1泄漏现象：对x1。奏样，得至叫(k),由于加窗作用，栅栏现象：Mx 1。墨样，x 1 (k)可能不在乂。的极值点，降低谐波的频率估计精度， 通过补零，可以有效弥补栅栏现象，但补零不能提高频谱分析的分辨率。2.4窗函数及加窗0 n M otherwise2n/M0 n M/2Bartlett窗：wn日2-2n/MM/2 n M0 otherwise鼻0.5 0.5cos2兀n/M0 n MHanning 窗： wn=10 oth

6、erwiseHamming 窗:10.54 - 0.46 cos 2 / M wn= 10 otherwise0.42 - 0.5cos2m / M + 0.08cos4m / M0 n MBlackman窗:wn=10 otherwise窗函数性能分析如下:窗函数名称卜主瓣宽度Q3dB带宽一最大旁瓣电平q矩形卜Nq0.89 *2诚 Nq-13 dB#Bartlett1E N*1.28*2 诚 Nq-25dBHanning.WWWWW.-Vr-E N#1.44*2/E-31dBHamming 】1.3*2招-41dBBlackman1.68*2 诚 Nq-57dB3窗函数选择与频谱重构3.1窗

7、函数的一般形式窗函数的一般功能为截断无限长信号为有限长信号，但性能较好的窗函数还 应当有良好的频域特征，但由于测不准原理，窗函数时域受限的同时不可能在频 域上得到单频点脉冲，必定也是频谱受限信号，由于余弦函数良好的频域特征， 一般情况下，我们基于余弦函数构造需要的窗函数。fa + a cos2兀n/M + a cos4兀n/M + a cos6兀n/M0 n Mwn = 01230 otherwise从复杂度出发，我们只考虑三阶以内的窗函数。3.2加窗DFT的重构算法3.2.1相关理论计算I a + a cos 2兀n / M + a cos 4兀n / M + a cos 6兀n / M0

8、n l X(k) 1- k o o Z-j a,(p(A) (p(0 ) +07im (-0 )2 -m20 Af +1m=00令k = X om X (k) 1= ke，a 7im 人2 m2m=0A sin(7i0 )k (-14t7 -13a -10a -5a 认4+(49。+ 36a +9a + 4a )X2 -36a=k00123012307112 (X2 1)(人2 4)(人2 -9)其中，令i + a + a + a = 0 01233.2.2参数选择原则1、简化插值重构算法原则为简化插值重构过程，我们令上式中高次项系数均为o,得到三阶窗函数参数选择公式如下： TOC o 1-5

9、 h z a + a + a + a =0,线性相位 0123a -a +a -a =1,线性相位012314。一13。一10。一 5。= 0 012349。+36。+9a +4。= 0V 01232、旁瓣最快衰减原则旁瓣衰减速度由窗函数的光滑程度有关，对窗函数求导，由窗函数导数的连 续性，可以得到该窗函数的最大可导阶次，此时，窗函数参数满足a + a + a + a = 0,线性相位a0-a11+ a22-a33= 1,线性相位a + 4a + 9a = 0a + 16a + 81a = 0我们的计算过程显示，最大旁瓣衰减速度的余弦窗函数正是我们需要的具有 最简插值重构算法的余弦窗函数，因此

10、，从理论上讲，上述计算得到的窗函数是 所有三阶窗函数中旁瓣衰减最快的。继续推导发现，更高阶窗函数同样能够同时满足以上两个原则。3、窗函数阶数选择原则窗函数阶数越高，旁瓣衰减越快，但同时主瓣宽度越宽，对于FFT来讲，可 能淹没高次谐波，因此，选取窗函数时阶数要适中。3.2.3三阶窗函数原则上采用数值分析算法，任意的窗函数都可以计算出频偏人，进而求出幅 度信息Ak，但为了算法实现的方便，我们定义自己的窗函数满足以下条件：a + a + a + a = 0,线性相位a0-a1+ a2-a 3=1，线性相位 0123-14a - 13a - 10a - 5a = 049a + 36a + 9a + 4

11、a = 0n a0 = 0.3125，a1 = -0.46875，a2 = 0.1875，a3 = -0.03125此时，DFT重构算法设计如下:I X(k +1)1I X(k)I1、 kA k4以-32、人= kk a +13、4、,兀X (k)X (1人 2)(4-人 2)(9-人 2)A = kkkkk11.25 sin(兀人)XkM中=中(k)- Z 1m +1FFT的结果是对连续频谱进行采样，因此需要合理设计频率分辨力，从而得 到窗函数的完整频谱。由于三阶窗函数的主瓣宽度为4个频率分辨单位，旁瓣宽度均为1个频率分 辨单位，因此Matlab仿真时取采样率为2048，计算点数为4096时

12、，可准确显 示出旁瓣电平。Matlab仿真得到窗函数形状如下:由上图可知，该三阶窗主瓣宽度为4个分辨单位，对于50Hz的电网基波来 讲，当分辨率为12.5Hz时，该窗函数刚好可以保证各次谐波之间不产生混叠； 最大旁瓣为-62.8dB。3.2.4二阶窗函数对二阶窗的情况有I x（k） | sin（兀 ）X （-6a - 4a - a ）X2 + 4ak 兀 0心（入2 -1）（X2 - 4） 0a +，2=0,线性相位a0-ai1+ a22 = 1，线性相位-6a - 4a - a = 01012n a0 = 0.375，a1 = -0.5，a2 = 0.125此时，DFT重构算法设计如下:1、

13、I X (k +1)1I X (k)12、3、4、_ 3o-2k 以+1k, 冗X (k)X (1 人 2)(4 人 2) A =kkkk1.5sin(KXk)Matlab仿真得到窗函数形状如下:Matlab仿真blackman窗函数形状如下:比较两图可知，blackman窗函数在对第一旁瓣的抑制能力上，要比我们设 计的二阶窗函数优异，但从第二旁瓣开始，我们设计的二阶窗函数快速衰减，远 比blackman窗函数性能更好，仿真也验证了该窗函数的最快旁瓣衰减特性。而且从重构算法的实现策略上，我们设计的二阶窗函数更为简单实用。3.2.5 一阶窗函数对一阶窗的情况有I X(k)1= 4 sin叫)七

14、J兀Xi(ki -1)Ja +a】=0,线性相位|a0-ai1 = 1，线性相位n a = 0.5，a = -0.5，此即为Hanning窗此时，DFT重构算法设计如下：1、I X (k +1)1I X (k)12、3、4、中=中(k)-Matlab仿真得到窗函数形状如下:从结果上来看，一阶窗函数中，旁瓣衰减最快的正是hanning窗2以-1，=kk 以+1A = 2kX (k )气(1 _气) ksin(兀人)3.2.6更高阶窗函数对于更高阶的窗函数，其旁瓣衰减更快，但主瓣宽度更宽，对于谐波分析而 言，可能会导致基波被拉伸而淹没高次谐波，因此并不适用，此处只作简要计算 分析(以5阶为例)。I

15、 X(k)I=气，诅叫知,兀m入2 - m 2m=0A sin(K0 )Xb 入8 + b 入6 + b 入4 + b 入2 -14400。= n入2 (入2 -1)(入2 - 4)(X2 - 9)(X2 -16)(入2 - 25)a + a + a + a + a + a = 0,线性相位 012345a - a + a - a + a - a = 1,线性相位 012345b =55a -54a -51a -46a -39a -30a = 0b = 1023a + 969a + 819a + 609a + 399a + 273a = 0b = -7645a - -6686a -4369a

16、-2164a -1261a -820a = 0b = 21076a + 14400a + 3600a +1600a + 900a + 576a = 0n a = 0.24609375, a = -0.41015625, a = 0.234375, a = -0.087890625,a4 = 0.01953125, a5 = -0.001953125此时，DFT重构算法设计如下：1、I X (k +1)1I X (k)12、3、4、6以-5k以+12兀X(k)X (1 -X 2)(4-X 2)(9-X 2)(16-X 2)(25-X 2) k kkkkk3543.75 sin(兀Xk)Matla

17、b仿真得到窗函数形状如下:0-20-40-60-80-100-120-140-160-1800102030405060708090100从性能上来讲，高阶窗函数第一旁瓣更低，而且旁瓣衰减也更快，但主瓣宽 度多达6个频率分辨单位，在频率分辨力较低的系统中，并不适用于进行谐波分析。4窗函数在谐波分析中的应用及性能假设电网单相电压表达式为u G）= 220 Xv；2cos（2&t + 0 ）+ 10cos（4Kf t + 0 ）+ 20cos（6Kf t + 0 ）+ 15cos（8Kf t + 0 ）A01020304其中基波f = 50Hz的幅度、相位分别为A = 220 x克、0 = 30。；

18、二次谐波 0112f0 = 100Hz的幅度、相位分别为A2 = 10，02=10。；三次谐波3f0 = 150Hz的幅度、相位分别为A3 = 20，0广45。；四次谐波4f0 = 200Hz的幅度、相位分别为A4 =15，04=60。当系统采样率为51200Hz，采样点数为4096点时，频率分辨力为12.5Hz，主瓣刚好落在第5个采样点上，高次谐波全部落在旁瓣电压为0处，故FFT过程 并不存在频谱泄露，采用任何窗函数效果都是一样的，也无需进行频谱重构。而实际电网中基波频率可能会出现偏离50Hz的情况，此时FFT过程存在频 谱的泄漏，且FFT结果无法准确反映电网频率，此时需要加窗处理并对频谱进行 重构，才能得到准确的频谱分析结果。现假设基波f = 49.5Hz，分别用矩形窗和第3节计算得到的一阶、二阶、 0三阶窗函数进行频谱分析及重构，得到结果如下：1、矩形窗：f0 = 50Hz谐波次数幅度（V）相位（）计算参数精确参数计算参数精确参数1308.3088311.12723.059930210.6987102.864810319.9582024.919845415.55761532.5669602、一阶窗：f0 = 49-5025Hz谐波次数幅度（V）相位（）计算参数精确参数计算参数精确参数1311.1605311.12729.959430