楊 超 李 波,2 胡緒權(quán),3 付志紅

基于解析插值離散時間傅里葉變換的精確頻率估計

楊 超1李 波1,2胡緒權(quán)1,3付志紅1

（1. 輸配電裝備及系統(tǒng)安全與新技術(shù)國家重點實驗室（重慶大學(xué)） 重慶 400044 2. 云南電網(wǎng)有限責(zé)任公司電力科學(xué)研究院 昆明 650217 3. 重慶璀陸探測技術(shù)有限公司 重慶 402660）

短時波動周期下，來自負(fù)頻率的頻譜泄漏是頻域插值算法估計誤差的主要來源。為消除頻譜泄漏影響以實現(xiàn)任意波動周期下實正弦信號頻率的精確估計，該文提出一種基于廣義離散時間傅里葉變換（DTFT）插值解析的頻率估計算法。該算法利用位于頻譜主瓣范圍內(nèi)的任意DTFT譜線構(gòu)建對應(yīng)實/虛部的線性方程組，進而通過方程組求解得到基于解析插值DTFT的無偏頻率估計公式。此外，考慮實際中加性白噪聲的影響，根據(jù)頻率估計的統(tǒng)計特性分析給出一種簡單且有效的迭代方法，從而確保頻率估計的方均誤差（MSE）在任意波動周期下能夠一致逼近無偏頻率估計的克拉美羅下屆（CRLB）。仿真和實驗結(jié)果驗證了所提算法優(yōu)于現(xiàn)有的DFT插值算法。

短時波動周期 頻譜泄漏 解析插值DTFT 迭代方法 克拉美羅下屆（CRLB）

0 引言

正弦信號的頻率估計是一個經(jīng)典且重要的信號處理問題，在科學(xué)和工程領(lǐng)域具有極為廣泛的應(yīng)用。在電力系統(tǒng)中，頻率作為最重要的電能質(zhì)量參數(shù)之一，其變化反映了發(fā)電和負(fù)載之間的動態(tài)平衡[1-3]，通過頻率及其變化的準(zhǔn)確估計可以實現(xiàn)系統(tǒng)的動態(tài)過程監(jiān)測。此外，在計量基準(zhǔn)測試和標(biāo)準(zhǔn)溯源等特定應(yīng)用中，對頻率估計精度往往具有更高的要求。

頻率估計算法總體上可分為時域和頻域兩類，頻域算法主要基于由快速離散傅里葉變換（Fast Discrete Fourier Transform,FDFT）實現(xiàn)的插值離散傅里葉變換（Interpolated Discrete Fourier Transform, IpDFT），在計算復(fù)雜度方面相對于通常需要矩陣求逆的時域方法更具優(yōu)勢（時域為(3)，頻域為(log2)[4]），因而在實時頻率估計應(yīng)用中表現(xiàn)出特別的吸引力。但是，過去大多數(shù)針對IpDFT頻率估計的研究都集中在單頻解析信號上[5-9]，對于實際中應(yīng)用更為廣泛的實正弦信號，受共軛負(fù)頻率的頻譜長泄漏影響，尤其是在被測信號對應(yīng)波動周期（Cycles in Record, CiR）較小時，IpDFT存在嚴(yán)重的估計偏差。

短時CiR頻譜泄漏問題已引起國內(nèi)外學(xué)者的高度關(guān)注并開展了大量研究，提出多種改進IpDFT算法，總體上可分為直接方法和迭代方法。直接方法力求在計及負(fù)頻率頻譜長泄漏基礎(chǔ)上直接得到頻率估計的解析解。文獻[10]提出了一種基于矩形窗的3點IpDFT（3IpDFT），利用變換構(gòu)建對應(yīng)譜線的線性方程組，進而求解得到頻率估計的解析解，但仿真結(jié)果顯示，即便在最優(yōu)情況下算法方均誤差（Mean Square Error, MSE）仍然無法逼近無偏頻率估計的克拉美羅下界（Cramer-Rao Lower Bound, CRLB）。文獻[11-12]在文獻[13]解析泄漏補償方法基礎(chǔ)上，提出了利用信號補零和DFT系數(shù)變換的2點IpDFT（2IpDFT），但補零操作和繁瑣的譜線選擇判據(jù)使得算法計算復(fù)雜度相對較高。文獻[14-15]基于最大旁瓣衰減窗離散等間隔頻譜近似線性比例特性，以求解頻譜線性方程組的方式得到了3點IpDFT（3IpDFT），但鑒于窗函數(shù)的數(shù)學(xué)復(fù)雜性增加[16]，事實上在高階情況下很難得到頻率估計的解析解，因而無法進一步消除窗函數(shù)頻譜近似導(dǎo)致的估計偏差。此外，相對矩形窗而言，采用具有更高等效噪聲帶寬（Equivalent Noise Band Width, ENBW）[17]的最大旁瓣衰減窗，算法MSE在噪聲擾動下將無法逼近CRLB。

