動態(tài)時頻譜分析、探測和跟蹤的隨機有限集法

王悅斌，蔣景飛，張建秋，＊

1.復旦大學 智慧網(wǎng)絡系統(tǒng)研究中心和電子工程系，上海 200433 2.電子信息控制重點實驗室，成都 610036

非平穩(wěn)信號的分析和參數(shù)估計，在雷達［1］、語音分析［2］和生物醫(yī)學［3］等領域中正發(fā)揮著越來越重要的作用。然而，傳統(tǒng)的譜估計算法只能分析出全局譜，無法揭示信號頻率的局部變化。為了獲取不同時刻信號分量的存在與否，以及各頻率成分隨時間變化情況，時頻分布（Time-Frequency Distribution，TFD）的分析方法應運而生。其中短時 傅 里 葉 變 換 （Short-Time Fourier Transform，STFT）、Wigner-Ville分布和S 變換［4］是3種典型的時頻分析算法。這些方法簡單易懂，不需要額外的參數(shù)模型，應用十分方便。

可是，對多分量非穩(wěn)態(tài)信號來說，這些傳統(tǒng)算法在頻率分辨率，非線性調(diào)頻特性，以及動態(tài)時頻譜（據(jù)文獻［5］的定義：若時頻譜中的信號分量隨著時間動態(tài)地出現(xiàn)和／或消失，則稱其為“動態(tài)時頻譜”）的分析等方面還存在局限性。例如，Wigner-Ville分布對于單分量信號具有較好的時頻分析能力，而對于多分量信號來說，則會引入較強的交叉項［6］；S變換雖然具有可變的時頻分辨率，但由于窗寬度是頻率的函數(shù)，因而在不同的信號分析中，會表現(xiàn)出不同的能量集中度，從而限制了應用范圍［7］。

為了提高對非平穩(wěn)信號的時頻表示（Time-Frequency Representations，TFRs）能力，近年來，人們研究出了許多新的算法。其中，廣義S變換（Generalized S-Transform，GST）［8］通過多參數(shù)的窗函數(shù)，為克服時頻分析方法所存在時頻分辨率單一的缺點提供了一條新途徑。然而，由于依賴于窗函數(shù)的選擇，GST無法與不同類型的信號相適應。不同于GST，文獻［9］提出的自適應 的 STFT （Adaptive Short-Time Fourier Transform，ASTFT）算法，可自適應非平穩(wěn)信號特性，即據(jù)信號調(diào)頻斜率的變化來不斷調(diào)節(jié)自身STFT窗的長度。可是，對于具有多個分量且調(diào)頻特性較復雜的時頻信號，該算法不能獲得較為可靠的調(diào)頻斜率變化信息。匹配解調(diào)變換（Matching Demodulation Transform，MDT）［10］可以獲得較高的瞬時頻率估計精度，并且能有效提高時頻表示的能量集中度，但該算法需要預先給定信號的全局參數(shù)模型，而大多數(shù)情況下瞬時頻率的模型是未知的，因此該算法實用性較低。文獻［11］提出了不需要信號全局參數(shù)模型的變分非線性調(diào)頻模式分解（Variational Nonlinear Chirp Mode Decomposition，VNCMD）算法。該算法對于非線性調(diào)頻信號的時頻分析，在具有較好性能的同時，也具有較高的時頻分辨率，并且可以解決信號頻率分量交叉問題。然而，該算法需要預知信號分量數(shù)目，因而難以解決動態(tài)出現(xiàn)和消失多分量時頻信號的分析問題。為解決這一問題，文獻［12］提出一種基于Capon譜估計的多分量信號檢測和追蹤算法。該算法首先利用Capon譜估計得到時頻譜的粗估計，進而估計出噪聲譜，再進行信號峰值檢測，最后進行匹配追蹤，獲得了去噪效果較好的時頻分析結(jié)果。然而，該算法需要依靠全局譜判斷信號分量的存在與否，因而實時性較差。最近，文獻［5］提出一種利用狄利克雷混合模型來描述和估計動態(tài)時頻譜的算法，進而給出了一種基于自適應狄利克雷混合模型的Rao-Blackwellised粒子濾波（Adaptive Dirichlet Process Mixture Model-based Rao-Blackwellised Particle Filter，ADPM-RBPF）算法。該算法實時性較好，提高了低信噪比下的時頻譜估計的精度。可是，由于該算法利用STFT得到的原始時頻譜作為觀測，時頻分辨率較差。

