基于Radon-STFT變換的含噪LFM信號(hào)子空間分解

圖1　Radon變換的積分路徑

　　若用信號(hào)x(t)的STFT取代函數(shù)f(t,ω)，則所得到的Radon變換就是信號(hào)x(t)的Radon-STFT變換(Radon-short-time-Fourier transform)，用符號(hào)RSTFTx(u,α)表示，即
　　定義1　信號(hào)x(t)的短時(shí)Fourier變換如式(3)所示.定義x(t)的Radon-STFT變換為

RSTFTx(u,α)=［STFTx(t,f)］=∫PQ線STFTx(ucosα-vsinα,usinα+vcosα)dv　(6)
=∫+∞-∞∫+∞-∞STFTx(u′cosα-v′sinα,u′sinα+v′cosα)δ(u-u′)du′dv′　(7)

當(dāng)u′=u時(shí)，Radon-STFT變換的積分路徑即為PQ線.
　　2.Radon-STFT變換的意義
　　為表述方便，采用截距f0和斜率r為參數(shù)表示STFT時(shí)頻平面中的直線.在求沿直線f=f0+rt的積分時(shí)，可將圖1中的積分路徑（直線PQ）參數(shù)(u,α)替換成(f0,r)，這兩對參數(shù)的關(guān)系為

　(8)

由式(7)求信號(hào)x(t)的Radon-STFT變換，并以參數(shù)(r,f0)表示其積分路徑，則有

　(9)

式(9)表明，若x(t)是參數(shù)為f0和r的LFM信號(hào)，則積分值最大；而當(dāng)參數(shù)偏離f0與/或r時(shí)，積分值迅速減小，即對一給定的LFM信號(hào)，其Radon-STFT變換會(huì)在相應(yīng)的參數(shù)(f0,r)處呈現(xiàn)尖峰.由于Radon變換是以旋轉(zhuǎn)角度α和半徑u為參數(shù)的，因此，通過檢測Radon-STFT變換的最大值位置，即可確定STFT時(shí)頻平面中的直線傾角.
　　誠然，除了用Radon-STFT變換檢測時(shí)頻平面中直線的傾角外，還可采用Radon-Wigner變換［1，2］(或Radon變換與其它Cohen類時(shí)頻分布相結(jié)合的變換).用這兩種變換檢測出的傾角并不一致（由STFT的冗余性引起），但滿足簡單的幾何關(guān)系.雖然Wigner分布的時(shí)頻分辨率達(dá)到了時(shí)頻分辨率理論值的下限，而這一優(yōu)良特性是STFT所無法企及的，但這種分辨率的差異幾乎不影響Radon變換對時(shí)頻平面中直線傾角的檢測.究其原因，乃是Radon變換本質(zhì)上是沿二維平面中不同直線方向上的線積分，LFM信號(hào)的STFT沿分布的長對稱軸方向上的幅值，比沿同一方向上的其它軸上的幅值要大，這就保證了Radon變換能正確地檢測出LFM信號(hào)在時(shí)頻平面上的傾角.當(dāng)信號(hào)受i.i.d的白噪聲污染時(shí)，只要SNR不是太低（譬如－30dB以下，但鮮見于實(shí)際問題中），仍能保證檢測結(jié)果的準(zhǔn)確性.這是由于i.i.d的白噪聲趨于均勻地布滿整個(gè)時(shí)頻平面，而信號(hào)的STFT則比較集中，信號(hào)的分布顯然是時(shí)頻平面中的占優(yōu)分布.作者在一個(gè)很寬的SNR范圍內(nèi)（10～-10dB）進(jìn)行了3×100×10＝3000次Monte-Carlo試驗(yàn)（對3種不同線性調(diào)頻率的信號(hào)在10～-10dB之間共100個(gè)SNR環(huán)境下進(jìn)行試驗(yàn)，在相同的SNR下做10次），結(jié)果全部正確.此即作者考慮構(gòu)造Radon-STFT變換的一個(gè)主要?jiǎng)右?就本文的子空間分解算法而言，Radon-STFT變換較之Radon-Wigner變換還有另外兩個(gè)優(yōu)勢，一則省掉了耗時(shí)的Wigner分布的計(jì)算，二則所采用的STFT時(shí)頻平面僅為Wigner分布時(shí)頻平面的一半，這對于計(jì)算量很大的Radon變換來說，無疑減輕了不少計(jì)算負(fù)擔(dān).
　　3.信號(hào)子空間的低秩逼近
　　令A(yù)x=STFTx(t,f)，As=STFTs(t,f)，An=STFTn(t,f)，則有Ax=As+An.Ax的奇異值分解表示為［4］

