fft插值算法

性函数 3.2 窗函数的特性与分析
窗类型 矩形窗 Hanning 窗 Hamming 窗 Blackman 窗 Blackmaan-Har ris 窗 FDMS 主瓣宽 度 4π/N 8π/N 8π/N 12π/N 16π/N 16π/N 旁瓣宽 度 4π/N 4π/N 4π/N 4π/N 4π/N 4π/N 旁瓣峰 值（dB） -13 -31 -41 -57 -92 -61 旁瓣衰减速 度 （dB/oct） 6 18 6 18 6 42
电子技术
电子技术 基于 DSP 的 FFT 加窗插值算法研究
王培慧 邹连明 宁波三星电气股份有限公司，浙江 宁波 315034
摘要： 快速傅里叶变换是最主要的电力谐波分析方法， 但是由于很难做到同步采样和整周期截断造成了频谱泄露和栅栏效应， 对精确地测量产生较大的影响。采用加窗以及插值修正算法是抑制频谱泄露和栅栏效应的有效方法。本文针对加窗插值算法 利用改进的四项余弦窗函数及其相应的插值修正公式，计算出实际频率点的幅值、相位。通过在 DSP 平台上仿真结果验证了 该窗函数及其修正公式算法的准确度优于汉宁窗和 Blackman 窗。 关键词：谐波分析；四项余弦窗函数；插值算法 中图分类号： TM714 文献标识码： A 文章编号：1671-5810（2015）62-0156-02 1 导言 随着电力电子设备的在电力系统中的广泛应用，电网谐 波污染日益严重，对谐波的精确分析对谐波的治理有重大意 义。快速傅里叶变换算法因其易于微机实现成为谐波分析的 主要方法。由于电网频率的波动，采样周期非整数倍的情况 下，各次谐波分量不一定正好位于频率分辨点上，而是有可 能位于两个频率分辨点之间，这样直接使用 FFT 算法对谐波 分析有严重的影响，不能准确的计算出电网谐波的准确分量。 随着对精确测量精度的更高要求，对信号施加合适的加窗函 数并用准确的插值公式修正后，信号的采样不再局限于整周 期采样，当信号频率波动造成非整数倍周期采样时，也可以 精确地计算出信号的幅值、相位和频率。 2 频谱泄露和栅栏效应 在从一个连续信号中抽取一个数据块的过程，如果不加 窗， 实际上就是默认加了一个矩形窗， 这样数据抽取的结果， 就是使得原来连续信号中集中在 f0 处的能量被分散到两个 较宽的频带中去了，这种现象称之为频谱能量泄漏。 频谱泄漏是由 FFT 算法中的假设造成的，FFT 算法中假 设离散时间序列可以精确地在整个时域进行周期延拓，而且 所有包含该离散时间序列的信号为周期函数，周期与时间序 列的长度相关。然而，实际计算机处理过程中，一方面不能 做到采样点数无限大，另一方面绝大多数情况下所处理的是 一个未知的平稳信号，不能保证采样点数为信号周期的整数 倍。 频谱泄露从严格意义上讲，可以分为两方面的泄露，长 范围频谱泄露和短范围频谱泄露。长范围频谱泄露是指由于 信号截断造成的信号频谱旁瓣之间的相互干扰，能量由集中 于中间一点扩散至整个频域，同时旁瓣之间相互干扰。短范 围泄露是对信号做 FFT 时，得到的是一系列离散的谱线，如 果信号中的频率成份位于谱线之间而不是正好落在谱线上， 此时就会造成幅值和相位上的偏差，离散的一条条谱线就象 一个栅栏，因此也称这种现象为栅栏效应。 3 窗函数特性及应用 加窗函数是对 FFT 算法进行改进的一种重要手段，选择 合适的窗能够有效抑制频谱泄漏，从而能够提高测量的精度。 加窗函数的实质就是对被分析信号在不同时刻加不同的权 值，以使截断所造成的影响尽可能的小。由于主瓣宽度决定 了被截断以后所得信号序列的频域分辨率，而旁瓣峰值有可 能湮没信号频谱中较小的成分。因此，一个理想的窗函数应 该具有最小的主瓣宽度和最小的第一旁瓣峰值（dB）以及最 大的旁瓣峰值衰减速度（dB/oct） ，从而使频域能量主要集 中在主瓣内， 我们便可得到接近于真实的频谱信号。 主瓣窄， 频率分辨率高，主瓣宽，通带和阻带之间的过渡带宽；第一 旁瓣峰值最小，衰减速度快，频谱泄露小，使得滤波器通带 与阻带之间的波动减小，影响输出信号的频谱分析精度；旁 瓣峰值衰减速度最大同样是为了减少频谱泄露。 但是根据 FFT 算法原理，FFT 产生离散的频谱，出现在 FFT 每个谱线的是在每个谱线上的连续卷积频谱。如果原始 信号的频谱成份与 FFT 中的谱线完全一致，这种情况下采样 数据的长度为信号周期的整数倍，频谱中只有主瓣。没有出 现旁瓣的原因是旁瓣正处在窗函数主瓣两侧采样频率间隔 处的零分量点。如果时间序列的长度不是周期的整数倍，窗 函数的连续频谱将偏离主瓣的中心，频率偏移量对应着信号 频率和 FFT 频率分辨率的差异，这个偏移导致了频谱中出现 旁瓣。所以窗函数的旁瓣特性直接影响着各频谱分量向相邻 频谱的泄漏宽度。实际上，能同时具备以上 3 个指标的理想 窗函数是不存在的，对窗函数进行选择时应满足以下基本要 求：时域为改善截断处的不连续状态（由于吉布斯现象而产 生的振荡） ； 频域为窗谱的主瓣要求窄而高， 以提高分辨率， 旁瓣应小，正负交替接近相等以减小泄漏和假频。 但同时也应注意到，不管是加任何窗函数还是增加采样 长度，都只可能在一定程度上抑频谱泄漏误差和栅栏效应， 将误差减小到可以接受的程度，而不能完全消除由时域截断 和离散化带来的误差。 3.1 窗函数的基本概念 设x （n） 是一个序列， w （n） 是一个长度为 N 的窗函数， 用 w（n）截断 x（n） ，得到 N 点序列 xn（n） ，即