迭代方法則常以適用于單頻解析信號的IpDFT為基礎(chǔ)，力求通過迭代不斷消除頻譜泄漏干擾，從而確保算法最終收斂于頻率估計的解析解。文獻[18]在最大旁瓣衰減窗頻譜近似基礎(chǔ)上，提出了一種基于頻譜系數(shù)加權(quán)的單步迭代IpDFT（3IpDFT），但以經(jīng)典IpDFT[19]作為其起始解析方式，將導(dǎo)致短時CiR下算法MSE中存在明顯的估計偏差。文獻[20]對AM算法[21]進行了有效改進，提出了一種利用位于DFT譜線中間的離散時間傅里葉變化（Discrete Time Fourier Transform, DTFT）系數(shù)插值的多步迭代算法（AM）。其基本思想與2IpDFT類似，即盡量采用幅度最大的頻譜系數(shù)以降低噪聲干擾影響，進而通過不斷重構(gòu)并消除頻譜中負(fù)頻率泄漏分量以逐漸降低頻率估計偏差，使得算法MSE最終一致收斂逼近于CRLB，是現(xiàn)有頻域估計算法中精度最高的[22]。但由于其核心兩點插值方法為有偏頻率估計，在估計偏差占主導(dǎo)情況下需要多次計算DTFT系數(shù)并進行迭代，這無疑大大增加了算法的計算成本。此外，采用固定0.5D（D為頻率分辨率）的DTFT系數(shù)還使得算法在CiR＜1時性能將急劇下降。

基于AM算法的基本思想，本文提出了一種基于廣義DTFT插值的精確頻率估計算法。一方面，通過準(zhǔn)確解析由任意位置的廣義DTFT譜線構(gòu)成的線性方程組，可以確保算法在任意CiR下得到無偏頻率估計；另一方面，通過統(tǒng)計特性分析導(dǎo)出用以優(yōu)化算法性能的簡單迭代方法，可以最大程度地抑制背景噪聲干擾以確保算法MSE能夠一致逼近CRLB，從而實現(xiàn)任意CiR下頻率的精確估計。

1 頻率估計的解析解

1.1 基本原理

假設(shè)以采樣頻率s進行等間隔均勻離散采樣，則點采樣序列()的數(shù)學(xué)模型可表示為

式中，=0, 1,…,-1?N；、和分別為信號()的頻率、幅值和初始相位；()為具有零均值和方差2的加性高斯白噪聲；/s//為頻域歸一化頻率，在時域中對應(yīng)于被測信號的波動周期個數(shù)，在頻域中則表示在頻譜中所處位置。據(jù)奈奎斯特采樣定律可知，max＜0.5，且無量綱，其單位通常采用bin標(biāo)示；=[]?N，符號[ ]表示四舍五入取整，由D=s/可知，非整余數(shù)?[-0.5, 0.5]，則點離散信號()的廣義DTFT可表示為

（2）

式中，j為復(fù)數(shù)單位；(?)為寬帶噪聲的DTFT；第二項為負(fù)頻率分量引起的頻譜長泄漏干擾。值得注意的是，此處參量可為任意[0,-1]范圍內(nèi)的實數(shù)，當(dāng)其為整數(shù)時，式（2）即為標(biāo)準(zhǔn)DFT，此時()的等間隔抽樣即可表示為(|?N)。

便于推導(dǎo)，首先忽略噪聲項(?)，通過計算可將式（2）轉(zhuǎn)化為

將式（2）改寫為式（3）所示方程具有兩個非常重要的作用：一方面，式（3）中含有3個未知參量、和，可知還需另外兩條DTFT譜線{(±)}以形成方程組進行解析；更為重要的一方面在于式（3）左右側(cè)為復(fù)數(shù)，自然地可分解為對應(yīng)實部和虛部的方程，這是求解方程組并導(dǎo)出解析插值頻率估計的關(guān)鍵。

1.2 齊次線性方程組

如式（3）左側(cè)所示，()e-jpk(1-N)/N與未知參量無關(guān)，僅由()及其位置索引決定，簡單起見采用 =()e-jpk(1-N)/N表示，則據(jù)式（3）可得，對應(yīng)實部的齊次線性方程組，以矩陣形式表示為

方程組式（4）中系數(shù)矩陣定義如下

其中

同理，對應(yīng)虛部的齊次線性方程組為

