圖1　Radon變換的積分路徑

　　若用信號x(t)的STFT取代函數(shù)f(t,ω)，則所得到的Radon變換就是信號x(t)的Radon-STFT變換(Radon-short-time-Fourier transform)，用符號RSTFTx(u,α)表示，即
　　定義1　信號x(t)的短時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)

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

　(8)

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

　(9)

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

Axk=UΣxkVH　(11)

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

wi=θi+δzi　(12)

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

　(13)

　　定義3　硬閾值變換定義為

(14)

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

　(15)

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

　(16)

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

　(17)

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

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

其中“?！睘橛吃O(shè)或操作算子，ISTFT為短時Fourier逆變換.
　　根據(jù)Parseval等式，有

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

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

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

其中表示Hilbert變換

　(21)

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

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

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

　(23)

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

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

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

　(25)

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

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

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

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