第37卷第1期2003年1月浙　江　大　学　学　报(工学版)Journal of Zhejiang U niversity (Engineering Science )V o l .37N o .1Jan .2003收稿日期:2002204215.作者简介:祁才君(1964-),男,江苏宜兴人,副教授,从事数字信号处理研究.E 2m ail :qcj @zju .edu .cn应用插值FFT 算法精确估计电网谐波参数祁才君,陈隆道,王小海(浙江大学电气工程学院,浙江杭州310027)摘　要:深入研究了插值快速傅里叶变换(FFT )算法在电网谐波参数估计中的应用.加窗宽度和窗函数的类型是影响插值FFT 算法分析精度的主要因素.通过对常用窗函数的特性分析,得出了加窗宽度关于分析精度的估计公式.电网信号的基波幅值远大于各次谐波幅值,分析表明,H anning 窗比较适合分析电网信号,同时给出了基于H anning 窗的电网谐波幅值、频率和相位的显式计算公式.仿真结果证明,应用上述分析结果,电网谐波幅度、频率和相位的估计达到了预期的分析精度.关键词:谐波分析;电力系统;电能质量;快速傅里叶变换中图分类号:TM 933 文献标识码:A 文章编号:10082973X (2003)0120112205H igh -accuracy esti m ation of electr ical harm on ic param eters by usi ng the i n terpolated FFT algor ith mQ I Cai 2jun ,CH EN L ong 2dao ,W AN G X iao 2hai(Colleg e of E lectrical E ng ineering ,Z hej iang U niversity ,H ang z hou 310027,Ch ina )Abstract :A tho rough research on the h igh 2accu racy esti m ati on fo r the electrical har m on ic param eters by u sing the In terpo lated FFT algo rithm w as conducted .T he w indow size and the w indow type are the m ain facto rs influencing the analysis p recisi on .B ased on the featu res of the comm on ly u sed w indow s ,a fo r m u la fo r the analysis p recisi on in relati on to the w indow size is given .Studies show that acco rding to the pow er signal’s characteristic that the fundam en tal am p litude is m uch h igher than that of the har m on ics H ann ing W indow is m o re su itab le fo r the analysis of electrical signals and the com pu ting fo r m u la fo r the har m on ic p aram eter esti m ati on are given as w ell.Si m u lati on resu lt dem on strates that by u sing the resu lts above the esti m ati on of har m on ic frequencies ,am p litudes and p hases reaches the exp ected p recisi on .Key words :har m on ic analysis ;pow er system ;pow er quality ;FFT 随着电子技术和电力电子器件的发展,电力电子设备的应用越来越广泛,这些非线性电力设备对电网产生了很大的污染.同时,以微处理器为核心的各种控制设备、自动化设备对电能质量十分敏感,瞬间的电能恶化,有可能导致程序混乱,从而产生重大的质量事故.显然电网参数的快速、准确监测对电能质量的治理具有十分重要的意义.本文将重点讨论电网谐波(包括基波)参数的精确估计方法.最直接、有效的电网谐波参数的分析方法是应用FFT 算法.但是,应用FFT 算法分析周期信号频谱的前提是要求对信号进行同步采样,即要求采样长度与信号周期成整数倍关系[1].但是,实际电网频率通常总会在额定频率(我国为50H z )附近波动.为了保证同步采样,方法之一是在采集系统中加入锁相同步技术,由硬件实现同步采样[2];另一种方法是首先通过均匀(大部分为非同步)采样,然后由一定的算法解决非同步问题[3～9].前一种方法的优点是信号处理比较简单.但由于锁相环响应速度较慢,不能及时跟踪信号频率的快速变化,从而不能实现真正意义上的同步采样.后一种方法具有一定的通用性,文[3]提出了采样数据的软件同步化算法,该算法首先由过零检测法[9]估计信号基频,然后通过插值对原始数据进行同步化处理,并作FFT分析,算法的精度对基频估计的准确性有很大的依赖性,同时由于插值过程带来的非线性,导致谐波的相位估计误差过大.文[7,8]提出了插值FFT算法,有效地改善了对谐波幅值,尤其是谐波相位的估计精度.文[9]进一步讨论了各种升余弦窗在插值FFT算法中对减少泄漏的作用.本文在详细分析插值FFT算法的基础上,提出了窗函数的选择准则,给出了加窗宽度的估计公式.同时根据电网波形的特点,对谐波参数的计算公式作了修正,得到了更为精确的谐波参数估计.1　插值FFT算法1.1　正弦(余弦)信号非同步采样时的频谱分析 为了分析方便,将简谐信号表示成如下的解析形式:x a(t)=A1・e j(2Πf1t+Υ1).(1)式中:A1为信号幅值;f1为信号频率;Υ1为信号初相.简谐信号经采样并矩形截断后的序列为x(n)=x a(nT s)・w R(n),n=0,1,…,N-1.(2)式中:x a(nT s)为x a(t)的无限长采样序列;T s为采样周期;w R(n)为矩形窗;N为采样点数.由此可得采样频率f s和采样长度t p为f s=1 T s,t p=N・T s.(3)设x a(nT s)的D T FT为X a(e jΞ),w R(n)的D T FT为W R(e jΞ),则不难得到:X a(e jΞ)=2Π∆(Ξ-Ξ1)e jΥ1,Ξ1=2Πf1 f s,W R(e jΞ)=W R(Ξ)e-j N-12Ξ,W R(Ξ)=sin(NΞ 2)N sin(Ξ 2).(4)根据Fou rier变换的性质,可以得到简谐信号采样序列x(n)的D T FT值X(e jΞ)为X(e jΞ)=12Π∫+Π-ΠX a(e jΗ)W R[e j(Ξ-Η)]dΗ= W R(Ξ-Ξ1) e j[Υ1-N-12Ξ].(5)当进一步应用D FT算法计算x(n)的频谱时,相当于在数字频域以2Π N的间隔(对应的模拟频域间隔为1 t p)对式(5)中的X(e jΞ)进行抽样.设模拟频率f1对应的数字频率为(k1+∆1)・2Π N,且k1为正整数,0≤∆1<1,则: (k1+∆1)=N・f1f s.(6)当∆1=0和∆1=0.5时,x(n)的D T FT幅频特性和D FT的幅频特性分别如图1所示.图1　同步和非同步采样后的离散频谱比较F ig.1　Single frequency signal’s spectrum samp led bysynch ronizati on and asynch ronis m图1中的连续曲线为D T FT的结果,符号“○”为D FT的抽样结果.由此可以得到以下结论:(1)简谐信号的理想频谱为单位冲激,而采样序列的D T FT谱与窗函数的谱一致,因而产生了泄漏;(2)当同步采样时,在各离散点上的D FT谱值与信号的理想谱值完全一致;(3)当非同步采样时,在各离散点上的D FT谱值与信号的理想谱值不一致,因而会产生泄漏.1.2　正弦(余弦)信号的插值FFT算法实际采样过程中,非同步采样很难避免,从而采样序列的频谱会产生泄漏.若原始信号为正弦(余弦)信号,则不论采样同步与否,可以对信号的D FT 算法结果作相应的处理,最终得到信号的真实频谱.设图1中信号的离散频谱在离散频点k1和k1+1处出现谱峰,由式(5)可得:X(e jΞ) Ξ=k12ΠN=W R(k1-(k1+∆1))2ΠN= sin(Π∆1) [sin(Π∆1 N) N].(7)311　第1期祁才君,等:应用插值FFT算法精确估计电网谐波参数X(e jΞ) Ξ=(k1+1)2ΠN=W R(k1+1-(k1+∆1))2ΠN=sin(Π∆1) [sin(Π(1-∆1) N) N].(8)在实际应用中采样点数N一般较大,同时考虑 ∆1<1,则sin(Π∆1) sin(Π∆1 N)=sin(Π∆1) [Π∆1 N+O(N-3)]≈sin(Π∆1) (Π∆1 N).(9)同理sin(Π∆1) sin(Π(1-∆1) N)=sin(Π∆1) [Π(1-∆1) N+O(N-3)]≈sin(Π∆1) [Π(1-∆1) N].(10)进一步,设Β1= X(e jΞ) Ξ=(k1+1)2ΠN X(e jΞ) Ξ=k12ΠN,则由式(9)和式(10)可推得∆1=Β11+Β1.(11)由此可得,信号的实际频率估计为f1=(k1+∆1) t p.(12)由式(9)或式(10)(建议用幅值大的式子)可推得信号的幅值A1和初相Υ1的估计分别为A1= X(e jΞ) Ξ=k12ΠN・(Π∆1) sin(Π∆1),(13)Υ1=ang le( X(e jΞ) Ξ=k12ΠN)-∆1Π(N-1) N.(14) 2　电网谐波参数的估计2.1　FFT算法应用中的泄漏问题实际电网信号往往并不是简谐信号,它具有如下特点:①含有丰富的谐波分量;②谐波分量的幅值一般仅为基波分量幅值的百分之几,或更小.当对电网信号进行非同步采样时,基波分量的频谱泄漏将严重影响2次、3次等谐波分量的频谱,从而导致谐波测量产生很大的误差.若相邻谐波之间的幅值相差过大,幅值大的谐波分量同样有可能淹没幅值小的谐波分量.若用矩形窗对信号进行截断,主瓣对旁瓣的泄漏最大达-13dB,对远隔2～10次旁瓣的泄漏也有约-30dB[1].假定2次谐波的幅值为基波的1.0%,即两者幅度相差-40dB.当采样长度为一个信号周期时,若直接应用D FT算法进行频谱分析,则基波对2次谐波的泄漏高达-13dB,完全淹没了-40 dB的实际谐波信号.为了减少D FT算法应用过程的频谱泄漏,可以改用其他窗函数对信号进行截断或增加采样长度.典型的窗函数为各种升余弦窗[1].其中,H ann ing窗的旁瓣最大泄漏为-31dB,对远隔旁瓣以-6dB 个的速度下降.H amm ing窗旁瓣最大泄漏为-42 dB,对远隔旁瓣以-2dB 个的速度下降.但是,这两种窗的主瓣半宽为矩形窗的2倍,对邻近一次谐波的泄漏十分严重,为了避免这个缺点,采样长度至少必须是信号周期的2倍以上.B lackm an窗旁瓣最大泄漏为-58dB,对远隔旁瓣以-6dB 个的速度下降.但是,该窗的主瓣半宽为矩形窗的3倍,对邻近一、二次谐波的泄漏十分严重,为了避免这个缺点,采样长度至少必须是信号周期的3倍以上.设窗函数的主瓣半宽为M・2Π 2,旁瓣最大泄漏为A1(d B),衰减速度为D(dB 个),电网谐波之间的幅度分瓣率要求为∃A(dB),则采样长度至少需要d个信号周期以上,且:d=(∃A-A1) D+M+0.5.(15)从各种窗函数的特性看,B lackm an窗虽然旁瓣衰减大,但其计算相对复杂.H amm ing窗的旁瓣衰减略比H ann ing窗大,但随旁瓣增加其衰减速度很慢.对电网信号的分析,若选用H ann ing窗,不但计算量较小,同时可以通过调节采样长度达到减少谐波间泄漏的目的.2.2　应用Hann i ng窗的插值FFT算法在满足采样定理的前提下,同时按式(15)选择H ann ing窗截断时的采样长度,则可以忽略电网各次谐波之间的频谱泄漏影响.利用D FT算法的线性性,可用插值FFT算法分别分析各次谐波的参数.H ann ing窗的窗函数表达式和D T FT结果为w H(n)=0.5-0.5co s2ΠnN-1R N(n),(16) n=0,1,…,N-1.W H(e jΞ)=W H(Ξ)e-jΞ(N-1) N,W H(Ξ)=0.5W R(Ξ)+0.25W RΞ-2ΠN-1+W RΞ+2ΠN-1.(17)为不失一般性,设第m次谐波信号为x m(t)=A m・e j(2Πf m t+Υm).并设谐波信号采样序列对应的离散频点为(k m +∆m)=N・f mf s,其中k m为整数;0≤∆m<1;f m是谐波频率;t p为采样长度.与矩形截断时的推导相类似,设加窗后采样序列频谱为X H(e jΞ),可得X H(e jΞ) Ξ=km2ΠN≈A msin(Π∆m)2∆m(1-∆2m)Π.(18)X H(e jΞ) Ξ=(km+1)2ΠN≈A msin(Π∆m)2∆m(1-∆m)(2-∆m)Π.(19)411浙　江　大　学　学　报(工学版) 第37卷　设Βm = X H (e j Ξ) Ξ=(k m +1)2ΠNXH(e j Ξ)Ξ=k m 2ΠN,则可得∆m =2Βm -11+Βm.(20)由式(17)和式(20)可得谐波幅值、频率、相位的估计式分别为 A m = X H (e j Ξ) Ξ=k m 2ΠN ・2Π∆m (1-∆2m )sin (Π∆m ),(21) f m =(k m +∆m ) t p ,(22)Υm =ang le X H (e j Ξ) Ξ=k m 2ΠN -∆m Π(N -1) N .(23)对电网波形而言,基波分量远大于各次谐波,谐波对基波的泄漏影响很小,而基波对谐波的泄漏影响很大,故基波频率的估计精度要远高于谐波频率的估计.式(20)～式(22)主要用作基波频率估计,谐波频率的估计可作如下修正:f m =m ・f 1,(24)∆m =m ・(k 1+∆1)-k m .(25)3　仿真结果及分析为了验证上述算法的准确性,假定电网波形可表达为x (t )=∑9m =1Amco s (2Πm f 1t +Υm ).(26)其中:基波频率设为49.5H z ,基波最大幅值为310V ,初相为5°.各次谐波的幅值按实际电网信号的特点设置,各谐波相位任意设定,具体参数如表1所示.表1　谐波参数的设置T ab .1　H ar monic param eters setup参数谐波 次23456789含有率 %0.53.00.45.00.32.00.21.0初相(°)10.020.0120.0150.0210.0240.0300.0330.0 由于最大的频谱泄漏是基波对2次谐波的影响,而表1中基波和2次谐波的幅值相差-46dB ,故若希望将泄漏限止在0.1%(即-60dB )以内,则∃A 至少为106dB .由式(15)可得H ann ing 窗截断时的窗函数长度至少为d =15个信号周期.若采样点数N 设定为1024,采样频率取f s =2500H z ,则以电网额定频率(50H z )为基准,实际采样信号的周期个数达约20个,符合限制泄漏要求.基于普通FFT 算法和基于H ann ing 窗的插值FFT 算法的电网参数估计结果如表2和表3所示.表2　频率参数的估计结果T ab .2　T he result of the pow er frequency esti m ati on 谐波次设定值普通FFT 算法估计值估计偏差 %插值FFT 算法估计值估计偏差 %149.551.273.57649.5010.002299.0102.543.57699.0020.0023148.5151.371.933148.5030.0024198.0202.642.343198.0040.0025247.5249.020.614247.5050.0026297.0297.850.286297.0060.0027346.5349.120.756346.5070.0028396.0397.950.492396.0080.0029445.5449.220.835445.5090.002表3　幅值参数的估计结果T ab .3　T he result of the amp litude esti m ati on 谐波 次设定值普通FFT 算法估计值估计偏差 %插值FFT 算法估计值估计偏差 %1310.0272.4-12.15309.8-0.071821.555.034225.31.5580.534239.309.3210.2269.3140.150341.241.80545.541.239-0.0576515.512.75-17.7215.510.031260.931.32041.970.9320.166576.206.4724.3826.2280.452580.621.245100.90.620-0.003193.102.252-27.353.1060.1904表4　相位参数的估计结果T ab .4　T he result of the phase esti m ati on 谐波 次设定值普通FFT 算法估计值估计偏差 %插值FFT 算法估计值估计偏差 %15.0174.233844.919-1.6143210.0299.228929.476-5.2381320.088.52342.619.74-1.08194120.0100.5-16.27119.6-0.29885150.084.79-43.47149.6-0.27016210.0199.4-5.070209.5-0.25367240.077.87-67.56239.4-0.23588300.0283.5-5.517299.3-0.21769330.04.658-98.59329.3-0.2206分析表2至表4的结果,可得出以下结论:(1)基于插值FFT 算法的谐波参数分析精度明显优于基于普通FFT 算法的谐波参数分析精度;(2)基于普通FFT 算法的谐波频率(尤其是基频)、谐波幅值分析精度达不到电能质量的国家标准,而基于插值FFT 算法的频率分析精度可以控制在0.01%以内,幅值分析精度可控制在0.5%以内,完全可以满足国家标准;(3)普通FFT 算法的各次谐波相位参数估计几乎没有参考价值,而插值FFT 算法的相位估计精度可达5%;(4)与文[8]相比,按公式(25)修正后,参数估计精度有了进一步地提高.511　第1期祁才君,等:应用插值FFT 算法精确估计电网谐波参数4　结　语本文提出的基于H ann ing窗的插值FFT算法具有较高的谐波参数估计精度.随着采样长度的增加,估计精度还可进一步提高,并且对幅值、相位、频率的估计都有比较简单的计算公式.显然该方法同样适合于其他周期信号或多谐信号的参数估计.本算法的不足之处是分析窗的宽度一般要达十几个信号周期,参数估计的实时性不够理想.另外,当信号中包含噪声时,如何提高参数估计准确度和精度还值得作进一步的研究.参考文献(References):[1]M ITRA S K.D ig it al signal processi ng-a co m puter-based approach[M].Beijing:T singhua U niversity P ress,2001.446-460.[2]周箭,陈隆道.新型电力参数测试与分析[J].中国电力,1998,31(11):18-20,36.ZHOU J ian,CH EN L ong2dao.A new pow er param eter testing and analyzing system[J].Electr ic Power,1998, 31(11):18-20,36.[3]陈隆道,钱照明,张圣训.周期域分析中非同步取样数据的同步化[J].电子学报,2001,29(7):950-953.CH EN L ong2dao,Q I AN Zhao2m ing,ZHAN G Sheng2 xun.Synch ronizati on algo rithm fo r asynch ronous sam2p led data set in peri od dom ain signal analysis[J].ACTA Electron ica SI N I CA,2001,29(7):950-953.[4]戴先中.进一步提高准同步采样谐波分析准确度的两种方法[J].仪器仪表学报,1992,13(4):350-357.DA I X ian2zhong.M ethods fo r i m p roving the accuracies of quasi2synch ronous samp ling har monics analysis[J].Ch i nese Journal of Sc ien tif ic I n stru m en t,1992,13(4): 350-357.[5]张建秋,沈毅.非整周期采样算法的误差分析[J].哈尔滨工业大学学报,1995,27(5):99-103.ZHAN G J ian2qiu,SH EN G Y i.T he erro r analysis of the non2integer2peri od samp ling algo rithm[J].Journal of Harbi n I n stitute of Technology,1995,27(5):99-103.[6]JA I N V K,COLL I N S W L,DAV IS D C.H igh2accura2cy analog m easurem ents via interpo lated FFT[J].I EEE Tran s I M,1979,28(2):113-122.[7]THOM A S G.Interpo lati on algo rithm s fo r discretefourier transfo r m of w eighed signals[J].I EEE Tran sI M,1983,32(2):350-355.[8]ZHAN G Fu2sheng,GEN G Zheng2xing,YUAN W ei.T he algo rithm of interpo lating w indow ed FFT fo r har2 monic analysis of electric pow er system[J].I EEE Tran s on Power D elivery,2001,16(2):160-164.[9]BEGOV I C M M,DJU R I C P M,DUNLA P E S,et al.F requency track ing in pow er netw o rk s in the p resenceof har monics[J].I EEE Tran s on Power D elivery, 1993,8(2):480-486.下期论文摘要预登共享总线多微处理器系统中分布式仲裁方案的研究及实现刘　峰,葛霁光(浙江大学生物医学信息研究所,浙江杭州310027)摘　要:共享总线多微处理器系统具有相当高的性能价格比,但随着处理器速度的不断提高,总线通信带宽和可靠性成为限制并行处理能力的主要瓶颈.总线仲裁对提高系统性能具有关键的作用.分布式一致优先权机制是一种新型的仲裁机制模型.在此基础上提出了基于时间变化的优先权仲裁事务缓冲的分布式仲裁方案.仲裁缓冲机制能够充分利用系统总线的并发特性来提高总线的利用率,此机制有利于并行处理系统中的实时和总线重负载应用环境.虚拟令牌环和时变优先权策略为处理器节点公平地共享总线资源提供了保障.多字节优先权标识使优先权更易于扩展.仿真工作验证了所提方案的正确性.关键词:并行处理;仲裁缓冲;虚拟令牌环;实时仲裁;优先权队列611浙　江　大　学　学　报(工学版) 第37卷　。

