本發(fā)明涉及信號譜估計領(lǐng)域，更具體地說，涉及一種基于自適應(yīng)濾波器幅頻響應(yīng)的高精度頻率估計方法。

背景技術(shù)：

信號譜估計是指對一個有限長的平穩(wěn)序列，估計其在整個頻率域內(nèi)的功率分布。信號譜估計在通信、雷達(dá)、聲吶、地震學(xué)等諸多領(lǐng)域中均有廣泛應(yīng)用。以雷達(dá)和聲吶系統(tǒng)為例，人們可以根據(jù)接收信號的頻譜對觀測范圍內(nèi)的目標(biāo)進(jìn)行定位。根據(jù)功率譜的分布情況，信號譜估計可以分為兩類即：連續(xù)譜估計和線譜估計，其中線譜估計的主要任務(wù)就是對包含正弦信號和噪聲的觀測信號進(jìn)行處理，以獲得正弦信號的頻率估計值。

已有的頻率估計方法包括：快速傅里葉變換(Fast Fourier Transformation,FFT)方法、多重信號分類(Multiple Signal Classification,MUSIC)方法、求根多重信號分類(Root-MUSIC)方法、信號參數(shù)估計旋轉(zhuǎn)不變技術(shù)(Estimating Signal Parameters via Rotational Invariance Techniques,ESPRIT)等。其中，F(xiàn)FT方法計算量小且易于實現(xiàn)，因此在實際工程中得到了廣泛的應(yīng)用，但是FFT方法的頻率分辨率不夠高；MUSIC方法、Root-MUSIC方法以及ESPRIT方法屬于子空間處理方法，它們在高信噪比情況下可獲得更高的頻率分辨率。但是，當(dāng)信噪比很低且信號采樣點數(shù)很少時，不管是常規(guī)的FFT方法，還是MUSIC、Root-MUSIC及ESPRIT等子空間方法，其估計性能均顯著惡化，所得到的頻率估計精度遠(yuǎn)遠(yuǎn)不能達(dá)到實際應(yīng)用的要求。因此，研究一種適用于信噪比很低、信號采樣點數(shù)很少等情況下的高精度頻率估計方法具有重要的實用價值。

技術(shù)實現(xiàn)要素：

為了有效解決已有頻率估計方法在信噪比很低且信號采樣點數(shù)很少等情況下頻率估計性能顯著惡化的問題，本發(fā)明提供了一種基于自適應(yīng)濾波器幅頻響應(yīng)的高精度頻率估計方法。在信噪比很低和信號采樣點數(shù)很少等情況下，本發(fā)明方法的頻率估計精度明顯優(yōu)于已有頻率估計方法。

本發(fā)明的基本思路是：首先，對觀測的連續(xù)時間信號進(jìn)行采樣得到離散時間信號，并由離散時間信號計算得到信號協(xié)方差矩陣；其次，由先驗信息確定第一個正弦信號的頻率區(qū)間；再次，由信號協(xié)方差矩陣和頻率區(qū)間計算自適應(yīng)濾波器的權(quán)系數(shù)向量，并使自適應(yīng)濾波器的幅頻響應(yīng)滿足以下特點：當(dāng)自適應(yīng)濾波器的中心頻率不等于第一個正弦信號的頻率時，自適應(yīng)濾波器的幅頻響應(yīng)在頻率區(qū)間內(nèi)存在明顯的零陷，而當(dāng)自適應(yīng)濾波器的中心頻率等于第一個正弦信號的頻率時，自適應(yīng)濾波器的幅頻響應(yīng)在頻率區(qū)間內(nèi)不存在明顯的零陷；然后，由自適應(yīng)濾波器的權(quán)系數(shù)向量計算自適應(yīng)濾波器在頻率區(qū)間內(nèi)的幅頻響應(yīng)；接著，在頻率區(qū)間內(nèi)搜索自適應(yīng)濾波器的幅頻響應(yīng)的最小值及最大值，并由上述最小值和最大值計算信號偽功率譜；最后，在頻率區(qū)間內(nèi)搜索信號偽功率譜的譜峰，并將信號偽功率譜的譜峰對應(yīng)的頻率作為第一個正弦信號的頻率估計值。若在觀測信號中正弦信號的個數(shù)大于1，則根據(jù)先驗信息確定其它正弦信號的頻率區(qū)間，并重復(fù)上述步驟即可得到其它正弦信號的頻率估計值。

本發(fā)明的技術(shù)方案是：一種基于自適應(yīng)濾波器幅頻響應(yīng)的高精度頻率估計方法，具體包括下述步驟：

第一步：計算信號協(xié)方差矩陣

