數字信號處理大實驗2 利用FFT估計信號的頻率

3.1 實驗目的

3.2 實驗內容與要求

3.3 實驗原理

3.3.1 基于時域求導-頻域乘法的n階導數積分法

3.3.2 基于頻域卷積的雙/多譜線插值法

3.3.3 基于譜峰和滑動平均的多譜線綜合插值方法

3.3.4 基于相鄰顯著譜線的滑動平均綜合插值方法

3.3.5 基于（2）的多步小量時移平均計算法（頻域操作）

3.3.6 基于（5）的多步小量頻移平均計算法（時域操作）

3.4 實驗結果

3.4.1 真實頻率估計

3.4.2 探究采樣點數N對于估計的相對精度（去量綱化）的影響

3.4.3 探究真實頻率對于估計的相對精度的影響

3.4.4 探究初相（采樣起始點）對于估計的相對精度的影響

3.4.5 說明、總結與分析

參考文獻

3.1 實驗目的

了解電磁通信等實際工程中的一些基本問題及解決思路；
培養查閱資料解決問題的能力；
掌握頻偏的概念，加深對于DFT、DFT的常用性質（時移/頻移等）的理解；
設計具備一定準確性、魯棒性的算法，估計單頻正弦型信號的真實頻率；
探究采樣率固定時，采樣區間長度、起始采樣點、被測信號真實頻率對于計算精確度和穩定性的影響。

3.2 實驗內容與要求

在無噪聲干擾條件下，用8000 Hz采樣率對某387Hz頻率正弦信號進行采樣，采集了77個樣點作為原始信號。

1、在知網、IEEE網站等查找頻率估計方法，利用至少兩種頻率估計方法編程求給定信號的頻率；

2、無噪聲、20dB、10dB、5dB、0dB、-5dB，每個信噪比通過至少100次的平均計算，得到頻偏的方均根；

3、所有算法使用相同點數的FFT [不大于2048點]；

4、繪出信噪比（對數）-估計均方誤差圖。

3.3 實驗原理

??????

本仿真實驗筆者一共實現了6種算法（實際只有5種，不包括（2）），包括：

基于時域求導-頻域乘法的n階導數積分法;
Quinn改進-實部虛部比值與綜合權重的雙譜線插值法；
基于譜峰和滑動平均的多譜線綜合插值方法；
基于相鄰顯著譜線的滑動平均綜合插值方法；
基于（2）的多步小量時移平均計算法（頻域操作）；
基于（5）的多步小量頻移平均計算法（時域操作）。

3.3.1 基于時域求導-頻域乘法的n階導數積分法

Figure 14??????? 華南理工大學電信學院寧更新老師課件 Chapter5-v3 P28

代碼實現

Figure 15??????? 3.3.1 n階導數-頻域范圍內求和的代碼實現

3.3.2 基于頻域卷積的雙/多譜線插值法

長度為N點的矩形窗的離散傅里葉變換

對單頻f0正弦型信號時域采樣，采樣序列為y[n]=A*cos(2*π*f0/fs*n+φ)，假定f0/fs是有理數，fs>2*f0；周期為N=m*fs/f0 (?m，s.t. ?m*fs/f0∈Z),該信號的DTFT——am=A/2*ejφ, aN-m=A/2*e-jφ,ak=0 for 0<=k<=N-1 and k≠m,N-m

使用N點矩形窗對正弦信號進行截取，需要保證頻譜分辨率一致，此時正弦信號的DTFT轉化為同樣點數的DFT，am和aN-m的序列號和N有關，可能成為小數；頻譜為DFT的圓周卷積，對N取余，因此得到的譜線的真實譜峰位置索引經常不是整數，這也是頻差δ的來歷；

由于實信號幅度頻譜的正負頻率對稱性，算法截取0到pi的頻率區間，實際由于雙峰卷積造成混疊，從而影響理論插值計算的適用性和準確性。

Figure 16??????? 由于正負頻率的對稱性造成的混疊，fs在4倍f0附近時雙卷積帶來的頻差計算誤差會較小

筆者的算法主要依據為理論插值計算，并在Quinn方法和文獻[1]的基礎上，通過引入虛部比值和三角函數權重的方式減小頻偏估計的誤差。截斷后信號s的DFT實部虛部分別如下圖所示。

Quinn雙譜線實部比值插值方法：if δ1≥0 and δ2≥0,δ=δ1;elseif δ1<0 and δ2<0,δ=δ2

基于虛部比值的雙譜線插值方法：if δ1≥0 and δ2≥0,δ=δ1;elseif δ1<0 and δ2<0,δ=δ2