針對動態(tài)出現(xiàn)和消失多分量非平穩(wěn)信號，本文提出了一種基于隨機有限集的時頻分析法。該算法利用時頻變換在每個時刻獲得的信號分量幅度的局部極值為測量的隨機有限集，然后將時頻信號各分量幅度和頻率的變化描述為多項式預測模型［13］，進而為每個信號分量建立起了一種新的線性高斯狀態(tài)空間模型。這樣，本文就將動態(tài)出現(xiàn)和消失的時頻譜分析、探測和跟蹤問題，歸納為可利用隨機有限集進行多目標追蹤的問題。借助于高斯混合概率假設密度（Gaussian Mixture Probability Hypothesis Density，GM-PHD）濾波器［14］，本文算法可以分析動態(tài)的時頻譜。對于GM-PHD濾波器，初始權(quán)重的賦值和權(quán)重的更新尤為關鍵。據(jù)時頻分布的特點，本文提出了一種初始權(quán)重賦值算法，結(jié)合提出的譜分量幅度和頻率的聯(lián)合似然函數(shù)，進而解決了權(quán)重的更新問題。分析和數(shù)值仿真結(jié)果均證明：所提算法可有效提升動態(tài)時頻譜的跟蹤精度，對微弱時頻譜分量的探測能力，以及對載頻差異的分析能力。

1 時頻分析的隨機有限集模型

對于一個由多個調(diào)頻調(diào)幅分量構(gòu)成的時頻信號，為了描述其分量隨時間動態(tài)出現(xiàn)和消失的情況，采用Kn來表示隨離散時間n變化的信號分量數(shù)目，這樣本文考慮的時頻信號的離散形式可描述為［5，12］

式中：yn為信號的觀測值；akn和fkn分別為第k個分量信號的幅度和頻率函數(shù)；θk0為第k個分量信號的初始相位；fs為采樣頻率；ηn為觀測噪聲。對時頻信號式（1）進行STFT時頻變換，可以得到信號的時頻分布為［15］

式中：n和τ分別為離散的時間和頻率；N 為STFT短時窗的長度；ωm為窗函數(shù)。由于時頻分布TFDN（n，τ）反映的是信號在（n，τ）時頻鄰域內(nèi)的能量，那么任意時刻信號分量的能量，就必定會集中在該時刻的瞬時頻率（Instantaneous Frequency，IF）附近。如果視每一時刻時頻分布TFDN（n，τ）的局部極大值即式（3）中的Zn為測量值，這意味著任意時刻的信號分量必定包含在該時刻的Zn中。這是因為任意時刻時頻分布TFDN（n，τ）的能量是信號能量和噪聲能量的疊加，對高斯白噪聲而言，其理想時頻分布的能量是某一常數(shù)，這樣當一個常數(shù)加上任一小的信號能量時，都將為高斯白噪聲的理想時頻分布能量增加一個局部極值。

考慮到微弱信號分量的幅值可能很小，而觀測噪聲可能是非高斯的，以及高斯觀測噪聲的樣本可能不是足夠大，因此本文將所有檢測到的局部極值均視為測量，如圖1所示。由于觀測噪聲的存在，以及時頻信號的分量可能隨機出現(xiàn)或消失，那么每一時刻的局部極值也將是隨機的。這樣，可把Zn視為一個測量隨機有限集（RFS）集合［16－17］：

式中：Jn為n時刻局部極值（測量值）的個數(shù)；zjn＝［］T為n時刻第j個局部極值幅度和頻率的測量值集合，其中j＝1，2，…，Jn；Ez為多分量時頻信號的測量空間。因此，n時刻的測量隨機集可以進一步表示為