方程組式（9）中系數(shù)矩陣的相關(guān)向量定義如下

在式（4）～式（12）中，需要明確指出的是?R+，其范圍控制在(0, 0.5]，使得所用DTFT譜線全部位于矩形窗譜的主瓣范圍內(nèi)，以確保具有高信噪比。

1.3 頻率估計

由于齊次線性方程組式（4）和式（9）必有非零解，則根據(jù)線性代數(shù)中克萊姆法則（Cramer’s Rule）可知方程組系數(shù)矩陣的行列式必然為零[22]，即

雖然據(jù)1.2節(jié)中各矩陣定義可知，上述方程都僅含有單一的未知變量，至此可以單獨求解式（13）或式（14）而直接得到的解析解。但是，當(dāng)所用DTFT譜線的實部或虛部接近零時，噪聲干擾將導(dǎo)致估計極為不準(zhǔn)確。為防止這種情況的出現(xiàn)，需結(jié)合式（13）和式（14）得到用于頻率估計的解析插值公式（具體求解過程見附錄第1節(jié)），結(jié)果如下

其中

2 頻率估計的統(tǒng)計特性分析

第1節(jié)通過研究無噪聲情況，得到了基于廣義插值DTFT的頻率估計解析式，本節(jié)將分析加性高斯白噪聲背景下頻率估計算法的性能并得到估計的漸近性質(zhì)，進而通過簡單迭代優(yōu)化算法性能。

據(jù)式（23）可知，頻率估計由決定，因而首先分析噪聲對估計的影響。令=( )e-jpk(1-N)/N，則背景噪聲干擾下的估計應(yīng)準(zhǔn)確表示為

其中

由此可得噪聲干擾下的估計偏差D為

其中

由此利用泰勒級數(shù)展開，式（23）所示的頻率估計解析式可改寫為

則噪聲干擾下頻率估計的統(tǒng)計偏差為

由白噪聲特性(1)=(2)=0可知，期望值(D)= 0，表明式（23）為頻率的無偏估計。

進一步考慮噪聲項之間的相關(guān)性，可得頻率估計的方差（詳細推導(dǎo)過程見附錄第2節(jié)）為

式中，[Re2(D )]如式（A20）所示。

進一步分析可知，任意頻點下式（31）的變化趨勢及其量級大小始終由|1-G2-G3|決定，因而通過分析|1-G2-G3|即可掌握頻率估計的統(tǒng)計方差變化。根據(jù)中對應(yīng)各定義，將式（2）和式（21）代入1-G2-G3即可簡化為

式中，=(/2)exp{j[+p(-1)/]}；′=2+；()定義為

通常情況下＞1，此時中的負(fù)頻率泄漏分量可忽略不計，在中心譜線位置搜索正確的前提下（即=），1-G2-G3可近似為

根據(jù)式（34）可知：

（1）＞0。當(dāng)Re()=0和Im()=1時，對應(yīng)|1-G2-G3|的最大值，此時單一頻點下式（31）僅包含Im(1)時接近其最小值（下界）；當(dāng)Re()=1和Im()=0時，對應(yīng)|1-G2-G3|的最小值，此時單一頻點下式（31）僅包含Re(1)時接近其最大值（上界）。

（2）≤0。當(dāng)Re()=1和Im()=0時，對應(yīng)|1-G2-G3|的最大值，此時單一頻點下式（31）僅包含Im(1)時接近其最小值（下界）；當(dāng)Re()=0和Im()=1時，對應(yīng)|1-G2-G3|的最小值，此時單一頻點下式（31）僅包含Re(1)時接近其最大值（上界）。

實際中由于信號相位為隨機分布，因此取其中值作為平均估計，即Re()=Im()=cos(p/2)。以= 1 024，=2，s=1 024Hz，SNR=40dB，?[2.55, 3.45]和=0.1為例，對應(yīng)每一個時頻率估計方差的理論分布如圖1所示。

圖1 i=0.1時頻率估計方差的范圍

圖1結(jié)果顯示，在=0.1且相位隨機分布前提下，所提頻率估計方法的方差在→0時接近最小值，且最小值十分接近無偏頻率估計的CRLB[11]。進一步地，考慮相位在[-p,p]范圍內(nèi)隨機變化，值在0.1～0.5范圍內(nèi)逐漸增大，頻率估計方差理論均值的最小值及其與CRLB之比見表1。

表1 頻率估計的方差最小值

Tab.1 The minimum variance of the frequency estimates

據(jù)表1結(jié)果可知，隨著逐漸增大，頻率估計方差最小值逐漸偏離CRLB，但差距十分微小。這表明，算法在≤0.5范圍內(nèi)性能十分穩(wěn)定且受值的影響極小，理論上，隨著→0頻率估計方差將完全逼近無偏頻率估計的CRLB。