首先，對觀測的連續(xù)時間信號進(jìn)行采樣得到離散時間信號{x(1),x(2),...,x(N)}，其中采樣頻率為f_s，單位為Hz，N表示離散時間信號的采樣點數(shù)。然后，由離散時間信號{x(1),x(2),...,x(N)}構(gòu)造信號向量x(k)，即：

x(k)＝[x(k-m+1),x(k-m+2),...,x(k)]^H

其中，m為信號向量x(k)的長度，m為正整數(shù)，且1＜m＜N，k也為正整數(shù)，且m≤k≤N，m的具體大小根據(jù)實際應(yīng)用情況確定，符號[·]^H表示共軛轉(zhuǎn)置操作。最后，由信號向量x(k)計算信號協(xié)方差矩陣R，即：

其中x^H(k)表示對信號向量x(k)的共軛轉(zhuǎn)置。

第二步：確定第一個正弦信號的頻率區(qū)間

根據(jù)先驗信息確定第一個正弦信號的頻率區(qū)間，并記為Φ₁＝[f₁^l,f₁^r]，其中f₁^l為頻率區(qū)間的下邊界，f₁^r為頻率區(qū)間的上邊界，單位均為Hz。第一個正弦信號的頻率落在頻率區(qū)間Φ₁之中。

第三步：計算自適應(yīng)濾波器的權(quán)系數(shù)向量

根據(jù)第一步得到的信號協(xié)方差矩陣R以及第二步得到的第一個正弦信號的頻率區(qū)間Φ₁，按照下式計算自適應(yīng)濾波器的權(quán)系數(shù)向量h，即：

其中，R^-1表示信號協(xié)方差矩陣R的逆，向量f_c為自適應(yīng)濾波器的中心頻率，f_c∈Φ₁表示f_c的取值范圍在頻率區(qū)間Φ₁內(nèi)。按照上述方式設(shè)計的自適應(yīng)濾波器具有以下特性：當(dāng)自適應(yīng)濾波器的中心頻率不等于第一個正弦信號的頻率時，自適應(yīng)濾波器的幅頻響應(yīng)在頻率區(qū)間內(nèi)存在明顯的零陷，而當(dāng)自適應(yīng)濾波器的中心頻率等于第一個正弦信號的頻率時，自適應(yīng)濾波器的幅頻響應(yīng)在頻率區(qū)間內(nèi)不存在明顯的零陷。

第四步：計算自適應(yīng)濾波器的幅頻響應(yīng)

根據(jù)第三步計算得到的自適應(yīng)濾波器的權(quán)系數(shù)向量h，計算自適應(yīng)濾波器在頻率區(qū)間Φ₁內(nèi)的幅頻響應(yīng)H(f)，即：

H(f)＝|a^H(f)h|,f∈Φ₁

其中，向量|·|表示求復(fù)數(shù)的幅度。自適應(yīng)濾波器在頻率區(qū)間Φ₁內(nèi)的幅頻響應(yīng)H(f)也可以由信號協(xié)方差矩陣R和自適應(yīng)濾波器的中心頻率f_c直接計算得到，即：

第五步：計算信號偽功率譜

首先，根據(jù)第四步計算得到的自適應(yīng)濾波器在頻率區(qū)間Φ₁內(nèi)的幅頻響應(yīng)H(f)，搜索其最大值和最小值，分別記為：

然后，由自適應(yīng)濾波器幅頻響應(yīng)最大值H_max(f_c)和自適應(yīng)濾波器幅頻響應(yīng)最小值H_min(f_c)計算信號偽功率譜P(f_c)，即：

信號偽功率譜P(f_c)也可以由信號協(xié)方差矩陣R和自適應(yīng)濾波器的中心頻率f_c直接計算得到，即：

第六步：計算第一個正弦信號的頻率估計值

根據(jù)第五步計算得到的信號偽功率譜P(f_c)，搜索其在頻率區(qū)間Φ₁內(nèi)的譜峰，并將信號偽功率譜P(f_c)的譜峰對應(yīng)的頻率作為第一個正弦信號的頻率估計值f₁，即：

若在觀測信號中正弦信號的個數(shù)大于1，則首先根據(jù)先驗信息確定其它正弦信號的頻率區(qū)間，然后重復(fù)上述步驟中的第三步至第六步，即可得到其它正弦信號的頻率估計值。

本發(fā)明的有益效果是：在同樣的信噪比及信號采樣點數(shù)等條件下，本發(fā)明方法的頻率估計精度優(yōu)于或接近已有頻率估計方法的頻率估計精度，特別是在信噪比很低、信號采樣點數(shù)很少等情況下，相比于已有的頻率估計方法，本發(fā)明方法能夠在正弦信號頻率附近得到更加明顯的譜峰，從而得到更加精確的頻率估計值。

附圖說明