式中：ηn為n時刻源于觀測噪聲和／或雜波的測量RFS；Θn（x）為n時刻源于時頻信號多分量狀態(tài)x的測量RFS；∪為并運算。視信號分量為目標，則多分量信號的目標狀態(tài)RFS就可描述為［16-17］

式中：Kn為n時刻目標個數(shù)（時頻信號分量的個數(shù)）；xkn為n時刻第k個目標的狀態(tài)，其中k＝1，2，…，Kn；Ex為多目標狀態(tài)空間。n時刻的多目標狀態(tài)隨機集可以進一步表示為

式中：Sn｜n－1（xn－1）為n－1時刻的目標狀態(tài)xn－1在n時刻的存活RFS；Γn為n時刻新生分量目標RFS。這樣，可推導出基于RFS的全體多目標聯(lián)合后驗概率密度的貝葉斯遞推公式［16－17］：

圖1 峰值檢測示意圖Fig.1 Diagram of peak detection

式中：Fn｜n－1（·｜·）為 多 目 標 轉(zhuǎn) 移 概 率 密 度；gn（·｜·）為多目標似然函數(shù)；Z1：n＝ ｛Z1，Z2，…，Zn｝為前n 個時刻所有的測量 RFS；pn｜n－1（·｜Z1：n－1）為多目標預測概率密度；pn｜n（·｜Z1：n）為多目標聯(lián)合后驗概率密度。

式（8）和式（9）的遞推公式需要在多目標狀態(tài)空間上計算集積分，計算量非常大，在實際計算中一般難以實現(xiàn)。在線性高斯條件下，基于RFS的GM-PHD，通過把目標的概率假設密度（Probability Hypothesis Density，PHD）近似為混合高斯形式，可以得到遞推形式式（8）和式（9）的閉合解［14］。針對式（4）和式（6）的 RFS變量，本文將為其提出一種新的線性高斯狀態(tài)空間模型，使得本文可以利用GM-PHD濾波計算目標數(shù)量并估計各個目標的狀態(tài)，進而得到了一種可分析動態(tài)出現(xiàn)和／或消失多分量時頻信號的新復算法。

2 多分量時頻信號分析

2.1 時頻分量信號的狀態(tài)空間模型

考慮到式（1）中信號分量的幅度akn和頻率fkn，既可能是線性，也可能是非線性的，為了使討論具有一般性，本文假設式（1）中的akn和fkn都是非線性函數(shù)，此時線性就可認為是非線性的特例。對于任一非線性函數(shù)，Weierstrass逼近定理［18］可知：任一有界閉區(qū)間上的連續(xù)函數(shù)都可以由該區(qū)間內(nèi)的多項式以任意精度逼近。這意味著對任意連續(xù)函數(shù)，如果將其閉區(qū)間進行分割，且讓每個分割出的小區(qū)間足夠小，以致在每一小區(qū)間內(nèi)，該函數(shù)都可以用一個低階多項式（例如一階多項式）以任一精度逼近。據(jù)這個原理，假設信號分量的瞬時頻率fkn，在每一分割的小區(qū)間內(nèi)可由L階多項式描述：

式中：p（l）為多項式的系數(shù)。如果想用前M個時刻的頻率值來預測當前時刻的頻率，即

顯然，式（11）是以h（m）（m＝1，2，…，M）為系數(shù)的有限沖激響應（FIR）濾波器。將式（10）代入式（11），消去多項式系數(shù)p（l），可得

由式（12）無法唯一確定系數(shù)h（m）。考慮到信號的觀測通常都是存在噪聲的，這樣就需要濾波器的噪聲增益：

具有最小值。在式（12）的約束條件下，利用拉格朗日法就可求得式（13）的最優(yōu)解為［13］：

1）當L＝1時

2）當L＝2時

式中：m＝1，2，…，M，其中M 表示濾波器的長度。當L為其他值時，計算結(jié)果詳見文獻［13］。

上面的分析表明：多項式預測模型式（11）與FIR預測濾波器在數(shù)學上嚴格等效。據(jù)上述多項式的預測模型，就可對時頻信號分量的瞬時幅度和頻率函數(shù)進行建模。定義第k個時頻信號分量在n時刻的狀態(tài)向量為

