基于解析插值離散時(shí)間傅里葉變換的精確頻率估計(jì)

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

基于解析插值離散時(shí)間傅里葉變換的精確頻率估計(jì)

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

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

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

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

0 引言

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

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

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

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

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

1 頻率估計(jì)的解析解

1.1 基本原理

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

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

（2）

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

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

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

1.2 齊次線性方程組

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

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

其中

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

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

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

1.3 頻率估計(jì)

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

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

其中

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

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

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

其中

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

其中

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

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

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

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

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

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

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

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

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

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

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

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

圖1 i=0.1時(shí)頻率估計(jì)方差的范圍

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

表1 頻率估計(jì)的方差最小值

Tab.1 The minimum variance of the frequency estimates

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

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

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

（4）=+1。

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

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

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

3 仿真結(jié)果

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

3.1 譜線間隔i的影響分析

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

圖2 頻率估計(jì)MSE

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

3.2 迭代有效性驗(yàn)證

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

圖3 i=0.1時(shí)頻率估計(jì)的MSE

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

3.3 性能比對(duì)

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

圖4 不同算法頻率估計(jì)性能比對(duì)

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

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

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

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

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

Tab.2 Comparison of the computational complexity of different algorithms

4 實(shí)驗(yàn)結(jié)果

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

圖5 實(shí)驗(yàn)平臺(tái)

圖6 實(shí)驗(yàn)結(jié)果

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

圖7 實(shí)測(cè)電力信號(hào)

圖8 實(shí)測(cè)電壓信號(hào)的頻率估計(jì)

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

5 結(jié)論

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

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

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

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

附 錄

1.的具體求解過程

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

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

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

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

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

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

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

式中，由于各參量方差及互相關(guān)具有相似性，因此，此處僅分別以的方差以及和的互相關(guān)為例進(jì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)頻率測(cè)量方法及應(yīng)用[J]. 電工技術(shù)學(xué)報(bào), 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í)數(shù)的非同步采樣下諧波功率測(cè)量算法[J]. 電工技術(shù)學(xué)報(bào), 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ì)及全泄漏影響的多點(diǎn)插值離散傅里葉變換校正方法[J]. 電工技術(shù)學(xué)報(bào), 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é)報(bào), 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] 翟曉軍, 周波. 一種改進(jìn)的插值FFT諧波分析算法[J]. 中國(guó)電機(jī)工程學(xué)報(bào), 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年生，博士，研究方向?yàn)閿?shù)字信號(hào)處理及其在精密測(cè)量、電能計(jì)量與電能質(zhì)量分析中的深度應(yīng)用。E-mail: 328383369@qq.com（通信作者）

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

2020-12-20

2021-04-27

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

（編輯 陳 誠(chéng)）