圖1為本發(fā)明所提供的基于自適應(yīng)濾波器幅頻響應(yīng)的高精度頻率估計方法的流程圖；

圖2為仿真實驗中不同的自適應(yīng)濾波器中心頻率對應(yīng)的自適應(yīng)濾波器幅頻響應(yīng)；

圖3-1為仿真實驗中當(dāng)信噪比SNR＝10dB，信號采樣點數(shù)N＝2000時，F(xiàn)FT頻率估計方法的功率譜。

圖3-2為仿真實驗中當(dāng)信噪比SNR＝10dB，信號采樣點數(shù)N＝2000時，MUSIC頻率估計方法的功率譜。

圖3-3為仿真實驗中當(dāng)信噪比SNR＝10dB，信號采樣點數(shù)N＝2000時，本發(fā)明方法在第一個正弦信號頻率附近的偽功率譜。

圖3-4為仿真實驗中當(dāng)信噪比SNR＝10dB，信號采樣點數(shù)N＝2000時，本發(fā)明方法在第二個正弦信號頻率附近的偽功率譜。

圖4-1為仿真實驗中當(dāng)信噪比SNR＝-20dB，信號采樣點數(shù)N＝500時，F(xiàn)FT頻率估計方法的功率譜。

圖4-2為仿真實驗中當(dāng)信噪比SNR＝-20dB，信號采樣點數(shù)N＝500時，MUSIC頻率估計方法的功率譜。

圖4-3為仿真實驗中當(dāng)信噪比SNR＝-20dB，信號采樣點數(shù)N＝500時，本發(fā)明方法在第一個正弦信號頻率附近的偽功率譜。

圖4-4為仿真實驗中當(dāng)信噪比SNR＝-20dB，信號采樣點數(shù)N＝500時，本發(fā)明方法在第二個正弦信號頻率附近的偽功率譜。

圖5為仿真實驗中在不同的信噪比下已有頻率估計方法和本發(fā)明方法的頻率估計均方根誤差。

圖6為仿真實驗中在不同的信號采樣點數(shù)下已有頻率估計方法和本發(fā)明方法的頻率估計均方根誤差。

具體實施方式

下面結(jié)合附圖對本發(fā)明提供的基于自適應(yīng)濾波器幅頻響應(yīng)的高精度頻率估計方法進(jìn)行詳細(xì)說明。

圖1是本發(fā)明所提供的基于自適應(yīng)濾波器幅頻響應(yīng)的高精度頻率估計方法的流程圖。該流程圖的第一步是計算信號協(xié)方差矩陣，首先，對觀測的連續(xù)時間信號進(jìn)行采樣得到離散時間信號，然后，由離散時間信號計算得到信號協(xié)方差矩陣R。第二步是確定第一個正弦信號的頻率區(qū)間，即由先驗信息確定第一個正弦信號的頻率區(qū)間Φ₁。第三步是計算自適應(yīng)濾波器的權(quán)系數(shù)向量，即由信號協(xié)方差矩陣R和頻率區(qū)間Φ₁計算自適應(yīng)濾波器的權(quán)系數(shù)向量h。第四步是計算自適應(yīng)濾波器的幅頻響應(yīng)，即由自適應(yīng)濾波器的權(quán)系數(shù)向量h計算自適應(yīng)濾波器在頻率區(qū)間Φ₁內(nèi)的幅頻響應(yīng)H(f)。第五步是計算信號偽功率譜，即由自適應(yīng)濾波器的幅頻響應(yīng)在頻率區(qū)間內(nèi)的最小值及最大值計算信號偽功率譜P(f_c)。第六步是計算第一個正弦信號的頻率估計值，即在頻率區(qū)間Φ₁內(nèi)搜索信號偽功率譜P(f_c)的譜峰，并將信號偽功率譜P(f_c)的譜峰對應(yīng)的頻率作為第一個正弦信號的頻率估計值f₁。若剩余未估計的正弦信號個數(shù)大于0，則根據(jù)先驗信息確定下一個正弦信號的頻率區(qū)間，然后重復(fù)執(zhí)行步驟三至步驟六，得到該正弦信號的頻率，以此類推，最終得到所有正弦信號的頻率。

圖2至圖6是利用本發(fā)明方法在具體實施方式下進(jìn)行仿真實驗的處理結(jié)果，其中仿真實驗的基本參數(shù)設(shè)置如下：

觀測信號包含正弦信號和噪聲，其中正弦信號的個數(shù)K＝2，正弦信號的頻率分別為f₁＝150Hz和f₂＝300Hz，噪聲為零均值加性高斯白噪聲。對觀測信號的采樣頻率為1000Hz，信號向量x(k)的長度m取為m＝100。在FFT方法中，設(shè)置FFT運(yùn)算點數(shù)為65536點，在MUSIC方法、Root-MUSIC方法、ESPRIT方法以及本發(fā)明方法中，設(shè)置正弦信號的個數(shù)為2，在本發(fā)明方法中先驗信息還包括：第一個正弦信號的頻率區(qū)間為Φ₁＝[145Hz,155Hz]，第二個正弦信號的頻率區(qū)間為Φ₂＝[295Hz,305Hz]。

