王悅斌，張建秋

(復(fù)旦大學(xué)智慧網(wǎng)絡(luò)系統(tǒng)研究中心和電子工程系·上?！?00433)

0 引 言

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

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

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

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

1 時頻分析方法回顧

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

(1)

(2)

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

1.1 短時傅里葉變換STFT

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

(3)

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

1.2 Capon時頻變換

(4)

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

(5)

(6)

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

2 非參數(shù)化時頻分析法

yn=Hx+ηn

(7)

(8)

式中，W為加權(quán)矩陣。該問題的解為[12]

(9)

對于如何構(gòu)造加權(quán)矩陣W，首先要考慮觀測信號的廣義協(xié)方差矩陣。假設(shè)信號和噪聲相互獨立，并且不同頻率成分的信號不相關(guān)，則測量信號yn的協(xié)方差矩陣為

(10)

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

2.1 加權(quán)SPICE 估計協(xié)方差

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

(11)

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

(12)

(13)

(14)

2.2 加窗非參數(shù)化時頻分析

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

(a)矩形窗(a) Rectangle window

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

(15)

(16)

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

(17)

此時，非矩形窗下的SPICE模型則變?yōu)?/p>

(18)

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

(19)

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

3 仿真試驗及結(jié)果分析

3.1 超分辨率譜估計

仿真實驗中的采樣點數(shù)設(shè)為N(N=120)，采樣頻率為1Hz。觀測中待估計信號數(shù)目為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數(shù)據(jù)的分段長度L(L=15)，加窗SPICE迭代次數(shù)為J(J=6)。

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

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

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

(b)卡彭(b) Capon

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

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

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

(b)卡彭(b) Capon

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

3.2 多分量信號的時頻分析

考慮到多分量時變非平穩(wěn)信號的一般形式，本文將信號的線性和非線性調(diào)頻、信號分量的動態(tài)出現(xiàn)和消失等情況聯(lián)合起來進行分析，仿真信號的具體形式如下

(20)

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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