另一更為重要且關(guān)鍵的點在于：該結(jié)果表明算法具有可迭代執(zhí)行的能力，即通過頻移（不斷迭代），將≠0時的頻率估計轉(zhuǎn)化為→0時的頻率估計，進一步優(yōu)化提升算法性能，使得≠0時的頻率估計方差能一致逼近于CRLB，從而最大程度地降低噪聲對算法的干擾影響。迭代執(zhí)行程序的具體執(zhí)行步驟如下所示：

（2）通過峰值搜索(κ)=argmax{|DFT()|}得到中心譜線，根據(jù)DTFT定義計算輔助譜線(κ±)。

（4）=+1。

（6）據(jù)DTFT定義計算(κ+1)及輔助譜線(κ+1)。

（9）否則，重復(fù)步驟（4）～步驟（8）。

值得注意的是，雖然理論上可取任意接近于零的數(shù)值，但當(dāng)極小時有限字長效應(yīng)可能導(dǎo)致式（23）的分子和分母接近于零，使得頻率估計可能產(chǎn)生無意義的結(jié)果[23]；且當(dāng)?shù)陀谀骋婚撝禃r，受限于采樣點數(shù)和計算精度限制，頻率估計方差幾乎保持恒定而無法隨著減小而降低，具體結(jié)果見3.1節(jié)。

3 仿真結(jié)果

為全面驗證算法的有效性和準(zhǔn)確性，分別進行三個不同的仿真測試。第一個為分析所用譜線的間隔對算法在非迭代情況下性能的影響；第二個為驗證算法在迭代情況下是否能夠一致逼近無偏頻率估計的CRLB；第三個為不同性能優(yōu)異的插值算法與本文算法的性能比對。

3.1 譜線間隔i的影響分析

噪聲干擾情況下，據(jù)附錄第2節(jié)可知，單一譜線及與其他譜線間的統(tǒng)計特性隨著譜線間隔的變化而變化，進而導(dǎo)致不同值下頻率估計方差存在差異。設(shè)定信號采樣點數(shù)=1 024，幅度=2，采樣頻率s=1 024Hz，信號初始相位在[-p,p]范圍內(nèi)隨機變化，對應(yīng)信號波動周期在[2.55, 3.45]以步長0.1逐漸變化；噪聲干擾設(shè)置為加性高斯白噪聲，對應(yīng)信噪比為SNR=40dB，單一信號條件下重復(fù)= 10 000次獨立蒙特卡洛仿真測試，則頻率估計的MSE與CRLB之比隨的變化趨勢如圖2所示。

圖2 頻率估計MSE

圖2結(jié)果顯示：①頻率估計的MSE與理論方差結(jié)果一致吻合，表明噪聲擾動下算法為嚴(yán)格無偏頻率估計。②當(dāng)≤0.5時，頻率估計MSE在單一頻率條件下幾乎保持恒定，且在→0時十分接近于CRLB；反之，則隨著的增大而逐漸增大。其原因在于當(dāng)≤0.5時，由于所用三條譜線都位于矩形窗主瓣范圍內(nèi)，因而其信噪比相對較高；反之，隨著旁瓣衰減，輔助譜線的信噪比則很低，受噪聲擾動影響較大。③當(dāng)≤0.5時，頻率估計MSE不受頻率偏移符號影響，即在＞0或＜0范圍內(nèi)的MSE在量級上幾乎相等；當(dāng)＞0.5時（特別是在→1時），在＞0范圍內(nèi)的頻率估計MSE要小于＜0范圍內(nèi)，其原因在于當(dāng)＜0時，量級較小的輔助譜線受負(fù)頻率的頻譜泄漏影響將不可忽略。

3.2 迭代有效性驗證

仍以3.1節(jié)中設(shè)定的信號為例，對應(yīng)信號波動周期在[2.55, 3.45]則以步長0.02逐漸變化。根據(jù)3.1節(jié)仿真結(jié)果，此處僅考慮≤0.5范圍內(nèi)取值時的頻率估計結(jié)果，即對比＞0.5時保持相對最優(yōu)；同時，考慮實際應(yīng)用中數(shù)值穩(wěn)定性，此處僅示出=0.1時的頻率估計結(jié)果，如圖3所示。

圖3 i=0.1時頻率估計的MSE

