圖1　Radon變換的積分路徑

　　若用信號x(t)的STFT取代函數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為參數表示STFT時頻平面中的直線.在求沿直線f=f0+rt的積分時，可將圖1中的積分路徑（直線PQ）參數(u,α)替換成(f0,r)，這兩對參數的關系為

　(8)

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

　(9)

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

　(10)

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

Axk=UΣxkVH　(11)

其中Σxk是通過在Σx內令除去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的有效秩，即成為問題的關鍵.不妨設估計出的有效秩為k，則σxk為一閾值，前面的k個奇異值對應于信號子空間，后面的r-k個奇異值對應于噪聲子空間，因此問題歸結為如何有效地估計閾值σxk，以保證重構信號的精度.
　　為解決此問題，我們需要兩個定理.為簡便計，假定需研究的模型為

wi=θi+δzi　(12)

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

　(13)

　　定義3　硬閾值變換定義為

(14)

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

　(15)

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

　(16)

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

　(17)

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

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

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

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

再根據定理1即得式(17)成立.定理2證畢.
　　定理2表明，對STFTx(t,f)作硬閾值變換(當然這種硬閾值變換是通過正交分解進行的)，然后通過短時Fourier逆變換重構出的信號的均方差與同階，從而保證了重構信號的精度.
　　以上定理成立的基礎是閾值δ，而這個閾值是由極小極大準則［5］得到的.實際應用中閾值的取法有多種；在不同的準則下會得到不同的閾值.后文介紹的秩1逼近，就是一種硬閾值變換.
　　4.算法步驟
　　根據前文的分析，本節給出基于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數據，而對應于負頻率的另一半STFT數據則由于采用了Hilbert變換而變為零值.用式(7)計算Radon-STFT變換，找出與最大值點相對應的角度α.通過簡單的三角函數關系可知，直線的傾角β=α-90°.參數β用于控制時頻平面的旋轉角度；
　　步驟4　將時頻平面STFTz(t,f)旋轉β°，使得信號的時頻分布平行于時間軸（這一操作稱為調平），并將旋轉后的時頻平面記為z.這一旋轉過程可表示為

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

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

　(23)

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

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

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

　(25)

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

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

圖2　含噪LFM信號子空間分解算法的計算過程及相應的濾波結果

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