式（17）也可視作零階多項式的預測模型，即用前一時刻的幅度值來預測當前時刻的幅度值。根據(jù)式（11）和式（17）的模型描述，式（16）的狀態(tài)方程就可以表示為

有了式（18）幅度和頻率的聯(lián)合狀態(tài)方程后，還需要建立信號分量幅度和頻率作為聯(lián)合參數(shù)的測量方程。假設第k個信號分量的測量為zn＝［za，n，zf，n］T，其中za，n和zf，n分別為幅度和頻率的測量值，則式（16）的測量方程可以表示為

式中：σn以0均值的加性噪聲來描述式（3）的不確定性，噪聲方差為R，觀測矩陣為H ＝

在式（18）和式（20）所描述多分量時頻信號的分析中，必須確定式（20）的測量zn屬于式（1）中哪一個信號分量或雜波，才能根據(jù)式（18）和式（20）的狀態(tài)空間模型進行時頻分析，2.2節(jié)將討論。

2.2 GM-PHD濾波的時頻分析法

通過時頻變換，本文將非線性觀測方程式（1）轉(zhuǎn)換成了式（20）的線性觀測方程，而2.1節(jié)對線性／非線性調(diào)頻函數(shù)，通過多項式預測方法，則將其變換成了線性的狀態(tài)方程式（18）。這樣式（18）和式（20）構(gòu)成的時頻分析狀態(tài)空間模型是線性的。此外，由于本文假設噪聲是高斯的，因此采用了GM-PHD濾波器。當式（1）觀測噪聲是非高斯時，則可以采用SMC-PHD［20］等濾波器，但建立的時頻分析狀態(tài)空間模型式（18）和式（20）有效，且將簡化SMC-PHD等濾波器的計算復雜度。

2.2.1 算法預測

假設n－1時刻時頻分量的后驗概率密度是高斯混合形式［14］

式中：Kn－1為n－1時刻的時頻分量個數(shù)；和分別為分量k在n－1時刻的狀態(tài)均值和方差；為分量k的權(quán)重（表示該分量為真實分量的概率）；Ν（x；a，b）為隨機矢量x服從均值為a、方差為b的正態(tài)分布。n時刻已存在分量的預測概率密度υS，n｜n－1（x）為

式中：pS，n為時頻分量存活下來的概率［14］，分量k的預測均值和預測方差可通過 Kalman預測得到。注意，n－1時刻不屬于任何時頻分量的測量值將被假設為n時刻的新生時頻分量。假設時頻分量的預測概率密度也是高斯混合形式：

式中：Kγ，n為n時刻假設時頻分量的個數(shù)。假設時頻分量k的初始權(quán)重ωkγ，n的賦值方式需要據(jù)應用場景來確定。在時頻分析中，考慮到測量值的瞬時幅度越大，其屬于真實時頻分量的概率就越大，而在信噪比高的環(huán)境下，雜波較少，那么測量值來自于真實時頻分量的概率也越大。綜合上述兩個因素，本文給出權(quán)重ωkγ，n的初始值為

2.2.2 算法更新

獲得式（1）的觀測后，對N 個觀測值進行時頻變換，例如STFT，時頻變換后的結(jié)果通過式（3）獲得n時刻測量值后，需要對權(quán)重進行相應的更新。對于每個測量，權(quán)重）根據(jù)貝葉斯公式更新為［14］

式中：pD，n為檢測到時頻分量的概率；）為雜波在上的概率密度，一般取常數(shù)1／V，其中V是測量空間的大小；分別為來自時頻分量k測量的預測均值和方差；Kn｜n－1＝Kn－1＋Kγ，n為n 時刻預測的時頻分量個數(shù)。當假設時頻分量的頻率和幅度相互獨立時，式（25）中的似然函數(shù)可以寫為