圖3結(jié)果顯示，噪聲干擾下雖然所提無偏頻率估計算法在頻偏偏離0時無法一致逼近CRLB，但伴隨迭代程序的執(zhí)行，僅需1次迭代頻率估計的MSE在頻偏?[-0.5, 0.5]全域內(nèi)都一致趨近于最優(yōu)估計結(jié)果（→0時），即在任意頻偏情況下都十分接近于CRLB。在不顯著增加計算復(fù)雜度的前提下，迭代執(zhí)行程序有效降低了算法對頻率偏移的敏感度，極大地提高了算法的魯棒性。對于≤0.5范圍內(nèi)其他任意值，據(jù)圖2和第2節(jié)中算法統(tǒng)計特性分析可知，迭代程序仍有效且在噪聲干擾下頻偏估計結(jié)果將十分接近，不再贅述。

3.3 性能比對

為進一步驗證算法在不同波動周期下的有效性和準(zhǔn)確性，設(shè)定信號波動周期在[0.55, 4.45]范圍內(nèi)逐漸變化，其他信號參數(shù)保持與3.1節(jié)和3.2節(jié)一致。由于頻域類頻率估計算法計算復(fù)雜度通常要遠遠小于時域類參數(shù)化算法，分別為(log2)和(3)，因此，僅選取目前為數(shù)不多的幾種考慮負(fù)頻率泄漏影響的頻域插值算法進行性能比對，包括基于頻譜系數(shù)加權(quán)的3IpDFT[18]、采用矩陣解析的3IpDFT[15]、基于補零和DFT系數(shù)變換的2IpDFT[11]、解析插值算法3IpDFT[10]以及改進的迭代插值DTFT算法AM[20]。不同插值算法的頻率估計MSE結(jié)果如圖4所示。圖4a為SNR=40dB時?[0.55, 4.45]范圍內(nèi)的頻率估計結(jié)果；圖4b為本文算法單步迭代和AM多步迭代收斂的頻率估計結(jié)果比對；圖4c為SNR?[-20, 100]dB范圍內(nèi)=1.25時的頻率估計結(jié)果。此外，圖4中以IpDTFT(,)表示本文算法，其中為譜線間隔，為迭代次數(shù)；鑒于在≤0.5范圍內(nèi)算法性能十分接近，此處仍僅示出=0.1時的仿真結(jié)果。

圖4 不同算法頻率估計性能比對

圖4a結(jié)果顯示：①非迭代情況下，本文算法和2IpDFT整體上要遠優(yōu)于其他參考比對算法。②當(dāng)信號波動周期≤1時，本文算法要優(yōu)于2IpDFT；當(dāng)信號波動周期＞1時，在偏移||≤0.25范圍內(nèi)兩種算法性能十分接近；在0.25＜||≤0.5范圍內(nèi)2IpDF要優(yōu)于本文算法。其原因在于，當(dāng)||≤0.25時兩種算法的解析方式極為接近，進一步分析可知，當(dāng)且僅當(dāng)=0.5時兩者具有幾乎相同的解析插值公式，由3.1節(jié)結(jié)果可知，兩種算法性能相當(dāng)；而當(dāng) 0.25＜||≤0.5時，由于2IpDFT算法通過補零改變了解析插值中所用譜線，使得其在該范圍內(nèi)的性能與||≤0.25時保持一致，因而優(yōu)于本文算法的非迭代結(jié)果。

圖4b結(jié)果顯示：①當(dāng)信號波動周期＞1時，通過多次不斷迭代收斂的AM算法與本文算法單次迭代結(jié)果幾乎保持一致，明顯受限于多次迭代計算，AM算法計算量將遠遠高于本文算法；②當(dāng)信號波動周期≤1時，受限于與輔助譜線間隔固定為=0.5，AM算法中所用譜線量級極小，受負(fù)頻率的頻譜泄漏影響遠大于本文算法（=0.1），且因未采用信噪比較高的中心峰值譜線，導(dǎo)致算法MSE要遠遠高于本文算法。

圖4c結(jié)果顯示：①在信號波動周期=1.25時，此時負(fù)頻率頻譜泄漏干擾相對占主導(dǎo)，非迭代AM算法和3IpDFT算法存在明顯的估計偏差，因而性能上無法一致逼近CRLB；同理，迭代AM算法在信噪比SNR＞80dB時，因此處迭代次數(shù)設(shè)定限制無法完全消除估計偏差影響，導(dǎo)致其MSE偏離CRLB。②總體上本文算法（包括迭代）和其他算法對比，因準(zhǔn)確計及負(fù)頻率泄漏使得頻率估計結(jié)果幾乎一致逼近CRLB，在不同SNR情況下本文算法的單次迭代結(jié)果依然保持最優(yōu)。

綜上可知，相較其他插值類算法，雖然本文算法在非迭代情況下未保持全域最優(yōu)，但在不明顯增加計算量的前提下，僅通過單步迭代以優(yōu)化性能，即能達到最優(yōu)的頻率估計且全域一致逼近CRLB。