摊
０ 引言
电力系 统 中的频率 是 变化 的 ．但是 数字 化变 电
于这种窗函数的高精度谐波分析方法 ，利用曲线拟 合 方法 推 导 出简洁 的修 正公式 。对含 ２ １ 次谐波 的复 杂信号进行仿谐波分析 ， 试验结果表明： 在非同步采 样和非整周期截断的条件下 ，本文算法能显著提高 基波频率 、 各次谐波幅值和相位的计算精度 ， 提高信 号 分析 准确 性 ， 适 合 于高 精度 谐波 分析 ｏ １ 。
ｉ ｎ ｔ ｅ ｒ ｐ ｏ ｌ ａ ｔ ｉ ｏ ｎ ｆ ａ ｓ ｔ Ｆ ｏ ｕ ｒ ｉ ｅ ｒ ｔ ｒ ａ ｎ ｓ ｆ ｏ ｘ ｍ （ Ｆ Ｆ Ｔ ）ｉ ｓ ｐ ｒ ｏ ｐ ｏ ｓ ｅ ｄ ， ｔ ｈ ｅ ａ ｐ ｐ ｌ ｉ ｃ ａ ｂ ｌ ｅ ｒ ｅ ｃ ｔ ｉ ｉ ｆ ｃ ａ ｔ ｉ ｏ ｎ ｆ ｏ ｒ ｍ ｕ ｌ ａ ｓ ｏ ｆ ｔ ｈ ｅ ｄ ｏ ｕ ｂ ｌ ｅ － ｓ ｐ ｅ ｃ ｔ ｒ ｕ ｍ － ｌ ｉ ｎ ｅ ｉ ｎ ｔ ｅ ｒ ｐ ｏ ｌ ａ ｔ ｉ ｏ ｎ
０ １
一
种 高精 度加 窗插值Ｆ Ｆ Ｔ 谐 波 分析方 法
汪 立， 刘 翠琳
（ 武汉 大 学 电气工程 学院 ， 湖北 武汉 ４ ３ ０ ０ ７ ２ ）
摘
要： 非 同步采样和 非整数周期截断造成的频谱泄漏会影响谐波测量结果 的准确性。提 出了一种高精度加 窗
Ｆ Ｆ ｒ 插值谐波分析方法。介 绍了一种余弦组合窗函数， 讨论 了该余弦组合 窗的特性 ， 并首次将该 窗函数运用在谐
ｉ ｓ ｏ ｂ ｔ ａ ｉ ｎ ｅ ｄ ｂ ｙ ｕ ｓ ｉ ｎ ｇ ｐ ｏ ｌ ｙ ｎ ｏ ｍｉ ａ ｌ ｃ ｕ ｒ ｖ ｅ ｉ ｆ ｔ ｆ ｕ ｎ ｃ ｔ ｉ ｏ ｎ ｓ ． Ｔ ｈ ｅ ｔ ｅ ｓ ｔ ｓ ｈ ｏ ｗｓ ｔ ｈ ａ ｔ ｔ ｈ ｅ ｎ ｅ ｗ ｍｅ ｔ ｈ ｏ ｄ ｈ ａ ｓ ｈ ｉ ｇ ｈ ｅ ｒ ｃ ａ ｌ ｃ ｕ ｌ ａ ｔ ｉ ｏ ｎ ａ ｃ ｃ ｕ ｒ ａ ｃ ｙ ａ ｎ ｄ ｉ ｓ ｍｏ ｒ ｅ ｐ ｒ ａ ｃ ｔ ｉ ｃ ａ ｌ ｉ ｎ ｔ ｈ ｅ ｅ ａ ｌ ｃ ｕ ｌ ａ ｔ ｉ ｏ ｎ ｏ ｆ ｆ ｒ ｅ ｑ ｕ ｅ ｎ ｃ ｙ ， ａ ｍｐ ｌ ｉ ｔ ｕ ｄ ｅ ａ ｎ ｄ ｐ ｈ ａ ｓ ｅ ．

