時頻信號的非參數加窗稀疏協方差迭代分析法*

王悅斌，張建秋

(復旦大學智慧網絡系統研究中心和電子工程系·上海·200433)

0 引 言

多分量非平穩信號的瞬時幅度和頻率估計在諸多領域(例如導航、制導、振動、語音分析、雷達系統和生物醫學等[1-4])扮演著重要的角色。然而，傳統的譜估計方法只能揭示全局的頻譜分布，無法分析出瞬時幅度和頻率隨時間變化的信號。鑒于此，時頻分析(Time-Frequency Analysis, TFA)方法可用來估計信號時頻兩域的聯合分布。其中最為人所熟知的是短時傅里葉變換(Short-Time Fourier Transform, STFT)[5]，該方法假設非平穩信號在較短的時間窗內是近似平穩的，這樣就可對短時窗內的信號做傅里葉變換，來獲得該時間窗內信號的局部頻率特性。

然而，在對非平穩信號進行STFT分析時，一般要求時間窗內的信號是平穩的。但是當信號具有非常強的非平穩性時，就要求時間窗必須很短，這樣其能量集中度必將下降。為了提高隨時間變化的頻譜的能量集中度，提出了Wigner分布(Wigner Distribution , WD)[6]分析法。WD分析法不需要任何窗函數，因此具有較好的時頻能量集中度，但該方法具有非常強的交叉項效應。此外，在基本S變換的基礎上，也發展出了廣義S變換[7]，它試圖采用多參數的窗函數，來靈活調節窗函數的形狀，進而為克服時頻分辨率單一的缺點提供了一條新途徑。

近幾年，研究者一直致力于改善時頻分析的性能。例如，針對STFT中時間分辨率和頻率分辨率必須進行折衷選擇的缺點，文獻[8]提出了自適應的STFT(Adaptive Short-Time Fourier Transform, ASTFT)算法。該算法能根據信號調頻斜率的變化不斷調節自身的STFT窗長度，從而使信號在調頻斜率較小時具有較高的頻率分辨率，在調頻斜率較大時有較高的時間分辨率。然而，當存在多個復雜的信號分量時，該算法則無法給出不同信號分量要求的不同最優分析窗長度。文獻[9]則提出了利用Capon方法對短時窗內的信號進行頻譜估計，以獲得窗內的局部頻率特性，并隨著窗的向前滑動，最終獲得了信號時頻兩域的聯合分布。

針對目前時頻分析方法普遍存在的時頻分辨率不足的缺點，本文基于傳統的時頻分析框架，給出了一種非參數加窗稀疏協方差迭代分析(Weighted Sparse Iterative Covariance-based Estimation，SPICE)[10]法。該方法首先給出了局部化窗內信號的非參數時頻模型，進而利用加權最小二乘(Weighted Least Square，WLS)法求解優化。將WLS的加權矩陣構建問題轉化為信號的廣義協方差矩陣擬合問題，可利用加窗SPICE方法來獲得短時窗內的局部頻率特性。最后，采用滑動的時間窗函數獲得時頻分布圖。結果表明：提出方法在噪聲抑制和能量集中度等方面均優于文獻報道的方法。

1 時頻分析方法回顧

本文考慮的多分量時變非平穩信號的離散形式可描述為[9]

(1)

(2)

對于形如式(1)的非平穩信號，其統計特性(包括該信號的頻譜特性)都會隨著時間變化而變化，這意味著在對其進行分析和處理時不能將時域和頻域截然分開。本文研究的時頻分析方法就是為了將信號時域和頻域特性結合起來以表征信號，并能對其進行分析和處理。

1.1 短時傅里葉變換STFT

在離散時間n和頻率m上的離散STFT表達式為[5]

(3)

測不準原理表明，寬的時間窗可以獲得高的頻率分辨率，但是卻可能無法檢測到頻率分量的快速變化；窄時間窗雖可以跟蹤到頻率分量的快速變化，卻降低了頻率的分辨率。兩種極端的情況是：當N=1時，SSTFT(n,m)=y(n)；當N等于總的數據長度時，SSTFT(n,m)=DFT{y(n)}，即退化到離散傅里葉變換。因此，如何在保證時間分辨率要求的同時，獲得最高的頻率分辨率，是本文的主要研究內容。