表2 不同算法計算復(fù)雜度比對

Tab.2 Comparison of the computational complexity of different algorithms

4 實驗結(jié)果

實驗驗證測試平臺如圖5所示，以精度為0.05級的DK-56B1三相交直流指示儀表檢定裝置作為功率信號源，額定輸出信號頻率為50Hz，信號幅值為5V；以24bit的IOtech數(shù)據(jù)記錄儀作為信號采集裝置，采樣頻率為100kHz。以截取信號采樣點數(shù)的改變對應(yīng)于變化，對應(yīng)于每個通過單點移動改變信號初始相位，并連續(xù)進行10 000次分析。此外，鑒于白噪聲擾動下IEEE 1057標(biāo)準(zhǔn)[24]建議的4參量正弦擬合算法（4PSF）提供了頻率的最大似然估計，其估計方差非常接近CRLB，因此以該算法估計結(jié)果作為基準(zhǔn)，其迭代次數(shù)設(shè)定為6。在?[1, 4]范圍內(nèi)不同算法的頻率估計MSE實驗結(jié)果如圖6所示。

圖5 實驗平臺

圖6 實驗結(jié)果

從圖6可看出，實驗結(jié)果與圖4仿真結(jié)果基本完全吻合。此處，雖然實驗中頻率分辨率隨著增大而增高，頻率估計的方差下界將逐漸降低，但本文算法的單步迭代結(jié)果在任意頻點上仍然與4PSF基準(zhǔn)結(jié)果基本保持一致，這也再次印證了本文算法的準(zhǔn)確性和迭代有效性。

圖7 實測電力信號

圖8 實測電壓信號的頻率估計

如圖8a所示，當(dāng)采樣點數(shù)較小即≤175時，實測電壓信號含有的一定量諧波泄漏干擾導(dǎo)致頻率估計波動偏大；隨著采樣點數(shù)增多，諧波泄漏干擾逐漸降低直至忽略不計，頻率估計趨于穩(wěn)定且能夠準(zhǔn)確反映負(fù)載運行狀態(tài)變化。圖8b所示頻率估計MSE隨變化而變化的趨勢與圖6實驗室測試結(jié)果基本一致，雖然在較小時受限于諧波泄漏干擾，導(dǎo)致與實驗室測試結(jié)果出現(xiàn)一定偏差，但總體上本文算法的單步迭代頻率估計結(jié)果依然保持最優(yōu)，而且極為接近于以4PSF算法為基準(zhǔn)近似的CRLB。

5 結(jié)論

針對短時波動周期下實正弦信號頻率精確估計問題，首先利用頻譜實/虛部線性方程組解析得到了基于廣義DTFT插值的無偏頻率估計公式，不僅擺脫了對窗函數(shù)旁瓣衰減特性的嚴(yán)重依賴，而且打破了以往傳統(tǒng)頻域插值算法僅限于單頻解析信號和離散DFT譜線插值的局限性；進而基于噪聲背景下頻率估計的統(tǒng)計性能分析導(dǎo)出了算法的迭代執(zhí)行方式，解決了算法在背景噪聲干擾下無法保持全域最優(yōu)且一致逼近CRLB的問題。

仿真和實驗結(jié)果顯示，一方面，得益于算法解析過程準(zhǔn)確計及來自負(fù)頻率的頻譜泄漏干擾，且采用間隔低至0.1D的DTFT輔助譜線，即便信號波動周期≤1時算法在未迭代情況下性能仍保持穩(wěn)定；另一方面，對于噪聲干擾影響，算法通過簡單有效的迭代執(zhí)行方式以優(yōu)化估計性能，在不顯著增加計算量和復(fù)雜度的前提下，僅需單步迭代即能達到最優(yōu)的頻率估計且全域一致逼近CRLB，在=0.1時理論上算法最小MSE與CRLB之比接近1.008 822。

對于諧波擾動的多頻信號頻率估計，在主瓣頻譜不存在混疊的前提下，可首先采用本文算法得到單頻點的頻率粗估計，進一步通過單頻點的粗估計實現(xiàn)頻譜重構(gòu)，逐步迭代消除諧波頻譜長泄漏的干擾影響，從而再次采用本文算法實現(xiàn)精確頻率估計。

這種針對多頻信號頻率估計的算法具體實現(xiàn)及相應(yīng)測試結(jié)果將在下一步研究工作中詳細探討。

附 錄

1.的具體求解過程

便于推導(dǎo)，首先將矩陣+/-和+/-中與位置索引相關(guān)的已知三角函數(shù)變量分別以a和b表示為