f (m1,1)
f (m1,n1)
矩阵表示的几点说明
▪ g (i, j)代表(i, j)点的灰度值，即亮度值。
▪ 以上数字化有以下几点说明：
▪ （1） 由于g (i, j)代表该点图像的光强度，而光
是能量的一种形式，故g (i, j)必须大于零，且为 有限值，即： 0＜g (i, j)＜∞。
▪ （2） 数字化采样一般是按正方形点阵取样的，
➢ 量化的准则不同，会导致不同的量化效果。 ➢ 从不同的角度将量化方法分成4类：
➢ （1）按量化级步长均匀性 ➢均匀量化和非均匀量化。
➢ （2）按量化对称性 ➢ 对称量化和非对称量化
➢ （3）按量化时采样点相互间的相关性分
➢无记忆和有记忆量化。
➢ （4）按量化时处理的采样点数分
➢标量量化和矢量量化。
采样定理
采样定理 一个频带限制在(0, fH)Hz内的时间连续信号m(t), 如果以大于 或等于2fH的采样速率对它进行等间隔采样，则可以由采样 序列无失真地恢复原始信号m(t)。 如采样间隔Ts=1/(2fH)，则Ts被称为奈奎斯特（Nyquist）间隔 。如采样频率fs=1/Ts，则fs被称为奈奎斯特（Nyquist）速率 。 采样定理告诉我们，若采样速率f < fs（或采样间隔T>Ts）， 则会产生失真，这种失真称为混叠失真。
一维向量量化
▪ 目的：对任一输入矢量X，在码书中寻找最佳匹
配码矢Xi。
▪ 常用的最佳匹配原则：寻求最小误差。
▪ 若码书尺寸为M，矢量X对应码矢Xi，信号矢量X
的概率密度函数为p(X)，则总的量化误差可表示
为
ae(X,Xi)p(X)
常用的误差有
i
：
均方绝对值误差（MAE）

e−
W jϕh R
(2π(
f
+ rh f1) /
fs )
(2)
式中WR(·)为矩形窗的频谱函数，表达式为
WR
(
f
)
=
sin(Nf sin( f
/ 2) / 2)
e− j( N −1)
f
/2
(3)
对式(2)离散抽样，并忽略负频点−f1处频峰的旁 瓣影响，得到信号加矩形窗后的离散傅里叶变换表
达式为
∑ X (kΔf ) =
谱，图中Mmag为幅值，ω为角频率。采用加窗插值 FFT进行谐波分析时，旁瓣电平小且渐近衰减速率
大的窗函数能抑制频谱泄漏的影响，提高谐波分析
的精度。由图 3 可见，5 项Rife-Vincent(I)窗具有最
理想的旁瓣特性，旁瓣峰值电平为−74.5 dB，旁瓣
渐进衰减速率达到 30 dB/oct。因此，本文采用 5 项
数周期截断后的频谱如图 2 所示。可看出，信号峰 值观测点出现偏差Δx，即短范围泄漏。同时，信号 时域截断时相邻频谱旁瓣间的相互干扰造成长范
围泄漏。
1.0
0.8
Ah/V
0.6
0.4
0.2
ki
ki'
0.0 01234567 k
图 1 同步采样和整周期截断后的频谱 Fig. 1 Spectrum under coherent sampling and
本文对信号在同步采样和整周期截断及非同
116
中国电机工程学报
第 29 卷
步采样和非整周期截断时的频谱进行分析，讨论 Rife-Vincent窗[20-21]的频谱特性，提出基于 5 项 Rife-Vincent(I)窗双谱线插值FFT的谐波分析方法， 利用曲线拟合函数推导修正公式。对弱信号分量的 提取和含 21 次谐波的复杂信号进行仿真的结果表 明：在非同步采样和非整周期截断的条件下，本文 算法能有效消除各次谐波相互干扰，提高信号分析 精度，适合于高准确度谐波分析。

— 45 —
低压电器（2010№9）
·电能质量·
同的影响，这就需要根据具体的情况选择合适的 窗函数，从而达到理想的测试效果。常用的窗函 数有矩形窗、三角形窗、汉宁窗、海明窗和布莱克 曼窗。
选择窗函数时通常要求满足：① 窗谱主瓣尽 可能窄，以获得较陡的过渡带；② 尽量减少窗谱 的最大旁瓣的相对幅度，也就是能量尽量集中于 主瓣，这样使肩峰和纹波减少，就可增大阻带的衰 减。但这两项要求是无法同时满足的，往往是增 加主瓣宽度以换取对旁瓣的抑制。
— 46 —
幅值，另外，从式中可看出当频率分辨率越高时， 得到的距离峰值频率附近的 2 条频谱（ 即幅值最 大和次大的谱线） 之间的频率间隔越小，经插值 后得到的结果越准确。根据 GB / T 15945—2008 《电能质量 电力系统频率偏差》的规定，电力系统 正常频率允许偏差值为 ± 0. 2 Hz，当系统容量较 小时，偏差可放宽到 ± 0. 5 Hz。因此，在考虑算法 实现时假设信号的基波频率为50 Hz，经 FFT 计算 以后，该基波频率对应的最大谱线幅值位于 40 ～ 60 Hz，可在该频率范围内通过搜寻最大值的方法 找到这条幅值最大的谱线，然后再通过对该幅值 最大的谱线相邻的谱线比较的方法找出次幅值最 大的曲线，最后通过这两条谱线计算出插值点的 位置，进而应用相应的插值公式准确地求出信号 中基波分量的幅值、频率和相位。在求得基波的 频率之后可通过倍频的方法来求得各次谐波的频 率（ 假设测出的基波频率为 f0 ，则第 n 次谐波的频 率为 nf0 ，根据频率分辨率的大小 Δf（1 / T） 可得到 第 n 次谐波在频谱上相应的最大谱线幅值和次大 谱线幅值的大概位置，通过幅值比较找出这两条 谱线，进而通过与窗函数对应的双峰谱线插值修 正算法，求出第 n 次谐波的精确幅值和相位） 。

一种加窗插值FFT谐波分析方法摘要：由于很难实现同步采样和整周期截断，因此，利用fft算法分析电网谐波信号时存在频谱泄露和栅栏效应，影响算法的分析精度。

加窗插值fft是抑制频谱泄露和消除栅栏效应的有效方法，在此提出一种基于3项3阶nuttall窗插值fft的谐波分析方法，推导了插值系数公式以及各次谐波的频率、幅值和相位的修正公式。

对该算法与hanning窗、blackman窗插值fft算法进行matlab 仿真对比研究，验证了该算法具有更高的分析精度。

关键词：谐波； fft；窗函数；插值；电力系统引言随着大量电力电子装置和非线性负载在电力系统中的广泛应用，使电网中产生了大量的高次谐波［12］，严重威胁电网的电能质量和用户设备的安全运行，因此谐波的准确测量具有重要意义。

快速傅里叶变换(fft)是最主要的电力谐波分析方法，但由于实际工程应用中很难实现同步采样和整周期截断，因此，fft方法存在频谱泄露和栅栏效应，影响谐波分析精度。

加窗插值fft算法是抑制频谱泄漏和消除栅栏效应的有效方法，其原理是通过加窗运算抑制频谱泄露、通过插值运算消除栅栏效应。

常用的窗函数有hanning 窗［35］、blackman窗［6］、blackman harris窗［78］等。

hanning窗的插值公式简单，计算量小，但分析精度较低；blackman和blackman harris窗插值fft算法的分析精度高，但插值公式复杂，计算量大。