1.2 Capon時頻變換

(4)

式中，gN=[gN(-N/2)，…，gN(N/2-1)]表示自適應數據的濾波器。該最優問題要求的是：盡可能保留頻點m處的能量，而抑制其他頻點的能量。上述約束優化問題借助于拉格朗日乘子法獲得的解為

(5)

(6)

Capon法分辨率雖然較FFT有所提升，但其準確性和分辨率極易受到數據分段的影響[12]，這限制了其分析較短時間窗內信號的能力。

2 非參數化時頻分析法

yn=Hx+ηn

(7)

(8)

式中，W為加權矩陣。該問題的解為[12]

(9)

對于如何構造加權矩陣W，首先要考慮觀測信號的廣義協方差矩陣。假設信號和噪聲相互獨立，并且不同頻率成分的信號不相關，則測量信號yn的協方差矩陣為

(10)

式中，IN×N表示N×N的單位陣。當W=R-1時，式(9)的均方誤差最小[12]。這樣，該問題即轉換為如何估計協方差R。

2.1 加權SPICE 估計協方差

對于如何估計測量信號yn的協方差，SPICE考慮的問題模型如下[13]

(11)

式中，第1項為信號的協方差矩陣，第2項則為噪聲的協方差矩陣。對比式(10)和式(11)，可以發現基于式的模型不僅可以用來實現稀疏譜估計，并且適應于非高斯噪聲的情況。而對于{pm}的估計，SPICE則是通過如下最小化協方差擬合準則而實現[10]

(12)

(13)

(14)

2.2 加窗非參數化時頻分析

對于式(3)的STFT，常用的非矩形窗有高斯窗(Gauss)、漢寧窗(Hanning)和哈明窗(Hamming)等。這些窗函數對于信號頻譜泄漏的大小及頻譜分辨率等方面的影響不同，因此可根據實際需求采用不同的窗函數。例如矩形窗主瓣窄，頻率識別精度高；高斯窗則旁瓣小，幅度識別精度高。

(a)矩形窗(a) Rectangle window

(b) 高斯窗(b) Gaussian window圖1 不同窗內的信號時域圖Fig.1 Signal of time domain in different window

(15)

(16)

由式(16)可以發現，測量信號的協方差矩陣依然具有式(11)的形式，這是因為噪聲協方差依然可以表示為

(17)

此時，非矩形窗下的SPICE模型則變為

(18)

由式(12)～(14)可推導出，任意窗函數下在離散時間n和頻率m上的非參數時頻分析結果為

(19)

式(19)表明：加窗非參數化時頻分析的結果包含了窗函數的信息，這也意味著其給出的時頻分析法不再局限于矩形窗。與式(14)相比，式(19)中的頻率導向矢量gm引入了窗內信號的非平穩性信息，因而在非矩形窗下的式(19)的分析結果更加精確。

3 仿真試驗及結果分析

3.1 超分辨率譜估計

仿真實驗中的采樣點數設為N(N=120)，采樣頻率為1Hz。觀測中待估計信號數目為K(K=4)，頻率位置分別處于f1=0.05Hz、f2=0.06Hz、f3=0.20Hz、f4=0.32Hz，幅度大小依次為1、1、1、0.5。在觀測中加入信噪比(Signal-to-Noise Ratio, SNR )為10dB的高斯白噪聲，分別用FFT、Capon和SPICE進行譜估計。離散化頻率間隔為0.002Hz，其中Capon數據的分段長度L(L=15)，加窗SPICE迭代次數為J(J=6)。

圖2為100次蒙特卡洛仿真下的平均譜估計結果曲線，可以看到，FFT譜估計結果(圖2(a))的主瓣寬度較寬，且對旁瓣的抑制效果較差，無法完全分離2個頻率差異較小的信號分量；Capon(圖2(b))對旁瓣的抑制效果有提升，但是依然無法完全分離出2個頻率差異較小的信號分量；由圖2(c)可以看出，SPICE算法不僅對旁瓣的抑制效果較好，且能夠基本分離出2個頻率差異較小的信號分量。

