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

楊 超 李 波,2 胡緒權,3 付志紅

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

楊 超1李 波1,2胡緒權1,3付志紅1

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

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

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

0 引言

正弦信號的頻率估計是一個經典且重要的信號處理問題，在科學和工程領域具有極為廣泛的應用。在電力系統中，頻率作為最重要的電能質量參數之一，其變化反映了發電和負載之間的動態平衡[1-3]，通過頻率及其變化的準確估計可以實現系統的動態過程監測。此外，在計量基準測試和標準溯源等特定應用中，對頻率估計精度往往具有更高的要求。

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

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

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

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

1 頻率估計的解析解

1.1 基本原理

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

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

（2）

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

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

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

1.2 齊次線性方程組

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

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

其中

同理，對應虛部的齊次線性方程組為

方程組式（9）中系數矩陣的相關向量定義如下

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

1.3 頻率估計

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

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

其中

2 頻率估計的統計特性分析

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

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

其中

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

其中

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

則噪聲干擾下頻率估計的統計偏差為

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

進一步考慮噪聲項之間的相關性，可得頻率估計的方差（詳細推導過程見附錄第2節）為

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

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

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

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

根據式（34）可知：

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

（2）≤0。當Re()=1和Im()=0時，對應|1-G2-G3|的最大值，此時單一頻點下式（31）僅包含Im(1)時接近其最小值（下界）；當Re()=0和Im()=1時，對應|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為例，對應每一個時頻率估計方差的理論分布如圖1所示。

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

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

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

Tab.1 The minimum variance of the frequency estimates

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

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

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

（4）=+1。

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

（9）否則，重復步驟（4）～步驟（8）。

值得注意的是，雖然理論上可取任意接近于零的數值，但當極小時有限字長效應可能導致式（23）的分子和分母接近于零，使得頻率估計可能產生無意義的結果[23]；且當低于某一閾值時，受限于采樣點數和計算精度限制，頻率估計方差幾乎保持恒定而無法隨著減小而降低，具體結果見3.1節。

3 仿真結果

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

3.1 譜線間隔i的影響分析

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

圖2 頻率估計MSE

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

3.2 迭代有效性驗證

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

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

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

3.3 性能比對

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

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

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

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

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

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

表2 不同算法計算復雜度比對

Tab.2 Comparison of the computational complexity of different algorithms

4 實驗結果

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

圖5 實驗平臺

圖6 實驗結果

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

圖7 實測電力信號

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

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

5 結論

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

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

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

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

附 錄

1.的具體求解過程

便于推導，首先將矩陣+/-和+/-中與位置索引相關的已知三角函數變量分別以a和b表示為

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

則根據行列式性質將式（A2）左側行列式分別進行簡化，即

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

2. 頻率估計的統計方差推導

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

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

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

其中

（A13）

其中

由此，據式（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] 張政, 溫和, 黎福海, 等. 多水平集單周期電力系統頻率測量方法及應用[J]. 電工技術學報, 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] 肖勇, 趙偉, 黃松嶺. 基于離散傅里葉級數的非同步采樣下諧波功率測量算法[J]. 電工技術學報, 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]. 電工技術學報, 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] 孫仲民, 何正友, 臧天磊. 一種混合卷積窗及其在諧波分析中的應用[J]. 電工技術學報, 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]. 中國電機工程學報, 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年生，博士，研究方向為數字信號處理及其在精密測量、電能計量與電能質量分析中的深度應用。E-mail: 328383369@qq.com（通信作者）

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

2020-12-20

2021-04-27

國家重點研發計劃課題（2018YFC0406904）和重慶市重點產業共性關鍵技術創新專項（重點研發項目）（CSTC2017ZDCY-ZDYFX0045）資助。

（編輯 陳 誠）