由此，式（13）和式（14）改寫為

則根據(jù)行列式性質(zhì)將式（A2）左側(cè)行列式分別進行簡化，即

再次利用行列式性質(zhì)，將式（A6）虛數(shù)部分乘以虛數(shù)單位j后與式（A5）實數(shù)部分進行合并，進而經(jīng)過展開計算后可得的解析解為

2. 頻率估計的統(tǒng)計方差推導(dǎo)

考慮中噪聲項之間的相關(guān)性，首先將D展 開，即

令=-，=-和=-，由于[Re(D)]=0，Re(D)的自相關(guān)與方差相等，即

式中，由于各參量方差及互相關(guān)具有相似性，因此，此處僅分別以的方差以及和的互相關(guān)為例進行推導(dǎo)，即的方差為

其中

（A13）

其中

由此，據(jù)式（A10）～式（A18）可得Re(D)的方差為

式中，P和N分別為

[1] Sun J, Aboutanios E, Smith D B, et al. Robust frequency, phase, and amplitude estimation in power Systems considering harmonics[J]. IEEE Transactions on Power Delivery, 2020, 35(3): 1158-1168.

[2] 張政, 溫和, 黎福海, 等. 多水平集單周期電力系統(tǒng)頻率測量方法及應(yīng)用[J]. 電工技術(shù)學(xué)報, 2017, 32(7): 119-127.

Zhang Zheng, Wen He, Li Fuhai, et al. Multi-level set single-cycle power system frequency measurement method and its application[J]. Transactions of China Electrotechnical Society, 2017, 32(7): 119-127.

[3] 肖勇, 趙偉, 黃松嶺. 基于離散傅里葉級數(shù)的非同步采樣下諧波功率測量算法[J]. 電工技術(shù)學(xué)報, 2018, 33(7): 1570-1578.

Xiao Yong, Zhao Wei, Huang Songling. Harmonic power measurement algorithm based on discrete fourier series in asynchronous sampling[J]. Transa- ctions of China Electrotechnical Society, 2018, 33(7): 1570-1578.

[4] Jafarpisheh B, Madani S M, Parvaresh F, et al. Power system frequency estimation using adaptive accelerated MUSIC[J]. IEEE Transactions on Instrumentation and Measurement, 2018, 67(11): 2592-2602.

[5] 楊超, 張淮清, 王耀, 等. 計及全泄漏影響的多點插值離散傅里葉變換校正方法[J]. 電工技術(shù)學(xué)報, 2020, 35(16): 3385-3395.

Yang Chao, Zhang Huaiqing, Wang Yao, et al. Multipoint interpolated discrete fourier transform correction method considering total leakage effect[J]. Transactions of China Electrotechnical Society, 2020, 35(16): 3385-3395.

[6] Shin D, Kwak C, Kim G. An efficient algorithm for frequency estimation from cosine-sum windowed DFT coefficients[J]. Signal Processing, 2020, 166: 107245.1-107245.9.

[7] 孫仲民, 何正友, 臧天磊. 一種混合卷積窗及其在諧波分析中的應(yīng)用[J]. 電工技術(shù)學(xué)報, 2016, 31(16): 207-214.

Sun Zhongmin, He Zhengyou, Zang Tianlei. A kind of hybrid convolution window and its application in harmonic analysis[J]. Transactions of China Electro- technical Society, 2016, 31(16): 207-214.

[8] Romano P, Paolone M. Enhanced interpolated-DFT for synchrophasor estimation in FPGAs: theory, implementation, and validation of a PMU prototype[J]. IEEE Transactions on Instrumentation and Measure- ment, 2014, 63(12): 2824-2836.

[9] 翟曉軍, 周波. 一種改進的插值FFT諧波分析算法[J]. 中國電機工程學(xué)報, 2016, 36(11): 2952-2958.

Zhai Xiaojun, Zhou Bo. An improved interpolated FFT algorithm for harmonic analysis[J]. Proceedings of the CSEE, 2016, 36(11): 2952-2958.

[10] Wang Kai, Wen He, Li Guoqing. Accurate frequency estimation by using three-point interpolated discrete fourier transform based on rectangular window[J]. IEEE Transactions on Industrial Informatics, 2021, 17(1): 73-81.

[11] Sun Yuxin, Zhuang Chungang, Xiong Zhenhua. A scale factor based interpolated DFT for chatter frequency estimation[J]. IEEE Transactions on Instrumentation and Measurement, 2015, 64(10): 2666-2678.

[12] Sun Yuxin, Zhuang Chungang, Xiong Zhenhua. A switch-based inter-polated DFT for the small number of acquired sine wave cycles[J]. IEEE Transactions on Instrumentation and Measurement, 2016, 65(4): 846-855.