圖2是仿真實驗中當(dāng)信噪比SNR＝10dB，信號采樣點數(shù)N＝2000時，不同的自適應(yīng)濾波器中心頻率對應(yīng)的自適應(yīng)濾波器幅頻響應(yīng)。圖2中橫軸表示頻率，單位為Hz，縱軸表示自適應(yīng)濾波器幅頻響應(yīng)的大小。在圖2中，當(dāng)自適應(yīng)濾波器中心頻率f_c＝150Hz時，即自適應(yīng)濾波器中心頻率f_c等于第一個正弦信號的頻率f₁時，自適應(yīng)濾波器幅頻響應(yīng)在頻率區(qū)間Φ₁內(nèi)近似為常數(shù)，且不存在明顯的零陷。當(dāng)自適應(yīng)濾波器中心頻率f_c＝149.9Hz或f_c＝149.8Hz時，即自適應(yīng)濾波器中心頻率f_c不等于第一個正弦信號的頻率f₁，此時自適應(yīng)濾波器幅頻響應(yīng)在頻率區(qū)間Φ₁內(nèi)存在明顯的零陷。

圖3-1至圖3-4是仿真實驗中當(dāng)信噪比SNR＝10dB，信號采樣點數(shù)N＝2000時，F(xiàn)FT頻率估計方法和MUSIC頻率估計方法的功率譜以及本發(fā)明方法在兩個正弦信號頻率附近的偽功率譜。在圖3-1至圖3-4中，橫軸表示頻率，單位為Hz，縱軸表示歸一化功率譜或歸一化偽功率譜，單位為dB。由圖3-1至圖3-4可見，F(xiàn)FT方法、MUSIC方法以及本發(fā)明方法在正弦信號的頻率附近均能得到明顯的譜峰，因此，F(xiàn)FT方法、MUSIC方法以及本發(fā)明方法在信噪比足夠大和信號采樣點數(shù)足夠多的條件下，都能夠得到精確的頻率估計值。

圖4-1至圖4-4是仿真實驗中當(dāng)信噪比SNR＝-20dB，信號采樣點數(shù)N＝500時，F(xiàn)FT頻率估計方法和MUSIC頻率估計方法的功率譜以及本發(fā)明方法在兩個正弦信號頻率附近的偽功率譜。在圖4-1至圖4-4中，橫軸表示頻率，單位為Hz，縱軸表示歸一化功率譜或歸一化偽功率譜，單位為dB。由圖4-1至圖4-4可見，由于信噪比很低且信號采樣點數(shù)很少，由FFT方法計算得到的功率譜沒有得到明顯的譜峰，而由MUSIC方法計算得到的功率譜在錯誤的頻率上得到了譜峰。但是，本發(fā)明方法仍然能夠在正弦信號的頻率附近得到較明顯的譜峰。因此，相對于FFT方法和MUSIC方法等已有頻率估計方法，本發(fā)明方法在信噪比很低和信號采樣點數(shù)很少等情況下能夠取得更加精確的頻率估計。

圖5為仿真實驗中當(dāng)信號采樣點數(shù)N＝500時FFT方法、MUSIC方法、Root-MUSIC方法及ESPRIT方法等已有頻率估計方法和本發(fā)明方法在不同信噪比下的頻率估計均方根誤差。圖5中橫軸表示信噪比，單位為dB，信噪比的變化范圍是-30dB～10dB，縱軸表示頻率估計均方根誤差，單位為Hz。由圖5可見，本發(fā)明方法的頻率估計精度在不同信噪比下均優(yōu)于或接近已有頻率估計方法，特別是當(dāng)信噪比很低時，例如當(dāng)信噪比SNR＜-10dB時，本發(fā)明方法的頻率估計均方根誤差明顯低于已有頻率估計方法。

圖6為仿真實驗中當(dāng)信噪比SNR＝-20dB時FFT方法、MUSIC方法、Root-MUSIC方法及ESPRIT方法等已有頻率估計方法和本發(fā)明方法在不同的信號采樣點數(shù)下的頻率估計均方根誤差。圖6中橫軸表示信號采樣點數(shù)，信號采樣點數(shù)的變化范圍是100～1000，縱軸表示頻率估計均方根誤差，單位為Hz。由圖6可見，本發(fā)明方法的頻率估計精度在不同信號采樣點數(shù)下均優(yōu)于已有頻率估計方法。