本文首先分析了fft算法频谱泄露的原因和3项3阶nuttall窗函数的特点，然后推导出其插值fft算法的计算公式，其插值系数具有简单的显式表达式，谐波的频率、幅值和相位的修正公式简单明了，易于实现。

同hanning窗和blackman窗插值fft算法的仿真对比研究结果表明，所提出算法更加适合于电力系统谐波的精确测量。

1fft频谱泄漏的原因离散傅里叶变换(dft)变换是针对有限长序列信号进行傅里叶变换的一种数值分析方法。

W ( ) Wr ( )
W () N / 2
（7） （8）
可见，Hanning 窗频谱（在主瓣下）具有严格的线性相位。而同样是余弦窗的 Hamming 窗 却不具有线性相位，这使得其频谱表达式难以简化，从而得不到实用的插值公式。余弦窗的 一般表达式为[1]
K 2 wK (n) (1)k ak cos kn , n 0, 1, , N 1 N k 0
其中，
（13）
Wr ( ) 0.5
sin( ) sin( / N )
cos( / N ) sin 2 ( / N ) 1 sin 2 ( / N )
（14）
此处，我们关心的频率范围是 1 1 。当 N 比较大时，利用
sin( / N ) / N
Fig. 1
图 1 N 16 时，近似表达式与精确表达式的比较 Comparison between approximate and exact expression (when N = 16)
图 2 是近似表达式的误差曲线。 从图 2 （b） 可以得到一个重要结论是： 即使在 N 不太大时， 插值算法也能达到相当高的精度。比如，当信号的频率接近 Nyquist 频率时，这时如果按照 消除频谱泄漏影响的原则选取截断长度内信号的周期数，则截断长度 N 不会很大，此时， 同样可以用插值算法来克服栅栏效应。 仿真部分将用一个检测电网载波信号的例子来验证这 种情况下的应用。
（9）
系数满足约束条件[5]
a
k 0
K
k
1
（10）
当 a0 1 时 ， w0 (n) 为 矩 形 窗 ； 当 a0 a1 0.5 时 ， w1 (n) 为 Hanning 窗 ； 当 当 a0 0.42, a1 0.50, a2 0.08 时， a0 0.54, a1 0.46 时， w1 (n) 为 Hamming 窗； w2 (n) 为 Blackman 窗。容易证明，除矩形窗外，相位线性的充分条件是

第２９卷第５期２０１２年１０月现　代　电　力Ｍｏｄｅｒｎ　Ｅｌｅｃｔｒｉｃ　ＰｏｗｅｒＶｏｌ．２９　Ｎｏ．５Ｏｃｔ．２０１２文章编号：１００７－２３２２（２０１２）０５－００２８－０４文献标识码：Ａ中图分类号：ＴＭ７１１基于加窗插值ＦＦＴ算法的间谐波检测方法研究王晓亮，李　娜（山东大学电气工程学院，山东济南　２５００６１）Ｒｅｓｅａｒｃｈ　ｏｎ　Ｉｎｔｅｒ－ｈａｒｍｏｎｉｃ　Ｄｅｔｅｃｔｉｏｎ　Ｕｓｉｎｇ　Ｗｉｎｄｏｗｓ　ａｎｄＩｎｔｅｒｐｏｌａｔｅｄ　ＦＦＴ　ＡｌｇｏｒｉｔｈｍＷＡＮＧ　Ｘｉａｏｌｉａｎｇ，ＬＩ　Ｎａ（Ｓｃｈｏｏｌ　ｏｆ　Ｅｌｅｃｔｒｉｃａｌ　Ｅｎｇｉｎｅｅｒｉｎｇ，Ｓｈａｎｄｏｎｇ　Ｕｎｉｖｅｒｓｉｔｙ，Ｊｉｎａｎ　２５００６１，Ｃｈｉｎａ）摘　要：电网中间谐波的存在，会对电能质量以及供电可靠性带来不利影响，故准确检测间谐波对电力系统稳定运行意义重大。

根据间谐波特性，在一般ＦＦＴ算法基础上，提出了基于加窗插值ＦＦＴ算法的间谐波检测方法。

通过分析对比不同窗函数的特点，选取检测精度较高的Ｈａｎｎｉｎｇ窗作为所加分析窗，同时确定所加窗函数的宽度及采样周期，可准确检测出系统中的谐波及间谐波。

在ＭＡＴＬＡＢ环境下仿真得到一般ＦＦＴ算法及加窗插值ＦＦＴ算法对谐波和间谐波的检测结果，通过对所得频率和幅值估计结果的对比分析可知，加窗插值ＦＦＴ算法检测精度更高、实用性更强。

关键词：间谐波；傅里叶变换；快速傅里叶变换；加窗插值Ａｂｓｔｒａｃｔ：Ｔｈｅ　ｉｎｔｅｒ－ｈａｒｍｏｎｉｃｓ　ｉｎ　ｐｏｗｅｒ　ｇｒｉｄ　ｍａｙ　ｃａｕｓｅ　ｎｅｇ－ａｔｉｖｅ　ｉｎｆｌｕｅｎｃｅ　ｏｎ　ｐｏｗｅｒ　ｑｕａｌｉｔｙ　ａｎｄ　ｐｏｗｅｒ　ｓｕｐｐｌｙ　ｒｅｌｉａｂｉｌｉ－ｔｙ，ｓｏ　ｔｈｅ　ａｃｃｕｒａｔｅ　ｄｅｔｅｃｔｉｎｇ　ｏｆ　ｉｎｔｅｒ－ｈａｒｍｏｎｉｃｓ　ｉｓ　ｉｍｐｏｒｔａｎｔｆｏｒ　ｔｈｅ　ｓｔａｂｌｅ　ｏｐｅｒａｔｉｏｎ　ｏｆ　ｔｈｅ　ｐｏｗｅｒ　ｓｙｓｔｅｍ．Ａｃｃｏｒｄｉｎｇ　ｔｏｔｈｅ　ｉｎｔｅｒ－ｈａｒｍｏｎｉｃｓ　ｃｈａｒａｃｔｅｒｉｓｔｉｃｓ，ｂａｓｅｄ　ｏｎ　ｔｈｅ　ＦＦＴ　ａｌｇｏ－ｒｉｔｈｍ，ｉｎｔｅｒｐｏｌａｔｉｏｎ　ＦＦＴ　ａｌｇｏｒｉｔｈｍ　ｗａｓ　ｐｒｏｐｏｓｅｄ　ｔｏ　ｄｅｔｅｃｔｉｎｔｅｒ－ｈａｒｍｏｎｉｃｓ．Ｔｈｒｏｕｇｈ　ｔｈｅ　ａｎａｌｙｓｉｓ　ｏｆ　ｔｈｅ　ｃｈａｒａｃｔｅｒｉｓｔｉｃｓｏｆ　ｄｉｆｆｅｒｅｎｔ　ｗｉｎｄｏｗ　ｆｕｎｃｔｉｏｎｓ，Ｈａｎｎｉｎｇ　ｗｉｎｄｏｗ　ｉｓ　ｓｅｌｅｃｔｅｄａｓ　ｔｈｅ　ａｎａｌｙｓｉｓ　ｗｉｎｄｏｗ，ａｎｄ　ｔｈｅ　ｗｉｎｄｏｗ－ｌｅｎｇｔｈ　ａｎｄ　ｔｈｅ　ｓａｍ－ｐｌｉｎｇ　ｐｅｒｉｏｄ　ａｒｅ　ｄｅｔｅｒｍｉｎｅｄ，ｗｈｉｃｈ　ｃａｎ　ｄｅｔｅｃｔ　ｈａｒｍｏｎｉｃｓａｎｄ　ｉｎｔｅｒ－ｈａｒｍｏｎｉｃｓ　ａｃｃｕｒａｔｅｌｙ．Ｂｅｓｉｄｅｓ，ＦＦＴ　ａｌｇｏｒｉｔｈｍ　ａｎｄｉｎｔｅｒｐｏｌａｔｉｏｎ　ＦＦＴ　ａｌｇｏｒｉｔｈｍ　ａｒｅ　ｓｉｍｕｌａｔｅｄ　ｂｙ　ＭＡＴＬＡＢ．Ｂｙｃｏｎｔｒａｓｔｉｎｇ　ａｎｄ　ａｎａｌｙｚｉｎｇ　ｔｈｅ　ｅｓｔｉｍａｔｅｄ　ｒｅｓｕｌｔｓ　ｏｆ　ｆｒｅｑｕｅｎｃｙａｎｄ　ａｍｐｌｉｔｕｄｅ，ｉｔ　ｃａｎ　ｂｅ　ｓｅｅｎ　ｔｈａｔ　ｔｈｅ　ｗｉｎｄｏｗｓ　ａｎｄ　ｉｎｔｅｒｐｏ－ｌａｔｅｄ　ＦＦＴ　ａｌｇｏｒｉｔｈｍ　ｃａｎ　ｄｅｔｅｃｔ　ｈａｒｍｏｎｉｃｓ　ａｎｄ　ｉｎｔｅｒ－ｈａｒｍｏｎ－ｉｃｓ　ｏｆ　ａ　ｇｉｖｅｎ　ｓｉｇｎａｌ　ｍｕｃｈ　ｍｏｒｅ　ａｃｃｕｒａｔｅｌｙ．Ｋｅｙ　ｗｏｒｄｓ：ｉｎｔｅｒ－ｈａｒｍｏｎｉｃｓ；Ｆｏｕｒｉｅｒ　ｔｒａｎｓｆｏｒｍ；ＦＦＴ；ｗｉｎ－ｄｏｗｓ　ａｎｄ　ｉｎｔｅｒｐｏｌａｔｅｄ　ＦＦＴ０　引　言随着大量非线性电气设备的引入，电力系统谐波和间谐波污染日渐严重。

