時(shí)頻信號(hào)的非參數(shù)加窗稀疏協(xié)方差迭代分析法*

王悅斌，張建秋

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

0 引 言

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

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

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

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

1 時(shí)頻分析方法回顧

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

(1)

(2)

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

1.1 短時(shí)傅里葉變換STFT

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

(3)

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

1.2 Capon時(shí)頻變換

(4)

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

(5)

(6)

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

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

yn=Hx+ηn

(7)

(8)

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

(9)

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

(10)

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

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

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

(11)

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

(12)

(13)

(14)

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

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

(a)矩形窗(a) Rectangle window

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

(15)

(16)

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

(17)

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

(18)

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

(19)

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

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

3.1 超分辨率譜估計(jì)

仿真實(shí)驗(yàn)中的采樣點(diǎn)數(shù)設(shè)為N(N=120)，采樣頻率為1Hz。觀測中待估計(jì)信號(hào)數(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進(jìn)行譜估計(jì)。離散化頻率間隔為0.002Hz，其中Capon數(shù)據(jù)的分段長度L(L=15)，加窗SPICE迭代次數(shù)為J(J=6)。

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

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

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

(b)卡彭(b) Capon

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

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

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

(b)卡彭(b) Capon

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

3.2 多分量信號(hào)的時(shí)頻分析

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

(20)

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

圖4 觀測信號(hào)的真實(shí)時(shí)頻分布圖Fig.4 Ground-truth time-frequency distribution of measured signal

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

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

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

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

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

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

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

4 結(jié) 論

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