事實上，測量zjn只可能來自某一個時頻分量或者來自雜波。這樣，為了降低計算復雜度，本文根據(jù)權(quán)重的大小將測量與時頻分量進行關聯(lián)式（27）表明：當n時刻所有時頻分量的權(quán)重都更新后，其中沒有與目標關聯(lián)過的測量值，將假設為下一時刻的新生時頻分量。然而由于雜波的存在，必須要有一個判別準則來決定假設時頻分量是否來自真實的新時頻分量，還是來自雜波。由于雜波能量較為分散，而時頻分量的能量相對較集中，而據(jù)式（25）可知，真實時頻分量的權(quán)重將接近于1，而虛假時頻分量（雜波）的更新權(quán)重將接近于0，應該舍去。這樣本文的判別準則如下：

1）若ωkn＞0.5，該假設分量為真實分量。

2）若ωkn＞1／V，該分量為假設分量。

3）否則，該假設分量為雜波（該準則也同樣適用于判斷時頻分量的消失）。

3 時頻分析的多目標追蹤算法

步驟1 對存在的時頻分量k進行Kalman預測，其中k＝1，2，…，Kn－1。

式中：Q為過程噪聲協(xié)方差矩陣。

步驟2 由式（27）可知，對在n－1時刻判別為來自雜波的Kγ，n個測量視為n時刻的假設時頻分量，并據(jù)式（24）賦予相應的初始權(quán)重。

步驟3 對存在時頻分量或假設時頻分量k的測量進行預測，其中k＝1，2，…，Kn－1＋ Kγ，n－1。

式中：R為測量噪聲方差。

步驟4 對于給定的zjn，其中j＝1，2，…，Jn，式（27）將該測量值與時頻分量相關聯(lián)，并獲得該分量權(quán)重的更新值；更新時頻分量的狀態(tài)方差。

1）若該測量由某一時頻分量k產(chǎn)生，更新該時頻分量的狀態(tài)均值mkn｜n

步驟5 根據(jù)判別準則，刪除權(quán)重較小的分量，接受權(quán)重較大的假設時頻分量為真實分量。

4 仿真實驗

為了使仿真實驗更具說服力，本文采用了文獻［5，12］的仿真實驗思路來設計動態(tài)時頻譜。仿真實驗分為3部分：① 仿真中的多分量線性調(diào)頻信號來自文獻［5］；② 仿真中包含微弱信號和非線性調(diào)頻信號的情況參考文獻［12］；③ 主要考慮載頻差異較小信號的時頻分析，以借此來說明所提算法的高分辨率性能。

4.1 多分量線性調(diào)頻信號

第1組仿真信號包含10個線性調(diào)頻分量，形式為

式中：αk和βk分別為第k個信號分量的初始頻率和調(diào)頻斜率；Ik（t）為第k個信號分量是否存在的指示函數(shù)，其形式為

加性復高斯白噪聲η（t）的方差為σ2η，信噪比（SNR）定義為［21］

仿真信號長度為100s，采樣頻率為512Hz，信噪比為－8dB。圖2分別給出了信號進行STFT、Capon 變 換［12］和迭 代 自 適 應 （Iterative Adaptive Approach，IAA）譜估計（見附錄 A）的時頻分布圖，其中STFT、Capon和IAA均采用矩形窗，窗長N＝128，窗的移動步長為N。對于圖2的含噪時頻圖，分別采用ADPM-RBPF算法［5］、基于 Capon的峰值檢測法［12］和本文算法進行分析。本文算法中狀態(tài)空間模型參數(shù)為L＝2，M ＝3，σa＝10－2，σfm＝10－4，其中m＝1，2，3。時頻分量存活下來的概率pS，n＝0.9，時頻分量被檢測到的概率pD，n＝0.9，S型曲線的陡峭程度參數(shù)δ＝20，權(quán)重因子ρ1＝2ρ2。

圖2 線性調(diào)頻信號時頻圖Fig.2 Time-frequency spectra of chirp signal

從圖3（a）可以看到，在低信噪比下，ADPMRBPF算法不能很好地去除雜波；從圖3（b）可以看到，基于Capon的峰值檢測法可以更好地去除雜波，但是該算法必須依靠更多的譜來判斷真實的信號分量和雜波［12］，因此實時性差。在相同的窗長下，本文算法的實時性與ADPM-RBPF相當，均優(yōu)于基于Capon的峰值檢測法。而圖3（c）和圖3（d）則表明：本文算法明顯優(yōu)于ADPM-RBPF算法，當短時譜是利用IAA估計得到時，本文算法效果最佳。