Figure 17??????? 公式推導[1]（含筆者的條件修正）、算法流程（部分）[1]

Y(m) 、Y(m-1) 和 Y(m+1) 實部的大小與峰值譜線的相位Φ(m)有關，該算法在峰值譜線的相位Φ(m)接近于0或π時，頻率校正精度較高，但當峰值譜線的相位Φ(m)接近于 ±π/2 時，信號頻譜分布與FFT系數的關系如圖1所示，FFT系數實部序列幅度很小，接近于零，在噪聲干擾下，根據實部序列索引的峰值譜線位置容易出錯，計算的頻偏δ誤差也較大，此時頻率估計精度較低[1]。

代碼實現（部分）

Figure 18??????? 3.3.2算法使用的三角函數綜合權重法，角度是幅度譜峰的極角（相位）

3.3.3 基于譜峰和滑動平均的多譜線綜合插值方法

Figure 19 滑動平均原理[2]

根據截斷后信號s的DFT表達式，可以推導出幅度譜譜峰和索引距離譜峰索引為k的譜線的實部/虛部比值——Figure18（文獻[2]式（5）），在以譜峰索引為中心且左右對稱、索引半徑為m的區間內，i表示當前的絕對值距離，則從1到m按照3.3.2的方法各可以計算并篩選出一個頻偏估計δi ，通過系數乘法配湊進行等權重估計，得到δ。

代碼實現（部分）

Figure 20??????? 3.3.3 算法的代碼實現（部分）

3.3.4 基于相鄰顯著譜線的滑動平均綜合插值方法

在以譜峰索引為中心且左右對稱、索引半徑為m的區間內，i表示當前的距離，則從-m到m按照3.3.2的方法各可以計算并篩選出一個頻偏估計δi ，通過系數加法配湊進行等權重估計，得到δ。

Figure 21??????? 算法3.3.4 關鍵代碼

3.3.5 基于（2）的多步小量時移平均計算法（頻域操作）

在頻域對于每一個點的DFT系數乘于ej2kπt0/N,t0是小時移量，本實驗設置為一個采樣周期（可以是小數倍的采樣周期，因為是頻域操作）；

這種算法不影響譜峰的位置，但是會影響信號的初相，因此可以一定程度上削減初相對于頻偏估計的影響。

Figure 22??????? 算法3.3.5（多步小量時移平均計算法）核心代碼

3.3.6 基于（5）的多步小量頻移平均計算法（時域操作）

在時域對截斷后信號s進行多步小量頻移，頻移量來自3.3.5算法的頻偏估計，通過凈頻偏存儲、多次自迭代計算，理論上頻偏量會越來越小，突破閾值后可以視作較準確的估計，如果在給定最大迭代次數內未能突破閾值，使用所有迭代過程中的平均頻偏估計值。

Figure 23??????? 算法3.3.6（多步小量時移x多步小量頻移）核心代碼實現

3.4 實驗結果

3.4.1 真實頻率估計

Figure 24 截斷后信號s，N=77

??? 真實信號頻率387Hz，幅度約為2.7, 初相約為137.8°；

??? 采樣頻率8KHz，頻率精度8000/77=103.90Hz。

??? 下面是部分噪聲強度配比以及部分算法的估計頻率的數值。

未摻雜高斯零均值白噪聲：

（1）PMA_ri_comb_interpolation（3.3.3）：324.0058Hz；

（2）Ts_two_p_ri_comb_interpolation (3.3.5):371.8061Hz。

Figure 25 ?????? 5種算法的準確度和門限測試，門限基本上在10dB左右（N=77，縱坐標為4lg（δRMSE））

在3.4.1的情況中，3.3.5小量時移算法效果最佳，信噪比SNR=30dB時去量綱化真實頻偏0.1449。

3.4.2 探究采樣點數N對于估計的相對精度（去量綱化）的影響

Figure 26 N=256，縱坐標4lg（δRMSE）- 2，初相、幅度同原信號

Figure 27??????? N=512，縱坐標4lg（δRMSE）- 2，初相、幅度同原信號

在3.4.2的情況中，3.3.3譜峰滑動平均算法效果最佳，如果把0.5的相對精度（絕對精度8000/256/2=15.625Hz）作為最低合格指標，那么對應到此處的對數坐標是-3.2，N=256、PMA的門限在-5dB左右，N=512、PMA門限在5dB左右；N=256、信噪比SNR=30dB時去量綱化真實頻偏0.136；N=512時，信噪比SNR=30dB時小量時移平均計算法精度稍高，去量綱化真實頻偏0.111。