專利名稱：快速高精度信號頻譜分析方法
技術(shù)領(lǐng)域：
本發(fā)明是關(guān)于一種對被測電信號進行頻譜分析的方法，更具體地說它是關(guān)于一種具有抗混效果的快速富里葉頻譜分析方法。
當前，數(shù)字計算機與數(shù)字信號處理器(DSP)等的迅速發(fā)展與普及，使得模擬信號也被離散成數(shù)字信號來處理。然而從信號的復原與信號的頻譜分析角度考慮，模擬信號用離散信號代替后會出現(xiàn)以下原理誤差。其一是非時限信號用有限時域的離散信號代替時產(chǎn)生的截尾誤差(Truncatederrors)，對周期信號這種誤差又表現(xiàn)為采樣周期與信號同期不同步時出現(xiàn)的泄漏誤差(Leakage errors);其二是非帶限信號時，以及帶限信號但采樣頻率低于香農(nóng)(Shannon)采樣定理要求時出現(xiàn)的頻譜混迭誤差(Aliasing errors);其三是離散信號頻譜周期性重復，不管原模擬信號包含多少頻率分量，采樣N點均只能確定N/2個帶有上述誤差的頻譜系數(shù)，并且頻率越高的分量，頻譜系數(shù)的誤差越大。
泄漏誤差已找到一些較有效的解決辦法，比如加窗技術(shù)與加窗插值技術(shù)等?；斓`差雖然已深入研究80余年，在各種介紹FFT原理的文獻中都被涉及到，但還只停留在如何確定誤差規(guī)模階段;如何消除這種誤差，特別是如何在信號快速處理過程中消除這種誤差尚待解決。
頻譜分析是信號處理的核心，電子系統(tǒng)的諧波分析即頻譜分析之一種。分析連續(xù)信號頻譜時，目前普遍采用快速富里葉變換(FFT)?，F(xiàn)有各種版本的FFT都只是離散富里葉變換(DFT)的一種快速算法，屬數(shù)字處理技術(shù)之列，存在著本發(fā)明要討論的頻譜混迭誤差與頻譜系數(shù)周期性重復的上述通病?；斓`差嚴重程度隨信號形式而異。以一周期內(nèi)的表達式為x(t)＝(2/T)t-1
的鋸齒波為例(見附表1)，在采樣周期與信號周期同步因而沒有泄漏誤差的情況下，一周期采樣64點時，現(xiàn)有FFT算出的第31次諧波的正弦系數(shù)只有理論值的0.074748，一周期采樣128點時，31次諧波的正弦系數(shù)也只有理論值的0.79914;余弦系數(shù)(理論值全為零)誤差更大?；斓`差的影響由此可見一斑。
為了消除混迭誤差，現(xiàn)有頻譜分析儀采取的對策是，用低通濾波器將信號中高于1/2采樣頻率的高頻分量濾除掉。這種方法雖能減少，所分析出的N/2個頻譜系數(shù)中的混迭誤差，但卻帶來其它誤差。原因是低通濾波器的幅頻特性不可能是矩形，相頻特性不可能是直線;信號通過后，高頻分量濾掉了，低頻分量的幅度與相位卻產(chǎn)生了失真，使得分析結(jié)果仍與真實情況不同。特別是低通濾波器帶來的相位誤差，在許多要求精確知道諧波相位的場合是不允許的。因此采用低通濾波器也不能很好解決現(xiàn)有FFT算法存在的精度不高，所分析的頻譜范圍不寬的問題。
面對這一情況，一種保留信號的高頻分量仍能壓制頻譜混迭誤差，從而可大大提高頻譜分析精度與擴大所分析的頻譜范圍，并能縮短分析時間的高精度快速頻譜分析方法，不僅對工程應用有重大意義，對數(shù)字處理技術(shù)與信號分析理論也有重大意義。
本發(fā)明的目的就是要提供一種上述的快速高精度頻譜分析方法。
下面闡述一下本發(fā)明快速高精度頻譜分析方法的數(shù)字原理。
滿足Dirichlet條件的函數(shù)x(t)，在[o，T]區(qū)間上可用Fourier級數(shù)表示為
式中
AK稱為K次諧波的復振幅，AK是其幅值，ψK是其相位，aK，bK分別為K次諧波的余弦系數(shù)與正弦系數(shù)。所謂頻譜分析，就是求出各次諧波對應的AK。
設△為很小的時間間隔，T＝N·△，x(t)在tn＝n△時的值為xn。用各個xn代替上式中的x(t)，即用求和代替積分時，就得到了現(xiàn)有快速頻譜分析方法FFT計算AK的公式
采樣定理從信號理論角度分析了AKF與AK的差別，從數(shù)學角度考慮，AKF是求和得到的，真正的AK應由積分得到，故用FFT方法得到的信號頻譜與真實頻譜之間必存在原理誤差。
為求出真正的AK，不能采用現(xiàn)有FFT采用的簡單求和方式，本發(fā)明采用精確計算(2)式所給的積分。
計算可分段進行
式中
一種可供選擇的計算AKn的方法是利用歐拉數(shù)值積分公式，參見安德列、安戈著、陸志剛等譯，人民郵電出版社出版的“電工、電信工程師數(shù)學(下)”，第十章。該公式計算積分時多次實行分部積分法，直至后面的分部積分結(jié)果可忽略。此法直接運用于(5)式時效果不佳，原因是被積函數(shù)中有周期因子exp(-jzkπt/T)，無論分部積分多少次，這個因子總原樣存在，無法判定后面部分是否可忽略。針對(5)式這樣一種特殊被積函數(shù)，我們對歐拉數(shù)值積分方法作了改進。我們先將被積函數(shù)中的xn(t)部分在[tn-△，tn+△]區(qū)間展開成臺勞級數(shù)，由于△很小，xn(t)可用它在該區(qū)間展開式的前m+1項來逼近Xn(t)≈Σm=0Mam(t-tn)m(6)]]>M的取值視△的取值而定，△大則逼近xn(t)所需的M值高，△小則M值低。將(6)式代入(5)式，這時再多次使用分部積分法，直至被積函數(shù)變?yōu)榱銥橹?，結(jié)果得
式中xn(n2l)(t)，xn(2l+1)(t)分別為xn(t)的2l階與(2l+1)階導數(shù)。現(xiàn)有FFT的頻譜系數(shù)計算方法，可看成是(6)式在M＝0時的特例，即用x(t)在tn處的采樣值代替[tn-△，tn+△]區(qū)間的xn(t)，因而帶來混迭誤差。為消除混迭誤差，必須考慮xn(t)展開式中的高式次項。一般△<<T，因此取xn(t)在該區(qū)間的級數(shù)展開式的前三項即可有相當高的精度。將(6)式中的各個系數(shù)am用t＝(n-1)△，n△，(n+1)△時的采樣值xn-1，xn，xn+1的函數(shù)來表示，并考慮截尾誤差，然后代入(7)式，計算整理得
式中，ZK= UK+ jVK(9)Z＊K為ZK的共軛復數(shù)，VK、UK與RK則是只與比值K/N有關(guān)的實常數(shù)，其值隨K的增加而單調(diào)下降。將(8)式代入(4)式，最終得
式中WN= e-J(2π)/(N) (11)F＝(X0-XN)/N (12)(10)式就是本文高精度頻譜分析的基本公式。
(10)式中上述的實常數(shù)UK、VK、RK的數(shù)學式如下πUX=3+2cos2(2KπN)2(2KπN)2-Sin2(2KπN)(2KπN)3(13)]]>VX=N2Kπ+Sin2(2KπN)2(2Kπ/N)2-1-cos2(2kπ/N)(2kπ/N)2(14)]]>
RX=4|Sin(2kπ/N)(2Kπ/N)2-COS(2Kπ/N)(2Kπ/N)2|(15)]]>(10)式所示的具有抗混(Anti-Aliasing)效果的AK精確表達式如不能實現(xiàn)快速計算，本文方法將難以推廣應用;而按(10)式快速計算AK時，如能盡可能利用現(xiàn)有FFT程序，則本文方法更易于推廣普及。下面我們來尋找滿足這種要求的算法，并將這種算法命名為FAFT(Fast Anti-Aliasing Fourier Transform)。
(10)式中的UK、RK、VK都是與采樣值無關(guān)的常數(shù)，可預先算好儲存?zhèn)溆?F也只與首尾兩點的采樣值有關(guān)，容易計算。關(guān)鍵是如何實現(xiàn)
的快速計算，為此我們令
k于是
為
＝2UKAeven(K)+RKAodd(K)WW-ZWF (18)這時如能同時實現(xiàn)Aeven(K)與Aodd(K)的快速計算就等于實現(xiàn)了AK的快速計算。由(16)，(17)兩式可知，Aeven(K)與Aodd(K)分別為N/2個偶次采樣數(shù)據(jù)與奇次采樣數(shù)據(jù)對應的復頻譜系數(shù)，為計算它們須先將采樣數(shù)據(jù)按奇、偶順序重新排列。這并非本文方法的額外要求，現(xiàn)有FFT也須如此重新排列數(shù)據(jù)時域抽取法(Decimation in Time)開始時進行，頻域抽取法(Decimation in Frequency)結(jié)尾時進行。所以這一步可直接利用現(xiàn)有FFT的有關(guān)程序?qū)崿F(xiàn)。
計算 even(K)與 odd(K)也可借用現(xiàn)有FFT程序，比如時域抽取法的Cooley-Tukey算法。仔細分析這種算法可知，對于N＝2S個被分析數(shù)據(jù)，該算法共循環(huán)S步，其中第(S-1)步變換出的數(shù)據(jù)，正好是這里的 even(k)與 odd(k)。所以利用現(xiàn)有FFT算法即可實現(xiàn) even(K)與 odd(k)的快速計算。第(S-1)步后，將現(xiàn)有FFT程序稍加改動，即不繼續(xù)計算 even(K)WKN與 odd(K)WKN，而是只計算后者，然后按(15)式要求合成 。這樣，借用幾乎沒作什么改動的現(xiàn)有FFT程序，就可一舉實現(xiàn)高精度AK的快速計算。
下面將結(jié)合附圖對本發(fā)明進行詳細說明。


圖1是實現(xiàn)本發(fā)明方法的頻譜分析儀的方框示意圖;
圖2是本發(fā)明信號頻譜分析方法的流程圖;
圖3是第一實驗的信號波形;
圖4是第二實驗的信號波形。
為了便于理解，首先對實現(xiàn)本發(fā)明方法的頻譜分析儀的方框圖作一介紹(參見圖1)。它主要包括數(shù)據(jù)采集電路(1)、存儲器(2)、中央處理單元CPU(3)、顯示器(4)、打印機(5)以及鍵盤(6)等。數(shù)據(jù)采集電路首先將要分析的信號(模擬形式)拾取，采樣并轉(zhuǎn)換成數(shù)字信號，該數(shù)字信號通過總線(7)送入存儲器(2)，中央處理單元CPU(3)將存儲器(2)中的數(shù)據(jù)進行處理并計算，最終得出各次諧波分量的系數(shù)AK值，并顯示在顯示器上或通過打印機打印出來。
下面將結(jié)合圖1和圖2對本發(fā)明的頻譜分析方法作一祥細介紹。
按照本發(fā)明的頻譜分析方法，它包括如下步驟1、初始數(shù)據(jù)采集數(shù)據(jù)采集電路(1)對輸入的被測信號進行N個點的采樣，并將采樣點的值變成數(shù)字信號存入存儲器(2)中。設該N個采樣值為X0，X1，X2……XN其中采樣點數(shù)N可以根據(jù)測量的精度要求而定，比如N為4、8、16……，實際操作時可以通過鍵盤來設定N的值。N值定了以后，按照N＝2S則S的值自然就定了，S為現(xiàn)有FFT中的計算步驟。
根據(jù)所存儲的采樣值計算FF＝(X0-XN)/N2、利用現(xiàn)有的FFT程序作奇偶排列即將存儲在存儲器(2)中的初始采樣數(shù)據(jù)按奇數(shù)和偶數(shù)重新排列，以便于以后的計算，X0，X2，X4……XN-2X1，X3，X5……XN-13、執(zhí)行現(xiàn)有FFT方法前(S-1)步，計算
4、完成現(xiàn)有FFT方法第S步的一半計算內(nèi)容，求出
5、合成 k按照下式求出 k: k＝2UKAeven(K)+RKA′odd(K)-ZKF
其中UK、RK、ZK是事先分別按照公式(13)、(15)和(9)計算好并存在存儲器中的。
6、顯示或打印將上面求出來的K次諧波系數(shù)AK顯示在顯示器(4)上或通過打印機(5)打印出來。
本發(fā)明的快速高精度信號頻譜分析方法(FAFT)用三項多項式逼近各區(qū)段的X3(t)，精確算出了Akn對應的積分值，精度必然比FFT高;這樣，采樣間隔取大一些仍能比FFT精度高。對于同一個被采樣信號，采樣間隔大就意味著采樣點數(shù)N少，而分析過程中所需的復數(shù)乘、加運算次數(shù)M與采樣點數(shù)N之間有以下關(guān)系M＝Nlog2N，N小即意味著計算時間短，因此FAFT與FFT相比，不僅可大大提高精度，還可縮短計算時間。另外，F(xiàn)AFT不是用離散信號頻譜代替連續(xù)信號頻譜，不受采樣定理限制。對帶限信號，即使用低于采樣定理要求的速率采樣，仍能有很高的精度。在硬件采樣速率已定的前題下，這等于擴大了可分析的頻譜范圍。而在精度要求一定的情況下，與FFT相比，F(xiàn)AFT可采用采樣速率低的采樣電路與位數(shù)低好A/D變換器，這就可降低硬件成本。
FFT的另一缺點，即采樣N點只能確定N/2個頻譜系數(shù)的缺點，是這一方法算出的
這一周期性重復性質(zhì)造成的。FAFT無此弊病。
由
even(K)，
odd(K)與WKN的遞推性質(zhì)
even(K+mN/2)＝
even(K) (19)
odd(K+mN/2)＝
even(K) (20)WK+mNN＝WNK(21)WK+N(m+1/2)N＝-WKN(22)
可得FAFT法算出的AK的遞推公式為
由于UK，RK，ZK無周期重復性質(zhì)，故FAFT法算出的AK+mN/2與AK之間不存在周期性重復關(guān)系，即利用N個采樣值可算出N/2個以上的頻譜系數(shù)。
FAFT方法可直接利用現(xiàn)有信號分析儀實現(xiàn)，也可另行設計更能充分利用FAFT優(yōu)點的硬件來實現(xiàn)。與現(xiàn)有信號分析儀相比，采用FAFT后不僅可降低對采樣速率與A/D精度的要求，還可省去數(shù)據(jù)采集電路中的抗混濾波器及有關(guān)部分，使硬件成本降低，結(jié)構(gòu)簡化。
為了證實FAFT的上述優(yōu)點，對已知頻譜系的兩種常見波形-鋸齒波X1(t)與余弦全波整流波形X2(t)，分別用FAFT與標準FFT[10]作了頻譜分析。X1(t)，的波形分別見圖3、4，它們在一周期內(nèi)的表達式與理論頻譜系數(shù)分別為X1(t)=(2/T)t·1=·2Σk=1∞1KπSinK2πTt=·2Σk=1∞bkSinK2πTt---(25)]]>式中，bK＝+ 1/(kπ) (26)X2(t)=|Cos2πTt|=2π+Σk=1∞(-1)k-12π(4K2-1)Cos2K2πTt]]>=a0+2Σk=1∞a2kCos2K2πTt---(27)]]>
式中，a0＝ 2/(π) (28)a2k=(-1)k-12π(4K2-1)(29)]]>分析結(jié)果見表1-表4。表中f1＝1/T為信號基波頻率，fS＝1/△為采樣頻率，K為諧波次數(shù)，fK＝Kf1為K次諧波頻率，aK，bK為分析值，aK，bK為理論值。為節(jié)省篇幅，每種比較只列出了一種波形的結(jié)果。
表1是同樣的采樣速率與A/D精度時兩種方法分析精度對比?？煽闯?，不僅fS相同時FAFT的精度遠比FFT的精度高，即使FAFT的fSFFT的fS低時，F(xiàn)AFT的精度也遠比FFT的精度高。
表2是同樣的A/D精度，不同的采樣速率時，兩種方法的分析精度對比。可看出，對波形1即使FAFT的fS有FFT的fS的1/64，分析出的前31次諧波系數(shù)的精度仍比FFT的精度高。
表3是同樣的分析精度時兩種方法所需的采樣速度與A/D精度對比?？煽闯?，對波形2在分析出的前62次諧波系數(shù)精度基本相同的情況下，F(xiàn)FT需采用每周期采樣256點的采樣電路與18位的A/D時，F(xiàn)AFT只需采用每周期采樣128點的采樣電路與10位的A/D。
表4是不同的采樣速度時，F(xiàn)AFT分析出的fK＝50fs、120fs、250fs的K次諧波的諧波系數(shù)及精度。這不僅說明FAFT可擴大諧波分析范圍，還說明此法的確不受采樣定理限制。按采樣定理，欲分析第3998次諧波，采樣頻率至少應為基頻的7996倍，即fS≥7996f1;FAFT方法，fS＝16f1即可以相當高的精度分析出3998次諧波及更高次數(shù)諧波的諧波系數(shù)。本例還證明與FFT相比，F(xiàn)AFT可大大提高分析速度。比如欲分析0～3998次諧波，F(xiàn)FT需處理7996個以上的采樣數(shù)據(jù)，F(xiàn)AFT只處理(16+1)個數(shù)據(jù)即可;這時兩種方法所需計算時間長短不言自明。
由于fS高則一周期內(nèi)采樣數(shù)據(jù)總數(shù)N大，所需的分析時間長，故FAFT方法可提高分析速度的優(yōu)點，由表1-表3也可得到證明。
權(quán)利要求
1.一種快速高精度信號頻譜分析方法，它包括下列步驟(1)初始數(shù)據(jù)采集即通過數(shù)據(jù)采集電路對輸入的被測信號進行N個點的采樣X0，X1，X2……XN，并將它們變成數(shù)字信號存入存儲器中，再計算F值[F=(Xo-Xx)/N](2)利用現(xiàn)有的FFT程序?qū)ι鲜鯪個采樣值作奇偶排列，即X0，X2，X4……XN-2X1，X3，X5……XN-1(3)按照現(xiàn)有FFT方法前(S-1)步計算Aeven(k)
(4)進行現(xiàn)有FFT方法第S步的一半計算內(nèi)容，求出A′odd(k)
其特征在于該方法還包括下述步驟(5)利用事先存儲的實常數(shù)Uk、Rk和Zk按下式合成Akf0=1/2πL0C]]>(6)顯示或打印f1=1/2πL1C]]>數(shù)據(jù)。
全文摘要
本發(fā)明在嚴格的理論分析基礎(chǔ)上，提出了一種新的快速高精度信號頻譜分析方法。它主要包括下列步驟初始數(shù)據(jù)采集；奇偶排列；按現(xiàn)有快速富利葉變換方法計算Aeven(k)和A＇odd(k)；利用事先存儲的實常數(shù)Uk、Rk和Zk合成Ak；將K次諧波的系數(shù)Ak顯示或打印。本發(fā)明不僅可提高頻譜分析精度，擴展頻譜分析范圍，縮短分析時間，還可降低硬件性能要求與簡化硬件結(jié)構(gòu)。
文檔編號G01R23/16GK1062793SQ9010602
公開日1992年7月15日 申請日期1990年12月25日 優(yōu)先權(quán)日1990年12月25日
發(fā)明者陳祥訓 申請人:能源部電力科學研究院