fft算法离散频谱校正FFT（Fast Fourier Transform）算法是一种快速计算离散傅里叶变换（Discrete Fourier Transform）的算法。

它的主要思想是通过对称性将N点的DFT分解为两个N/2点的DFT，再通过重组得到最终结果。

该算法的时间复杂度为O(NlogN)，相较于传统的DFT算法，其计算速度更快，因此广泛应用在信号处理、图像处理、通信等领域中。

离散频谱校正是指在频域中对信号进行处理，以消除或校正频谱中的不良效应。

在进行频域处理时，可能会出现混叠效应（频谱重叠）或频率偏移等问题，这会导致信号的失真或干扰。

离散频谱校正的目的就是通过一系列算法和技术，对频谱进行调整和修正，以恢复信号的原始特性。

离散频谱校正的方法有很多种，下面将简要介绍几种常见的方法。

1. 频谱外插频谱外插是一种常见的频谱校正方法，它通过在频谱中插入一定数量的零值来改变信号的频谱特性。

这样可以使频谱变得更加平滑，并且减小混叠效应。

频谱外插在FFT算法中很容易实现，只需要将原始信号补零到2的幂次方长度即可。

2. 频谱滤波频谱滤波是指通过滤波器对频谱进行处理，以去除或衰减不需要的频率分量。

常见的滤波器包括低通滤波器、高通滤波器、带通滤波器等。

滤波器可以选择不同的截止频率、通带宽度和阻带宽度，以满足不同的要求。

3. 频谱修正频谱修正是一种校正频谱幅度和相位的方法。

通常在进行频域分析时，频率响应对于不同频率的信号可能有不同的增益和相位差，这就需要进行补偿和修正。

频谱修正的方法包括经验修正和数学模型修正等，可以根据信号的特性进行选择。

4. 非线性变换非线性变换是一种通过对频谱进行非线性操作，以改变频谱特性的方法。

常见的非线性变换包括幂律变换、对数变换、绝对值变换等。

非线性变换可以改变频谱的动态范围和分辨率，从而提取出信号的细节或增强信号的特征。

5. 频率域插值频率域插值是指通过对频谱进行插值，以增加频率的分辨率或减小频率的间隔。

一种FFT插值正弦波快速频率估计算法对被噪声污染的正弦波信号进行频率估计是信号参数估计中的经典问题，目前国内外已提出不少方法。

文献给出了在高斯白噪声中对正弦波信号频率进行最大似然估计算法，该算法能够达到卡拉美-罗限(CRB)，但计算量大，实现困难。

FFT频率估计方法具有速度快、便于实时处理的特性而得到了广泛应用。

但FFT频率估计方法得到的是离散频率值，当信号频率与FFT离散频率不重合时，由于FFT的栅栏效应，信号的实际频率应位于两条谱线之间。

显然仅仅利用FFT幅度最大值估计信号频率难以满足精度要求，因此各种插值算法应运而生。

文献给出了Rife算法，在对输入信号进行一次FFT运算后，利用最大谱线及其相邻的一根次大谱线进行插值来确定真实频率位置。

当信号的真实频率处于两相邻量化频率之间的中心区域时，Rife算法精度很高，但是在FFT量化频率附近的误差却较大。

文献提出了一种修正Rife算法，通过对信号进行频移，使新信号的频率位于两个相邻量化频率点的中心区域，然后再利用Rife算法进行频率估计。

文献提出了基于傅里叶系数插值迭代的频率估计方法，该方法能够有效提高精度，但需要多次串行迭代，不利于发挥FPGA并行处理的优势。

本文分析了以上3种算法的特点，并以之为基础结合FPGA的并行处理优势，提出了一种利用信号FFT插值系数的幅度和相位信息来构造频率修正项的新算法。

1 基于FFT插值的正弦波频率估计法1．1 算法原理单一频率正弦信号表示为：式中：A，f0，分别为正弦信号的幅度、频率和初相；fs为采样频率。

目前基于FFT的正弦信号频率估计分为2个过程来实现：粗测频和精测频。

粗测频通过直接观察FFT幅谱最大值点m来完成，受观测时长T的限制，误差范围为l／(2T)。

假设为信号频率的真实值，为信号频率与其FFT幅度最大处对应频率的相对偏差，m，与的关系如式(2)所示：考虑到FPGA并行计算的特点，利用流水线结构同时计算多个Xm+p，Xm+p-1值，将串行迭代变为并行迭代，其运算步骤归纳如下：。

基于插值FFT 算法的间谐波参数估计Interharmonics Estimation Based on Interpolation FFT Algorithm祁才君　王小海(浙江大学电气工程学院　310027)Qi Caijun Wang Xiaohai (Zhejiang University 310027　China ) 摘要　间谐波是频率介于两个谐波之间的信号。

间谐波除具有一般谐波信号的特性外,还会严重影响现有谐波补偿装置,使谐波补偿失败,因此准确检测间谐波的参数具有十分重要的意义。

讨论了间谐波的特点及检测方法,提出了基于加窗插值FFT 算法的间谐波参数估计;讨论了窗函数的选择规则,推导了分析窗宽度的估计公式,以及基于Hanning 窗的间谐波频率、幅值和相位的显式估计公式。

仿真结果证明,该算法对电网间谐波和谐波的幅度、频率和相位的估计在一定条件下具有很高的精度。

关键词:间谐波　谐波分析　电能质量　FFT 算法中图分类号:TM711A bstract Interharmonics are defined as the spectrum components that have frequencies between two har -monics .Comparing to har monics ,interharmonics have serious impacts on usual harmonics compensation de -vices .It is very important to measure interharmonics exactly .This paper discusses interhar monics 'character -istics ,measurement methods and proposes interharmonics estimation based on interpolation FFT algorithm .Some explicit estimation formula for the tapered -window -length and interharmonics 'parameters are given .The simulating result sho ws that the estimation of the interhar monics ',or har monics ',frequenc y ,a mplitude and phase reached a reasonable accurac y .Keywords :Interharmonics ,harmonics analysis ,po wer quality ,FFT algorithm祁才君　男,1964年生,副教授,主要从事数字信号处理,计算机实时控制,ASIC 电路设计等方面的研究。

ｉｔ ｒ ｏ ａｉ ｎ ｍｅ ｈ ｄ ｉ ｇｖ ｎ ａ ｄ ｃ ｍｐ ｒ ｄ ｔ ｈ ｏ ｎ ｅ ｐ ｌ ｔ ｔｏ ｓ ｉｅ ｎ ｏ ａ ｅ ｈ ｔ ｅ ｃ ｍｍｏ ｎ ｅ ｅｎ ｈ ｉ ｌｔ ｎ ｒ ｓ ｌ ｈ ｗ ｔ ａ ｈ ｍｐ ｏ ｅ ｏ ｎ ｏ ｅ ｈ ｒ ｉ ．Ｔ ｅ ｓ ｍｕ ａｉ ｅ ｕ ｔ ｓ ｏ ｈ ｔ ｔ ｅ ｉ ｒｖ ｄ ｏ ｓ ｌ ｏ ｔ ｍ ａ ｒ ｃ ｓ ｌ ｅ ｒ ｕ ｎ ｙ ｃ ａ ａ ｔ ｒ ｔ ｆｔ ｉ ｎ ｗｉ ｕ ｄ ｉ ｇ ｗｏ ｋ ｏ ｄ ｈ ｐ ｌ ｔ ｉ ｒｖ ｄ ａ ｇ ｒ ｈ ｃ ｎ ｐ ｅ ｉｅｙ ｇ ｔｔｅ ｆ ｑ ｅ ｃ ｈ ｒ ｃ ｅ ｓｉｓ ｏ ｅ ｓｇ ａ ｔ ｏ ｔａ ｄ ｎ ｒ ｌ ａ ．Ｔ ｅ ａ ｐ ｉａｉ ｎ ｏ ｍｐ ｏ ｅ ｉ ｈ ｅ ｉ ｃ ｈ ｌ ｈ ｃ ｏ ｆ ａｇ ｒ ｈ ｏ ｅ ａ ｔ ａ ｎ ｌ ｓ ｆｈ ｄ ｏ ｌ ｃｒ ｅ ｉ ｒ ｔ ｎ ｓｇ ａ ｌｏ ｓ ｏ ｓ ｔ ａ ｅ ｒ ｓ ｌ ｅ ｃ ｎ ｉｔ ｎ ｉ ｈ ｅ ｌ ｌ ｏ ｔ ｍ ｎ ｔ ｃｕ ｌａ ａｙ ｉ ｏ ｙ ｒ ｅｅ ｔｉ ｓ ｔｖ ｂ ａｉ ｉ ｎ ｌａ ｓ ｈ ｗ ｔｔ ｅ ｕ ｔ ａ ｏ ｓｓｅ ｔｗ ｔ ｔ ｅ ｒ ａ ｉ ｈ ｓ ｃ ｏ ｈ ｈ ｓｒ ｈ
Ａｂ ｔ ａ ｔ ｓｒ ｃ ：Ｒｅ ａ ｄ ｎ ｈ ｔ ｔ ｅ ａ ｔ ａ ｉｎ ｌ ｗｈｌ ｔ ｆｅ ｕ ｎ ｙ ｓ ｅ ｔｕ ｉ ａ ａｙ ｅ ｓ ａｌ ｏ ｓ ＇ ｍｅ ｔ ｔ ｅ ｄ ｍａ ｄ ｏ ｇ ｒ ｉｇ ｔ ａ ｈ ｃｕ ｌｓｇ ａ ｉ ｉ ｒ ｑ ｅ ｃ ｐ ｃ ｒ ｍ ｓ ｎ ｚ ｄ ｕ ｕ ｌ ｄ ｅ ｎｔ ｅ ｓ ｌ ｙ ｅ ｈ ｅ ｎ ｆ
（ｃ ｏｌｆ ｎｒｙａｄＰ ｗ ｒ ｎ ｉｅｒ ｇＨ ａｈ ｎ ｎ ｅｓｙ ｆ ｃ ｎ ｅ ｎ ｅｈ ｏ ｇ， ｈ ｎ ３ ０ ４ Ｈ ｂ ｉＣ ｉａ Ｓｈ ｏｏ ｅ ｎ ｏ ｅ Ｅ ｇ ｅ ｎ ， ｕ ｚｏｇＵ ｉ ｒｉ ｉ ｃ ｄＴ ｃ ｎｌｙ Ｗｕ ａ ０７ ， ｕ ｅ ｈｎ） Ｅ ｇ ｎ ｉ ｖ ｔｏＳ ｅ ａ ｏ ４ ，