單時頻分量背景下的均方根誤差已不再適用于對本文算法的性能評價，因為需要考慮狀態(tài)估計值和真實值之間的對應關系。為此，本文利用OSPA（Optimal SubPattern Assignment）距離［22］來評價多分量估計的誤差。對于瞬時頻率真實值珚F＝［珚f1，珚f2，…，珚f珡K］和估計值＾F＝［＾f1，＾f2，…，＾fK＾］，當珡K≤＾K時，OSPA距離定義為［22］

式中：

式中：參數(shù)c為目標狀態(tài)估計誤差的閾值，同時可調(diào)節(jié)集合勢的估計誤差比重。本文采用p＝2階，c＝0.1的評價指標。

圖4給出了3種算法時頻譜估計的平均OSPA曲線，每種算法均進行100次蒙特卡羅仿真。通過圖4可以看到，本文算法瞬時頻率的估計精度優(yōu)于ADPM-RBPF算法和基于Capon的峰值檢測法，在低信噪比（SNR＜－6dB）下，本文算法優(yōu)勢更加明顯。這是因為本文算法為調(diào)頻函數(shù)建立一個可靠的狀態(tài)轉(zhuǎn)移模型（18），這個新模型使得在低信噪比下，本文算法能依據(jù)該模型魯棒地跟蹤時頻分量。需要強調(diào)的是，基于IAA短時譜的本文算法估計性能最佳，這是因為IAA算法優(yōu)于STFT算法。

圖3 線性調(diào)頻信號的時頻重構(gòu)圖Fig.3 Time-frequency reconstruction spectra of chirp signal

圖4 IF估計誤差Fig.4 IF estimation error

4.2 復雜多分量時頻信號

實際中，信號分量幅度大小往往不同，調(diào)頻方式多樣。為了驗證本文算法在各種復雜情況下的魯棒性和有效性，本文將非線性調(diào)頻特性、信號分量的動態(tài)出現(xiàn)和消失以及微弱信號分量都結(jié)合起來進行分析。第2組仿真信號形式為

仿真信號長度為100s，采樣頻率為512Hz，信噪比為0dB。圖5給出了信號進行STFT變換、Capon變換和IAA譜估計的時頻分布圖。從圖中可以看出，信號分量k＝4較強，而信號分量k＝3非常微弱，其局部信噪比相當于－10.46dB，檢測難度變大。

圖6給出了本文算法與對比算法的時頻重構(gòu)圖，可以發(fā)現(xiàn)，ADPM-RBPF算法（見圖6（a））不能完整重構(gòu)出信號分量k＝3，這是因為該算法需根據(jù)功率譜大小來進行隨機采樣獲得測量值，從而導致微弱信號分量被檢測到概率很小。基于Capon的峰值檢測法（見圖6（b））在微弱信號以及調(diào)頻斜率較大區(qū)域均出現(xiàn)漏檢情況，這是因為該算法采用的假設檢驗在信號能量較低區(qū)域（低信噪比區(qū)域）漏報率高的緣故。本文算法則可利用建立的預測模型式（18）來預測非線性調(diào)頻分量在下一個采樣時刻可能的近似值，同時也考慮到了時頻分量未被檢測到的概率（見式（25）），因此不論基于STFT或IAA短時譜，均可以完整重構(gòu)出全部4個信號分量。

圖7給出了ADPM-RBPF算法、基于Capon的峰值檢測法，以及本文算法估計出的信號分量數(shù)目，每種算法均進行100次蒙特卡羅仿真。當微弱信號分量存在時，ADPM-RBPF算法估計誤差較大；而當調(diào)頻斜率較大或微弱信號分量存在時，基于Capon的峰值檢測法對于信號分量數(shù)目的估計誤差較大。本文算法分析結(jié)果較為理想，這是因為本文算法考慮到了目標未被檢測到的情況。

