基于Radon-STFT变换的含噪LFM信号子空间分解

基于Radon-STFT变换的含噪LFM信号子空间分解

由于线性调频信号占有非常宽的频带，用奇异值分解就不能将含噪线性调频信号分解成信号子空间和噪声子空间.针对这一缺陷，本文提出了一种基于时频面旋转的含噪线性调频信号的子空间分解算法.文中分析了算法的性质，并提出了“伪信号子空间”的概念和用于检测直线倾角的Radon-STFT变换.理论分析和仿真实验的结果表明了这种子空间分解方法对一类含噪线性调频信号是可行的.
　　关键词:奇异值分解;伪信号子空间;时频分布;Radon-STFT变换;时频面旋转

Subspace Decomposition for Noisy LFM Signal Using Radon-STFT Transform

ZOU Hong-xing,ZHOU Xiao-bo,LI Yan-da
(Dept.of Automation,Tsinghua University,State Key Laboratory of Intelligent Technology and Systems,Beijing 100084,China)

　　Abstract:Since the LFM signal occupies a wide band in frequency domain,it's impossible to use singular value decomposition to separate the noisy LFM signal into signal subspace and noise subspace.To counter this drawback,a new subspace decomposition algorithm based on the rotation of time-frequency plane is presented in this paper along with its corresponding performance.A new concept,namely,the“pseudo signal subspace”,and a new transform for detecting the tilting angle called Radon-STFT transform are proposed.Theoretical predictions and simulation results indicate that the strategies advocated are feasible for denoising a class of LFM signals.
　　Key words:singular value decomposition；pseudo signal subspace；time-frequency distribution；Radon-STFT transform；rotation of time-frequency plane

一、引　　言
　　解线调法是一种专门估计线性调频(linear frequency modulation，简称LFM)信号参数（即起始频率f0和线性调频率r）的方法［1］，然而该方法失却了LFM信号的调幅(包络)信息.文献［2］提出了一种基于Radon-Wigner变换的LFM信号的滤波方法，但该方法存在两个明显的缺陷：(1)用于抽取单频正弦波的窄带滤波器的带宽对噪声抑制特性影响甚大，信号的带宽愈大，则重构出的信号中噪声含量愈显著；(2)窄带滤波器的设计需要在设计者的干预和指导下完成，算法无法自动实现这一设计过程.受文献［2］的启发，并为了克服其缺陷，作者曾提出一种增强LFM信号的时变滤波算法［3］，其核心思想是：首先，求出含噪LFM信号的瞬时频率，并由瞬时频率确定信号在时频平面上的倾斜角度；其次，根据这一倾斜角度，旋转含噪LFM信号的短时Fourier变换(STFT)，使得LFM信号成份的时频分布与时间轴平行；然后对旋转后的时频分布进行奇异值分解(SVD)，用低有效秩分离开信号和噪声子空间；最后反向旋转信号子空间，经短时Fourier逆变换，得出滤除掉噪声的LFM信号.这种确定旋转角度方法的不足之处在于瞬时频率的估计受噪声影响较大；当信噪比(SNR)过低时，很难正确地估计出瞬时频率.显然，若旋转角度估计不准确，则旋转后的时频平面的奇异值衰减就会变慢，用低秩逼近就会丢失信号的许多细节信息，导致重构出的信号产生较大的失真.此外，当量测信号中存在多个具有不同调频率的LFM信号时，就不可能通过求瞬时频率的途径拟合出旋转角度.本文提出一种Radon-STFT变换，该变换能鲁棒(robust)地确定LFM信号在时频平面上的倾斜角度，因而可以有效地克服文［3］的缺憾.

二、适用于含噪LFM信号的子空间分解算法
　　本文考虑具有如下形式的叠加波形

x(t)=s(t)+n(t)　(1)

其中n(t)是i.i.d的加性实值白噪声，s(t)是实值LFM信号

　(2)

上式中gμi，λi(t)是高斯包络.定义x(t)的STFT为

STFTx(t,f)=∫+∞-∞x(u)h*(t-u)e-j2πfudu　(3)

其中h为归一化窗函数，且满足‖h‖2=1.由于STFT是一种线性变换，显然有

STFTx(t,f)=STFTs(t,f)+STFTn(t,f)　(4)

成立.
　　1.Radon-STFT变换的定义
　　为准确检测LFM信号在时频平面上的倾斜角度，本文仿照Radon-Wigner变换，提出了Radon-STFT变换.Radon变换是一种积分投影变换，任意二维函数f(t,ω)的Radon变换为［1，2］

［f(t,ω)］=∫PQ线f(ucosα-vsinα,usinα+vcosα)dv　(5)

这里表示Radon变换，其积分路径如图1所示.

图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信号来说，有一个明显的时移.作者目前还不清楚造成这一现象的机理，这有待于进一步地研究.