测频技术主要有：过零检测技术、相位差分技 术以及现代谱估计技术和 FFT 技术，其中，FFT 是 目前测频的主流手段，但由于非同步采样或非整周
期采样的制约，采用 FFT 会不可避免地产生栅栏效 应，使得测量误差增加，影响后续进一步的解调和 分析。为解决这一问题，文献［5］运用了关于矩形 窗的插值思想；为抑制多谐波相互影响，文献［6］又 提出了加权多点插值方法，有效提高了插值精度； 文献［7］提出的补零技术只能减小频率偏移量的大 小，而且增加数据长度会增加计算量的开销，增加 处理时长；文献［8］提出了基于快速傅里叶变换和 离散时间傅里叶变换相结合的测频计算方法，对确 定的粗频率范围内，通过对频域进行一定程度的放
总第 287 期 2018 年第 5 期
舰船电子工程 Ship舰Ele船ctr电onic子En工gin程eering
Vol. 38 No. 5 37
基 于 汉 宁 窗 的 插 值 FFT 算 法 在 11 号 数 据 链 测 频 中 的应用∗
李晓龙 陈建青 刘 波
（92512 部队 宁波 315113）
Key Words Link11，Hanning window，frequency domain interpolation Class Number TN915.96
1 引言
随着网络中心战的发展，数据链在现代战争中 发挥的作用也日益突出。其中，Link11 作为用于交 换战术数据的数据链，是 C4ISR 系统的关键构成部 分。因此，对 Link11 数据链信号 的 ［1~4］ 检测和识别 有着重要的意义。
示： W (ω) = 0.5 sin(ωN/2)e-jωN 2 {e jω/2/ sin(ω/2) - 0.5
2 关于加窗插值 FFT 的理论分析

第37卷第4期电力科学与工程V ol. 37, No. 4 2021年4月Electric Power Science and Engineering Apr., 2021 doi: 10.3969/j.ISSN.1672-0792.2021.04.004汉宁双窗全相位FFT三谱线插值检测谐波算法彭咏龙，李蕊，马锡浩，李亚斌（华北电力大学电气与电子工程学院，河北保定071003）摘要：在非同步采样情况下，利用快速傅里叶变换（FFT）进行电力系统谐波分析时，会带来频谱泄漏现象和栅栏效应，影响了信号的量测精度。

为此，提出了一种汉宁双窗全相位FFT 三谱线插值检测谐波算法。

该算法原理是：在汉宁双窗全相位FFT分析的基础上，利用基波频点附近的3条相邻谱线幅值作比，计算出频率校正量，并由此估计出谐波信号的幅值；然后，结合全相位FFT分析的相位不变性，将采样点处幅值最大的谱线相位作为信号的初相。

仿真实验表明，与其他插值算法相比，该算法可以更有效地降低谐波参数检测误差，减少白噪声干扰的影响。

关键词：频谱泄漏；全相位FFT；校正算法；三谱线插值中图分类号：TM935 文献标识码：A 文章编号：1672-0792(2021)04-0025-05Harmonic Detection Algorithm Based on Three-spectrum-lineInterpolation of Double Hanning Windows All-phase FFTPENG Yonglong, LI Rui, MA Xihao, LI Yabin(School of Electrical and Electronic Engineering, North China Electric Power University, Baoding 071003, China)Abstract:Harmonic analysis of power systems based on fast Fourier transformation (FFT) produces spectrum leakage and fence effect in non-synchronous sampling, which brings errors to harmonic parameter measurement. To solve this problem, a harmonic detection algorithm based on double Hanning windows all-phase FFT three-spectrum-line interpolation is proposed. The principle of the algorithm is that based on the analysis of double Hanning windows all-phase FFT analysis, the frequency correction is calculated by comparing the three adjacent spectral line amplitudes near the fundamental frequency points and the amplitude of harmonic signal is estimated. And then, based on the invariance of all-phase FFT analysis, the phase of the main spectral line with the largest amplitude at the sampling point is used as the initial phase of the harmonic signal. The simulation experiment shows that compared with other收稿日期：2020-11-21作者简介：彭咏龙（1966—），男，副教授，研究方向为电力电子在电力系统中的应用；李蕊（1995—），女，硕士研究生，研究方向为基于FPGA的电力谐波检测装置。

法等 ， 基于参数模型 的 Ｍｕｉ 算法Ⅲ 、 Ｒ模 型算法 以及最 ｓ ｃ １Ａ ］ 大似然估计算法等 现代谱估 计 的方法＿ ， 以精确估 计正 ２可 ］
弦信 号 频 率 ， 算 法 复 杂 ， 算 量 大 ， 时 处 理 有 难 度 。 而 但 计 实
对 Ｒｉ 算法 和相 位差 算 法进 行 了讨 论 和研 究 ， 做 了改 ｆ ｅ 并
谱 线 对 应 的频 率 之 间 的 偏 差 有 关 ， 以 利 用 相 位 信 息 来 进 可
行频率估 计 ， 相位差法 是通过 对同一 信号进 行不 同长度 或
连续两段的傅里叶变换 ， 先 校正 相位 ， 后 再校 正频率 ， 首 然 从 而 实 现 利 用 Ｆ Ｔ 的 相 位 提 高 频 率 估 计 精 度 。本 文 主 要 Ｆ
１ 引 言
数字信号处理 中对淹没在噪声 中的正弦波 的频 率估计 是一项重要 的内容 ， 速高精度频率估计 已经应用于通 信 、 快 雷达 、 事 、 动 、 军 振 地震 监测 、 纳 以及 其 它 电子 测 试 测 量 技 声
术 等 众 多 领 域 ， 越 来 越 多 的 学 者 的关 注 和 重 视 ， 此 研 引起 因 究 频 率 估 计 具 有 重 要 的理 论 意 义 和 应 用 价 值 。 目前 国 内外 已经 提 出 了不 少 方 法 ， 要 分 为 时 域 、 域 及 时一 分 析 算 主 频 频
精度 。常用 的基 于 Ｆ Ｔ 的插 值 算 法有 双线 幅 度 Ｒｉ Ｆ ｆ ｅ算
法 ［ ｊ相 位 差 法 ｌ ］ 结 合 幅 度 与 相 位 的 双 线 Ｑｕｎ ３ 、 ６ 、 ｉｎ算
法ｌ ＿ ｌ 。其 中， ｉ 算法是借 助第二谱 线与最 大谱 线 的幅 等 Ｒｆ ｅ 度 比值来估 计 信 号 的 实 际频 率在 两 条 谱线 之 间 的位 置 。 Ｆ Ｔ最大值处 的相位与 信号 的实际频 率和 Ｆ Ｔ幅度最 大 Ｆ Ｆ

1.#include "stdafx.h"2.#include "Fft.h"3.#include "math.h"4.5.const double MOD_MAX = 65535.0;6.const double DISP_MAX = 1.0/255.0;7.extern FILE *fp;8.9.10.11.12.13.//FFT运算必须参数14.int fft_point,fft_order,fft_divide,fft_window,fft_scale;15.bool fft_cover;16.float filter[7];//FIR滤波参数17.18.extern FILE *fpIandQ;19.extern bool m_bIqWrite;20.double prFilter[256],piFilter[256];21.22.23./**************************************************************************************24. 0 RectangleWindow矩形窗25. FFT变换结果为对称型，矩形窗是使信号突然截断，旁瓣会很大，且衰减较慢，旁瓣的第一个负26.峰值为主瓣的21％，第一个正峰值为主瓣的12.6％，第二个负峰值为主瓣的9％，效果一般，泄27.漏较大。

28.**************************************************************************************/29.double WINAPI RectangleWindow(int t)30.{31.double wt;32. wt=1.0;33.return wt;34.}35.36./**************************************************************************************37. 1 TriangleWindow三角窗,也称费杰（Fejer）窗，Bartlett38.**************************************************************************************/39.double WINAPI TriangleWindow(int t)40.{41.42.double wt;43. wt=1-t/fft_point;44.return wt;45.}46.47./**************************************************************************************48. 2 HanningWindwo汉宁窗，即升余弦窗49.**************************************************************************************/50.double WINAPI HanningWindow(int t)51.{52.double wt;53. wt=(1-cos(2*PI*t/fft_point))/2;54.return wt;55.}56.57./**************************************************************************************58. 3 HammingWindow海明窗，即改进的升余弦窗59.**************************************************************************************/60.double WINAPI HammingWindow(int t)61.{62.double wt;63. wt=0.54-0.46*cos(2*PI*t/fft_point);64.return wt;65.}66.67./**************************************************************************************68. 4 BlackmanWindow布来克曼窗，即三阶升余弦窗69.**************************************************************************************/70.double WINAPI BlackmanWindow(int t)71.{72.double wt;73. wt=0.42-0.5*cos(2*PI*t/fft_point)+0.08*cos(4*PI*t/fft_point);74.return wt;75.}76.***********78. 5 CosgradeWindow余弦坡度窗79.**************************************************************************************/80.double WINAPI CosgradeWindow(int t)81.{82.double wt;83.if(t= int(4*fft_point/5) )84. wt=1.0;85.else86. wt=(1+cos(5*PI*t/fft_point))/2;87.return wt;88.}89.90.91./**************************************************************************************92. 6 ParzenWindow帕曾窗93.**************************************************************************************/94.double WINAPI ParzenWindow(int t)95.{96.double wt;97.if(t= int(fft_point/2) )98. wt=1-6*pow(t/fft_point,2)+6*pow(t/fft_point,3);99.else100. wt=2*pow(1-t/fft_point,3);101.return wt;102.}103.104./*************************************************************************** ***********105. 7 ExponentWindow指数窗106.**************************************************************************** **********/107.d ouble WINAPI ExponentWindow(int t)108.{109.double wt,alf;110. alf=0.0001;111. wt=exp((-1.0)*alf*t);112.return wt;113.}114.***********116. 8 GaussWindow高斯窗117.**************************************************************************** **********/118.d ouble WINAPI GaussWindow(int t)119.{120.double wt,alf;121. alf=0.0001;122. wt=exp(-1*alf*t*t);123.return wt;124.}125.126./*************************************************************************** ***********127. 9 KaiserWindow凯泽窗,128.凯泽在1966（1974）发现，利用第一类零阶修正（变形）贝赛尔函数可以构成一种近似最佳的窗129.最主要参数：beat－可同时调整主瓣宽度和旁瓣，beat越大，窗越窄，频谱旁瓣越小，但主瓣相应增加130. beta=0相当与矩形窗131.**************************************************************************** **********/132.d ouble WINAPI KaiserWindow(int t)133.{134.double wt,alfa,beta;135. alfa=(fft_point-1)/2;136. beta=0;137.//Bessel零阶第一类贝塞尔函数138. wt=Jim_Bessel0_R(0,beta*sqrt(1-pow((t-alfa)/alfa,2)))/Jim_Bessel0_R(0,be ta);139.return wt;140.}141.142./*************************************************************************** ***********143. 10 ChebWindow契比雪夫窗144.**************************************************************************** **********/145.d ouble WINAPI ChebWindow(int t)146.{147.double wt;148. wt=1.0;149.return wt;150.}151.152./*************************************************************************** ***********153. 11 BartlettWindow巴特利特窗154.**************************************************************************** **********/155.d ouble WINAPI BartlettWindow(int t)156.{157.double wt;158. wt=1.0;159.return wt;160.}161.162.b ool WINAPI fftInit(int point,int order,int divide,bool cover,int scale) 163.{164.165. fft_point=point;166. fft_order=order;167. fft_divide=divide;168. fft_cover=cover;169. fft_scale=scale;170.//采用8阶的FIR滤波171. Jim_FirFilter(7,3,50,7800/2.0,7800.0,1,filter);172.return TRUE;173.}174.175.176.B OOL WINAPI fftEnd()177.{178.return TRUE;179.}180.181.182.183.184.//傅立叶变换185.B OOL WINAPI fftTransform(VOID *inBuf,VOID *outBuf,int windows,int scale,int nProbe,int nWork)186.{187.BYTE *srcPtr = (BYTE*)inBuf;//数据源1024字节188.BYTE *dstPtr = (BYTE*)outBuf;//数据果256字节189.//double *testPtr =(double *)testBuf;190.WORD Ivalue,Qvalue;191.//WORD IvalueRev,QvalueRev;192. unsigned char flagQ,flagI;193.//double alfa,e;194.int i,j,Iorg=0,Qorg=0;195.double mod = 0;196.if(nProbe==LEFT_PROBE)197. {198. flagQ=0x30,flagI=0x20;199. }200.else201. {202. flagQ=0x10,flagI=0x0;203. }204.205.//每次取4个字节数据，分离I/Q分量，判断db(12)，0表示I分量，1表示Q分量206.if (((srcPtr[1]&0xf0) ==flagQ) && ((srcPtr[3]&0xf0) ==flagI)) 207. {208. Iorg = 2;209. Qorg = 0;210. }211.if (((srcPtr[1]&0xf0) ==flagI) && ((srcPtr[3]&0xf0) ==flagQ)) 212. {213. Iorg = 0;214. Qorg = 2;215. }216.double pr[256],pi[256],w;217.218.219.//计算自相关函数和互相关函数220.//double Rii=0.0,Rqq=0.0,Riq=0.0;221.222.for ( i=0;i<fft_point;i++ )223. {224.225.//更改存储的偏移地址分离I/Q，226. Ivalue = *((short*)(srcPtr + i*4 + Iorg));227. Qvalue = *((short*)(srcPtr + i*4 + Qorg));228.//~位非运算，低12位位有效数据229. pr[i]=double(Ivalue & 0xfff);230. pi[i]=double(Qvalue & 0xfff);231.232./*Rii=Rii+double(pr[i]*pr[i]);233. Rqq=Rqq+double(pi[i]*pi[i]);234. Riq=Riq+double(pr[i]*pi[i]);*/235. }236.237./*Rii=Rii/fft_point;238. Rqq=Rqq/fft_point;239. Riq=Riq/fft_point;240.241. alfa=asin(Riq/sqrt(Rii*Rqq));242. e=sqrt(Rqq/Rii)-1;*/243.244.for ( i=0;i<fft_point;i++ )245. {246.//防止谱泄漏，进行加窗处理247.switch(windows)248. {249.case WND_RECTANGLE: w=RectangleWindow(i); break; 250.case WND_TRIANGLE: w=TriangleWindow(i); break; 251.case WND_HANNING: w=HanningWindow(i); break; 252.case WND_HAMMING: w=HammingWindow(i); break; 253.case WND_BLACKMAN: w=BlackmanWindow(i); break; 254.case WND_COSGRADE: w=CosgradeWindow(i); break; 255.case WND_PARZEN: w=ParzenWindow(i); break; 256.case WND_EXPONENT: w=ExponentWindow(i); break; 257.case WND_GAUSS: w=GaussWindow(i); break; 258.case WND_KAISER: w=KaiserWindow(i); break; 259.case WND_CHEB: w=ChebWindow(i); break; 260.case WND_BARTLETT: w=BartlettWindow(i); break; 261.default: w=RectangleWindow(i); break; 262. }263.264./*pr[i]=((1+e)*cos(alfa)*pr[i])/((1+e)*cos(alfa))*w; 265. pi[i]=(pi[i]-(1+e)*sin(alfa)*pr[i])/((1+e)*cos(alfa))*w;*/ 266.267. pr[i]=pr[i]*w;268. pi[i]=pi[i]*w;269. prFilter[i]=piFilter[i]=0.0;270.for(j=0;j<7;j++)271. {272.if(i+j>255)273. {274. prFilter[i]=prFilter[i];275. piFilter[i]=piFilter[i];276. }277.else278. {279. prFilter[i]=prFilter[i]+pr[i+j]*w*filter[j];280. piFilter[i]=piFilter[i]+pi[i+j]*w*filter[j];281. }282. }283.//fprintf(fp,"%d,%f\n",i,filter[i]);284.//pr[i]=pr[i]*filter[i];285.//pi[i]=pi[i]*filter[i];286. }287.//fclose(fpIandQ);288.289.//fft处理290.//Jim_FFT(pr,pi,fft_point,fft_order,0,1);291. Jim_FFT(prFilter,piFilter,fft_point,fft_order,0,1);292./*for(int t=0;t<256;t++)293. {294. fprintf(fpIandQ,"%d\n",int(prFilter[t]));295. }*/296.//归一化297.//Jim_Unitary(pr,dstPtr,scale,nProbe);298. Jim_Unitary(prFilter,dstPtr,scale,nWork);299.300.301.return TRUE;302.}303.304./*************************************************************************** ***********305. I/Q相位校正算法 by Jim Fang at 2007306.**************************************************************************** **********/307.v oid WINAPI Jim_Pharev(void *pIvalue,void *pQvalue,void *pRev)308.{309.double *iPtr =(double *)pIvalue;310.double *qPtr =(double *)pQvalue;311.double *rPtr =(double *)pRev;312.//计算自相关函数和互相关函数313.double Rii=0.0,Rqq=0.0,Riq=0.0;314.int i;315.for ( i=0;i<fft_point;i++)316. {317. Rii=Rii+double(iPtr[i]*iPtr[i]);318. Rqq=Rqq+double(qPtr[i]*qPtr[i]);319. Riq=Riq+double(iPtr[i]*qPtr[i]);320. }321. rPtr[0]=sqrt(Rqq/Rii)-1;322. rPtr[1]=asin(Riq/sqrt(Rii*Rqq));323.}324.325./*************************************************************************** ***********326.能量归一化算法 by Jim Fang at 2007327. FFT的源数据为WORD型，范围在±32768之间，用于显示能量的色阶范围为0-255328. 1.消除奇异点329. 2.找出峰值330. 3.归一化331.**************************************************************************** **********/332.v oid WINAPI Jim_Unitary(void *pSrcData,BYTE *pDstData,int nScale,int nWork)333.{334.//平滑处理后取最大值作为255，进行归一化。