　(10)

這里H表示共軛轉(zhuǎn)置，Σx=diag(σx1，σx2，…，σxr)，相應(yīng)的奇異值滿足σx1≥σx2≥…≥σxr≥0，r=rank(Ax).信號(hào)／噪聲子空間分解，就是用Ax的低秩矩陣Axk去逼近As,

Axk=UΣxkVH　(11)

其中Σxk是通過在Σx內(nèi)令除去k個(gè)最大的奇異值以外的所有其它奇異值都等于零后得到的對角陣.記Σs=diag(σs1，σs2，…，σsr)為信號(hào)子空間的奇異值，Σn=diag(σn1，σn2，…，σnr)為噪聲子空間的奇異值，并且有σxi=σsi+σni，1≤i≤r.若噪聲是i.i.d的，則σn1=σn2=…=σnr.為了有效地恢復(fù)信號(hào)的時(shí)頻分布，就需要從σxi中去除σni，因此如何估計(jì)σni，或等效地，如何估計(jì)Σx的有效秩，即成為問題的關(guān)鍵.不妨設(shè)估計(jì)出的有效秩為k，則σxk為一閾值，前面的k個(gè)奇異值對應(yīng)于信號(hào)子空間，后面的r-k個(gè)奇異值對應(yīng)于噪聲子空間，因此問題歸結(jié)為如何有效地估計(jì)閾值σxk，以保證重構(gòu)信號(hào)的精度.
　　為解決此問題，我們需要兩個(gè)定理.為簡便計(jì)，假定需研究的模型為

wi=θi+δzi　(12)

其中δ＞0為噪聲方差，zi服從標(biāo)準(zhǔn)正態(tài)分布，即zi～N(0,1).
　　定義軟閾值變換與硬閾值變換，其中閾值τ＞0.
　　定義2　軟閾值變換定義為

　(13)

　　定義3　硬閾值變換定義為

(14)

前文述及的低有效秩重構(gòu)信號(hào)的出發(fā)點(diǎn)顯然是硬閾值變換，因此本文僅討論在硬閾值變換下的信號(hào)重構(gòu)問題.有關(guān)軟閾值變換的理論分析，本文不擬介紹.
　　定義4　估計(jì)量的測度或稱損失函數(shù)定義為

　(15)

其中M為采樣點(diǎn)數(shù)，‖.‖22,M為L2范數(shù).
　　根據(jù)Bickel［5］和Donoho［6］的工作，可知如下定理成立.
　　定理1　設(shè)lM為一逼近到的閾值序列，則硬閾值估計(jì)為i=wiχ|wi|＞lMδ，其中χ為特征函數(shù)，lM～2logM，并且

　(16)

這里θ=(θ1，…，θM)，為θ的估計(jì).
　　設(shè)(t)為x(t)的估計(jì).對x(t)的短時(shí)Fourier變換STFTx(t,f)作硬閾值變換，然后通過短時(shí)Fourier逆變換得重構(gòu)信號(hào)(t)，因之有如下定理.
　　定理2

　(17)

　　證明　設(shè)為對STFTx(t,f)所作的硬閾值變換結(jié)果，則

=ISTFT。。STFTx(t,f)　(18)

其中“。”為映設(shè)或操作算子，ISTFT為短時(shí)Fourier逆變換.
　　根據(jù)Parseval等式，有

E‖-x‖22,M=E‖-θ‖22,M　(19)

再根據(jù)定理1即得式(17)成立.定理2證畢.
　　定理2表明，對STFTx(t,f)作硬閾值變換(當(dāng)然這種硬閾值變換是通過正交分解進(jìn)行的)，然后通過短時(shí)Fourier逆變換重構(gòu)出的信號(hào)的均方差與同階，從而保證了重構(gòu)信號(hào)的精度.
　　以上定理成立的基礎(chǔ)是閾值δ，而這個(gè)閾值是由極小極大準(zhǔn)則［5］得到的.實(shí)際應(yīng)用中閾值的取法有多種；在不同的準(zhǔn)則下會(huì)得到不同的閾值.后文介紹的秩1逼近，就是一種硬閾值變換.
　　4.算法步驟
　　根據(jù)前文的分析，本節(jié)給出基于Radon-STFT變換的含噪LFM信號(hào)的子空間分解算法.
　　步驟1　計(jì)算實(shí)值含噪信號(hào)x(t)的解析信號(hào)z(t)

