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


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