減小信號觀測數量至N(N=60)，保持其他實驗條件不變，分別采用這3種算法進行譜估計。

(a)傅里葉變換(a) FFT

(b)卡彭(b) Capon

(c) 加窗SPICE(c) Windowed SPICE圖2 多分量信號譜估計結果Fig.2 Spectrum estimations of multi-components signal

由仿真條件可知，理論上最小頻率分辨為1/N(約為0.167)，大于信號分量的最小頻率間隔0.01。圖3為100次蒙特卡洛仿真下的平均譜估計結果曲線，可以看到，FFT和Capon均無法區分出2個頻率差異很小的信號分量；SPICE算法依然可以檢測出4個信號分量，分離效果和抑制旁瓣的能力更好，譜估計結果達到最佳。

(a)傅里葉變換(a) FFT

(b)卡彭(b) Capon

(c) 加窗SPICE(c) Windowed SPICE圖3 超分辨率譜估計結果Fig.3 Super-resolution spectrum estimations

3.2 多分量信號的時頻分析

考慮到多分量時變非平穩信號的一般形式，本文將信號的線性和非線性調頻、信號分量的動態出現和消失等情況聯合起來進行分析，仿真信號的具體形式如下

(20)

仿真信號的采樣頻率為512Hz，信號長度為10s，信噪比為0dB。該測量信號真實的時頻分布如圖4所示，信號分量間出現了一段最小載頻差異為10Hz的區域。分別采用STFT、基于Capon的時頻分析法和本文方法對仿真信號進行分析，如圖3所示。所有算法均采用矩形窗，窗的長度為N(N=64)，窗的移動步長為N/4。離散化頻率間隔為1Hz，其中Capon數據分段長度為L(L=15)，加窗SPICE的迭代次數為J(J=6)。

圖4 觀測信號的真實時頻分布圖Fig.4 Ground-truth time-frequency distribution of measured signal

(a) 短時傅里葉變換(a) STFT

(c)加窗SPICE估計法(c) Windowed SPICE method圖5 矩形窗的時頻估計結果Fig.5 Time-frequency distribution estimations with rectangle window

可以看出，STFT算法(圖5(a))的時頻脊線較粗，頻率分辨率較差，且旁瓣抑制能力較差；基于Capon的時頻譜(圖5(b))的能量集中度較STFT有所提高，但其噪聲抑制能力較差。本文方法的時頻譜如圖5(c)所示，其噪聲抑制能力較好，時頻脊線更細，能量集中度最佳。

保持其他實驗條件不變，將上述試驗中的矩形窗改為高斯窗，得到如圖6的實驗結果。可以看出，STFT算法(圖6(a))的時頻脊線變粗，頻率分辨率較差，且旁瓣抑制能力得到增強；基于Capon的時頻譜(圖6(b))噪聲抑制的能力提高，但能量集中度下降，基本與STFT相當。本文方法的時頻譜如圖6(c)所示，其噪聲抑制能力較好，時頻脊線更細，能量集中度依舊最高。

(a)短時傅里葉變換(a) STFT

(b)卡彭估計法(b) Capon method

(c)加窗SPICE估計法(c) Windowed SPICE method圖6 高斯窗的時頻估計結果Fig.6 Time-frequency distribution estimations with Gaussian window

4 結 論

針對多分量非平穩時變信號，本文提出了一種非參數化稀疏協方差迭代的時頻分析法。提出方法同樣基于傳統時頻分析的框架，即利用局部化矩形窗函數來分析信號的局部頻率特性。對于窗內信號，本文采用非參數化時頻模型，并借助WLS算法，獲得了時頻譜估計結果。其中，針對WLS加權矩陣的構建，本文利用加權SPICE對廣義噪聲協方差矩陣進行了迭代估計。此外，本文考慮了非矩形窗時的非參數化時頻模型的修正問題，適用于各種類型的窗信號。提出方法對于噪聲譜有很好的抑制效果，同時兼具更高的頻率分辨率，這一優越性在仿真對比試驗中得到了充分體現。