z(t)=x(t)+［x(t)］　(20)

其中表示Hilbert變換

　(21)

上式中p.v.表示積分的Cauchy主值［7］；
　　步驟2　用式(3)計(jì)算z(t)的短時(shí)Fourier變換STFTz(t,f)；
　　步驟3　令Bz=STFTz(t,f/2)，即存儲(chǔ)一半的STFT數(shù)據(jù)，而對應(yīng)于負(fù)頻率的另一半STFT數(shù)據(jù)則由于采用了Hilbert變換而變?yōu)榱阒?用式(7)計(jì)算Radon-STFT變換，找出與最大值點(diǎn)相對應(yīng)的角度α.通過簡單的三角函數(shù)關(guān)系可知，直線的傾角β=α-90°.參數(shù)β用于控制時(shí)頻平面的旋轉(zhuǎn)角度；
　　步驟4　將時(shí)頻平面STFTz(t,f)旋轉(zhuǎn)β°，使得信號(hào)的時(shí)頻分布平行于時(shí)間軸（這一操作稱為調(diào)平），并將旋轉(zhuǎn)后的時(shí)頻平面記為z.這一旋轉(zhuǎn)過程可表示為

z=Tβ［STFTz(t,f)］　(22)

其中Tβ表示旋轉(zhuǎn)β°；
　　步驟5　用式(10)計(jì)算z的SVD，即
　　步驟6　偽子空間分解.令有效秩k=1，z2=z3=…=zM=0，用式(11)重構(gòu)出偽信號(hào)子空間

　(23)

由于z1并非真正的信號(hào)子空間，作者姑且將其稱為“偽信號(hào)子空間”，相應(yīng)的子空間分解也稱為“偽子空間分解”.采用有效秩1逼近z是基于以下的事實(shí)：在旋轉(zhuǎn)后的時(shí)頻平面z上，被旋轉(zhuǎn)成單音信號(hào)的分布相對于背景噪聲來說仍是占優(yōu)分布，其能量非常集中，而白噪聲的能量則趨于均勻彌散于原先的半個(gè)時(shí)頻平面.反映到SVD的奇異值上，就是奇異值衰減訊速，第一個(gè)奇異值和第二個(gè)奇異值之間有較大的躍變，因而信號(hào)的能量基本上反映到第一個(gè)奇異值上，用秩1就足夠重構(gòu)出偽信號(hào)子空間——定理2保證了z1逼近z的精度；
　　步驟7　將z1再反向旋轉(zhuǎn)β°，得到傾斜的時(shí)頻分布，并將反向旋轉(zhuǎn)后的時(shí)頻分布存貯為STFT′z(t,f).這一過程可表示為

STFT′z(t,f)=T-β［z1］　(24)

其中T-β表示反向旋轉(zhuǎn)β°；
　　步驟8　時(shí)頻面截取.由于旋轉(zhuǎn)運(yùn)算會(huì)使時(shí)頻平面增大為原時(shí)頻平面的(cosβ+sinβ)倍，兩次旋轉(zhuǎn)運(yùn)算則增大為(cosβ+sinβ)2倍，因此需要將旋轉(zhuǎn)變換后的時(shí)頻平面STFT′z(t,f)截取成與原時(shí)頻平面STFTz(t,f)相同大小的時(shí)頻平面.設(shè)STFTz(t,f)為M×M的矩陣，STFT′z(t,f)為l×l的矩陣，則截取從第(l-M)/2行到(l+M)/2行和第(l-M)/2列到(l+M)/2列的STFT′z(t,f)，就為所需要的信號(hào)子空間，并將其記為STFT(t,f).事實(shí)上，被截掉的時(shí)頻平面均為零值.
　　得到信號(hào)子空間之后，剩下的問題是如何由STFT(t,f)重構(gòu)出LFM信號(hào)s(t)的估計(jì)(t).如果STFT(t,f)確實(shí)是某一物理信號(hào)的STFT，則可用濾波器組求和或重疊相加法［8］精確地重構(gòu)出信號(hào)(t)，否則這兩種方法的應(yīng)用就會(huì)受到限制.利用最小均方誤差準(zhǔn)則，得［8］