4.3 載頻差異較小的時頻信號

為了獲得較高的時間分辨率，必須減小時間窗的長度。然而，對于載頻差異較小且變化劇烈的時頻信號，STFT很難同時滿足時間和頻率分辨率的要求，而IAA譜估計算法具有高時頻分辨率和魯棒性，因而適合載頻差異較小的時頻分析問題。第3組仿真中，時頻信號包含2組載頻差異較小的線性和非線性調(diào)頻信號為

圖6 復雜信號時頻重構(gòu)圖Fig.6 Time-frequency reconstruction spectra of complicated signal

圖7 復雜信號分量個數(shù)估計Fig.7 Components estimation of complicated signal

仿真信號最小載頻差異約為0.008（歸一化），信號長度為100s，采樣頻率為512Hz，信噪比為0dB。圖8給出了信號進行STFT變換、Capon變換和IAA譜估計的時頻分布圖，3種算法均采用矩形窗，窗長N＝128，窗的移動步長為N。可以看到，在信號載頻差異較小的區(qū)域，STFT和Capon下信號分量均存在嚴重的重疊。

從圖9（a）～圖9（c）中可以看出，當信號分量間載頻差異較小時，ADPM-RBPF算法、基于Capon的峰值檢測法以及基于STFT的本文算法均無法分析出兩個載頻差異較小的信號分量。相反，如圖9（d）所示，基于IAA時頻譜的本文方法即使在載頻差異較小區(qū)域也可獲得較為理想的分析結(jié)果。圖10給出了每種算法都進行100次蒙特卡羅仿真的信號分量數(shù)目的估計結(jié)果。可以看到，只有基于IAA時頻譜的本文算法具有較為理想的估計效果。

圖8 載頻差異較小信號時頻圖Fig.8 Time-frequency spectra in close modes

圖9 載頻差異較小信號時頻重構(gòu)圖Fig.9Time-frequency reconstructionspectrainclosemodes

圖10 載頻差異較小信號分量個數(shù)估計Fig.10 Components estimation of signal in close modes

5 結(jié) 論

針對由多個信號分量構(gòu)成的時頻信號且其分量動態(tài)出現(xiàn)和消失的時頻譜分析問題，本文提出了一種分析、探測和跟蹤動態(tài)時頻譜的隨機有限集法。

1）該算法首先將觀測到的信號進行時頻變換，進而視每個時刻的時頻譜幅度的局部極值為測量。分析表明：每個時刻測量的集合是一個隨機有限集。

2）基于提出的狀態(tài)和測量隨機有限集模型，本文給出了利用高斯混合概率假設密度濾波器，并結(jié)合多項式預測模型以及提出的初始權(quán)重賦值算法，來進行動態(tài)時頻譜的分析、探測和跟蹤的算法。

3）仿真結(jié)果均證明：所提算法在分析精確度、魯棒性、對微弱時頻譜分量的探測能力、以及對載頻差異的分析能力等諸方面均優(yōu)于文獻中報道的算法。

式中：y＝為時間窗內(nèi)的測量數(shù)據(jù)，N為時間窗的寬度；I為離散化頻率 格 點 的 個 數(shù)；為離散化頻點τ對應的頻率矢量，tn為采樣時刻；eτ為離散化頻點τ對應的幅值；ε為均值為0且與信號相互獨立的高斯白噪聲。

對式（A1）所描述的信號模型，IAA算法的時頻譜估計結(jié)果為

式中：IAAy（n，τ）為n時刻IAA算法估計得到的譜；WJτ為第J次迭代中針對頻點τ的加權(quán)矩陣。IAA算法采用迭代技術，即每次迭代開始時都利用了上一次譜估計結(jié)果來構(gòu)建加權(quán)矩陣

式中：為第j次迭代中針對頻點τ的加權(quán)矩陣；為第j次迭代開始時針對頻點τ估計的信號 協(xié)方 差 矩 陣 ；為 第j－1次 迭 代 的 譜估計結(jié)果，其中初始譜估計結(jié)果由周期圖法獲得［23］。