[13] Hlupic N, Basch D, Fertalj K. A note on optimality of analytical leakage compensation at boundary fre- quencies[J]. IEEE Transactions on Instrumentation and Measurement, 2010, 59(2): 301-308.

[14] Borkowski J, Kania D, Mroczka J. Interpolated-DFT- based fast and accurate frequency estimation for the control of power[J]. IEEE Transactions on Industrial Electronics, 2014, 61(12): 7026-7034.

[15] Chen Kuifu, Cao Xu, Li Yangfeng. Sine wave fitting to short records initialized with the frequency retrieved from Hanning windowed FFT spectrum[J]. Measurement, 2009, 42(1): 127-135.

[16] Murakami T, Wang W. An analytical solution to Jacobsen estimator for windowed signals[C]//IEEE International Conference on Acoustics Speech and Signal Processing (ICASSP), Barcelona, Spain, 2020: 5950-5954.

[17] Belega D, Petri D. Effect of noise and harmonics on sine-wave frequency estimation by interpolated DFT algorithms based on few observed cycles[J]. Signal Processing, 2017, 140: 207-218.

[18] Belega D, Petri D, Dallet D. Frequency estimation of a sinusoidal signal via a three-point interpolated DFT method with high image component interference rejection capability[J]. Digital Signal Processing, 2014, 24: 162-169.

[19] Candan C. A method for fine resolution frequency estimation from three DFT samples[J]. IEEE Signal Process, Lett, 2011, 18(6): 351-354.

[20] Ye S, Sun J, Aboutanios E. On the estimation of the parameters of a real sinusoid in noise[J]. IEEE Signal Processing Letters, 2017, 24(5): 638-642.

[21] Aboutanios E, Mulgrew B. Iterative frequency estimation by interpolation on Fourier coefficients[J]. IEEE Transactions on Signal Processing, 2005, 53(4): 1237-1242.

[22] Horn R A, Johnson C R. Matrix analysis[M]. Cambridge: Cambridge University Pre, 1985.

[23] Lei Fan, Guoqing Qi. Frequency estimator of sinusoid based on interpolation of three DFT spectral lines[J]. Signal Processing, 2018, 144: 52-60.

[24] IEEE standard for digitizing waveform recorders[S]. New Jersey: Piscataway, 2007.

Accurate Frequency Estimation Based on Analytical Interpolated Discrete Time Fourier Transform

11,21,31

（1. State Key Laboratory of Power Transmission Equipment & System Security and New Technology Chongqing University Chongqing 400044 China 2. Electric Power Research Institute of Yunnan Power Grid Co. Ltd Kunming 650217 China 3. Chongqing Triloop Prospecting Technology Co. Ltd Chongqing 402660 China）

When the cycles in record (CiR) of the signal is small, the spectral leakage from the image component is the main source of estimation errors in the frequency-domain interpolation. To eliminate the influence of spectral leakage and accurately estimate the frequency of a real sinusoidal signal with any CiR, an analytic interpolation algorithm based on generalized DTFT is proposed. The algorithm uses any DTFT spectral lines located in the main lobe of the spectrum to construct a linear equation set corresponding to the real and imaginary parts, and then an unbiased frequency estimator based on analytic interpolated DTFT is obtained. Besides, considering the influence of additive white noise in practice, a simple and effective iterative method is presented based on the analysis of the statistical characteristics of the frequency estimator. The iterative method ensures that the mean square errors (MSE) of the frequency estimation can consistently achieve the Cramer-Rao lower bound (CRLB) of the unbiased frequency estimation under any CiR. Simulation and experimental results confirm that the proposed algorithm is superior to other existing interpolated DFT algorithms.

Small cycles in record, spectral leakage, analytically interpolated discrete time fourier transform (DTFT), iterative method, Cramer-Rao lower bound (CRLB)

10.19595/j.cnki.1000-6753.tces.201667

TM935

楊 超 男，1986年生，博士，研究方向為數(shù)字信號處理及其在精密測量、電能計量與電能質(zhì)量分析中的深度應(yīng)用。E-mail: 328383369@qq.com（通信作者）

李 波 男，1982年生，碩士，高級工程師，主要從事電氣測量以及智能電網(wǎng)測控研究工作。E-mail: libo2010@163.com

2020-12-20

2021-04-27

國家重點研發(fā)計劃課題（2018YFC0406904）和重慶市重點產(chǎn)業(yè)共性關(guān)鍵技術(shù)創(chuàng)新專項（重點研發(fā)項目）（CSTC2017ZDCY-ZDYFX0045）資助。

（編輯 陳 誠）