　(25)

與基于Cohen類雙線性時(shí)頻分布的信號(hào)重構(gòu)方法相比，式(25)所示的信號(hào)重構(gòu)方法省去了相位優(yōu)化過程.這是因?yàn)镾TFT保留了原信號(hào)的相位信息，而Cohen類的雙線性時(shí)頻分布則舍棄了相位信息.這也是作者考慮構(gòu)造Radon-STFT變換的另一個(gè)原因.

三、仿真實(shí)例
　　實(shí)驗(yàn)信號(hào)由式（1）和式（2）產(chǎn)生，參見圖2(f)的上部，信號(hào)的長度為256點(diǎn)，LFM信號(hào)的起始頻率為0.05Hz，終止頻率為0.3Hz，SNR為3dB.圖2示出了本文提出的子空間分解算法的計(jì)算過程及相應(yīng)的時(shí)變?yōu)V波結(jié)果.不失一般性，信號(hào)的采樣頻率均歸一化成1Hz.圖2(a)是含噪LFM信號(hào)的STFT；圖2(b)是其Radon-STFT變換，從圖中可以明顯看出，LFM信號(hào)在104°(對應(yīng)于直線的傾角β＝104°－90°＝14°)處有一尖峰；圖2(c)為調(diào)平后的時(shí)頻平面；圖2(d)為偽信號(hào)子空間，從該圖中可以發(fā)現(xiàn)，噪聲背景已基本濾除；圖2(e)顯示了經(jīng)反向旋轉(zhuǎn)和截取后的時(shí)頻平面，此即為所需要的信號(hào)子空間；圖2(f)的下部示出了用信號(hào)子空間重構(gòu)的LFM信號(hào)，即時(shí)變?yōu)V波結(jié)果.作為對比，作者還設(shè)計(jì)了4階最小均方誤差自適應(yīng)FIR濾波器(其中學(xué)習(xí)步長定為相關(guān)噪聲自相關(guān)矩陣的最大特征值的倒數(shù)的0.1倍，遺忘因子置為0.0)，相應(yīng)的濾波結(jié)果在圖2(f)的中部給出.對于自適應(yīng)濾波來說，學(xué)習(xí)步長對其影響很大，要獲得滿意的濾波效果，需采用較小的學(xué)習(xí)步長，但學(xué)習(xí)時(shí)間明顯增加.從圖中可以看出，用本文提出的算法所得到的濾波結(jié)果，具有良好的保幅特性.注意到圖中自適應(yīng)濾波結(jié)果的起始部分有較大的誤差.

圖2　含噪LFM信號(hào)子空間分解算法的計(jì)算過程及相應(yīng)的濾波結(jié)果

四、結(jié)束語
　　LFM信號(hào)是一大類人工信號(hào)（譬如雷達(dá)信號(hào)）的基礎(chǔ)，對淹沒于噪聲中的LFM信號(hào)作濾波處理，具有廣泛的現(xiàn)實(shí)意義.本文給出了應(yīng)用于單個(gè)含噪LFM信號(hào)的子空間分解算法，不難推測，該算法很容易被推廣到多個(gè)含噪LFM信號(hào)的情形.需要指出的是，對于多個(gè)LFM信號(hào)來說，相應(yīng)于每一個(gè)LFM信號(hào)，Radon-STFT變換均有一峰值與之對應(yīng).然而要自動(dòng)地確定峰值的數(shù)目，實(shí)非易事（除非給出一閾值），此時(shí)設(shè)計(jì)者的干預(yù)仍是不可或缺的.幸運(yùn)的是，相對于噪聲背景來說，LFM信號(hào)在時(shí)頻平面上非常集中，很容易從時(shí)頻平面上確定LFM信號(hào)的數(shù)目，因而也就有了Radon-STFT變換的峰值數(shù)目.另外，一個(gè)不容忽視的事實(shí)是，圖2(f)中的時(shí)變?yōu)V波結(jié)果相對于不含噪聲的LFM信號(hào)來說，有一個(gè)明顯的時(shí)移.作者目前還不清楚造成這一現(xiàn)象的機(jī)理，這有待于進(jìn)一步地研究.