基于matlab的信号时频分析仿真

文章编号:1002-6886(2005)05-0001-03

・1・

基于matlab的信号时频分析仿真

戴幻尧

(空军航空大学,吉林　长春　130022)

摘要:MATLAB软件在多个研究领域都有着广泛的应用。其中,它的频谱分析和滤波器的分析设计功能很强,从而使数字信号处理变得十分简单、直观。本文介绍了时频分析基础理论及一些应用,运用MATLAB语言实现了旨在构造一种时间和频率的密度函数,以揭示信号中所包含的频率分量及其演化特性的wigner2ville分布。关键词:时频分析　wigner2ville分布　matlab　解析信号

ApplicationandStudyofMatlabonTime2frequencyAnalysisofSignal

DAIHuan2yao

Abstract:ThesoftwareofMATLABhasgotextensiveapplicationinseveralresearchesrealm.Amongthem,itsfrequencychartanalysisanddesignfunctionofthefilterwiththeanalysisisverystrong,makearithmeticdigitalsignalhandledtobecomeverybrief,intuitionistic.Thistextintroducedtime2frequencyanalysisfoundationtheoriesandsomeapplied,makeuseoftheMATLABlanguagetoconstructakindoftimeandfrequencydensityfunction,aimsatexplainingsignalcontainingfrequencyweightanditsevolvecharacteristicbywigner2villedistribute.Keywords:time2frequencyanalysis;wigner2villedistribute;matlab;analyticssignal

1　引言

时频分析能清楚地揭示信号的时变频谱特性,是对时变、非平稳信号进行分析与处理的有力工具。分离谱处理在信号展开过程中使用频率窗,这种给信号加窗的方法正是时频变换的基本手段,因此可以利用时频分析的思想,通过研究信号在相空间上的特性,了解其时频局部化信息及能量分布,从而为确定信号频带宽度提供依据。本文在时频分析的基础上,提出了确定超声信号的频带范围及频率窗个数的方法。

WD具有很高的时间和频率分辨率,同时又有其它一些优良的性质。WD对单个chirp信号其结果是非常理想的,但对于两个信号的和,由于存在交叉项,严重地影响了信号的时频表示。关于交叉项的消除(或减弱)方法,可以先将实信号转变成解析信号,再进行WD分析,可消除这种频谱正负部分之间的交叉干扰项。文章给出了实验加以论证。而随着MATLAB软件尤其是MATLAB的信号处理工具箱Simu2link仿真工具的不断完善,不仅使时频分析设计有了可能,而且还可以使设计达到最优化。

域信号映射到一个二维的时域平面,全面反映观测信号的时频联合特征。短时傅立叶变换反映了这一思想,对于时变信号,采用某一滑动窗函数截取信号,并认为这些信号是准平稳的,然后,再分别对其进行傅立叶变换,构成时变信号的时变谱。短时傅立叶变换是一种常用的时—频域分析方法,其基本思想是在傅立叶变换的基础上实现时域的局部化。由于傅立叶变换在时域和频域的对偶关系,所以短时傅立叶变换可以从时域、频域来描述,相应的短时傅立叶变换可以从傅立叶变换以及频域滤波的观点来考虑。

短时傅立叶变换的定义(1):

∞

STFTx(n,w)=

m=-∞

∑x(m)w(n-

w)e-

jwm

(1)

式中:w(n)是一个窗函数,其作用是取出;在x(n)某时刻附近的一小段信号进行傅立叶变换,当n变化时,窗函数随n移动,从而得到信号频谱随时间n变化的规律,此时的傅立叶变换是一个二维域(n,w)的函数。

2.2维格纳变换:WD(WignerDistribution)

信号f(t)的WD定义为:

Wf(t,w)=

2　时频分析基础

2.1短时傅立叶变换STFT[8]

从历史上看,信号的时频分析用的最多的是短时傅立叶变换,这种变换的基本思想是用一个窗函数乘时间信号,该窗函数的时宽足够窄,使取出的信号可以被看成是平稳的,然后进行的傅立叶变换可以反映该时宽中的频谱变化规律,如果让窗函数随时间轴移动,可以得到信号频谱随时间变化的规律。

对于时变信号,了解不同时刻附近的频域特征是至关重要的。因此,人们采用时间—频率描述时变信号,将一维的时

τ)ff(t+∫2

-∞

3

(t-

τ-)e

2

jwt

τ(时间域)d(2)

θθθ1)ejtdθ(频率域)(3)X(Ω+)X3(Ω-π2m=+∞22

其中,f(t)是实际信号s(t)的解析信号,积分限从-∞到+∞。WVD变换在分析信号时具有很好的能量集中特性,但存在的问题是由于它是一种非线性变换,因此它在分析多分量信号时将产生交叉干扰项(Cross-terms)和人工产物(Artifacts)。WD存在两个主要问题:一是有负值,这就影响了WD作为信号能量密度的物理解释;二是对两个信号的

WDx(Ω,t)=

∫

　作者简介:戴幻尧(1982—),长春市人,2004年于空军航空大学获得学士学位,现在空军航空大学攻读电子对抗硕士学位,研究方向是雷达侦

察和信号处理。

　收稿日期:2005-4-12

・2・

和,其WD存在交叉项,交叉项出现在自主项中间,且呈振荡状态。由于交叉项的存在,势必会影响到时频表示的分辨率。所以人们提出了许多降低干扰的核设计方法和分布。

现代机械　2005年第5期　　

contour(F,T,abs(TF));　　对比实信号x(t)=2cos(jπkt2),k=4,0≤T≤5,信号带宽:fc=kT=30,采样频率fS

=4fc采样点数N=TfS=400,信号加了高斯白噪声,用短时

3　STFT和WD变换程序设计实例3.1STFT的MATLAB设计

jkt

设计过程:考虑如下信号:解析信号x(t)=eπ,技术指标为:k=6,0≤T≤5,信号带宽:fc=kT=30,采样频率

fS=4fc,采样点数N=TfS=600,用短时傅立叶变换分析其在时域、频域特性,经过STFT和WD运算后的时频分布结果特点,从而将他们识别。

设计的源程序如下:k=4;T=5;fc=k3T;fs=43fc;Ts=1/fs;N=T/Ts;

x=zeros(1,N);t=0:N-1;

x=exp(j3k3pi3(t3Ts).^2);subplot(2,2,1);

plot(t3Ts,real(x));X=fft(x);

X=fftshift(X);subplot(2,2,2);

plot((t-N/2)3fs/N,abs(X));Nw=20;L=Nw/2;

Tn=(N-Nw)/L+1;nfft=32;

TF=zeros(Tn,nfft);fori=1:Tn

　　xw=x((i-1)310+1:i310+10);　　temp=fft(xw,nfft);　　temp=fftshift(temp);　　TF(i,:)=temp;end

subplot(2,2,3);

fnew=((1:nfft)-nfft/2)3fs/nfft;tnew=(1:Tn)3L3Ts;

[F,T]=meshgrid(fnew,tnew);mesh(F,T,abs(TF));subplot(2,2,4);

傅立叶变换分析其特性。

2

　　可以看出,由于实信号的频谱包括信号中的每个物理频率对应的正负频率部分,致使其所含信号分量为有物理意义的信号分量的两倍。如果直接采用实信号进行时频分析,我们将得到包括所有正负频率分量信号项的干扰。当有噪声存在时,时频谱变差。在时间相干性不强的目标及海量信息将被消弱,当信号环境更加复杂时,干扰较强的情况下,真个频带内部都会有射频干扰形成的谱峰存在,难于设置相应的频域滤波器,而且也会造成其他有用信息过多的损失。

3.2WD变换的MATLAB设计

πkt2),用wigner-ville分布分析解析信号x(t)=exp(j

技术参数:k=6,0≤T≤4;信号带宽:fc=kT=24;采样频率:fS=4fc;采样点数N=TfS=384,信号加了复高斯白噪声,用WD仿真并分析其特性。

　　从输出的时频平面中可以清晰的看出频率随时间线形变化,其变化率为k=6,可见整个仿真分析是正确的。经过仿真对比实验,对于带宽为fc的模拟信号进行FFT分析,采样率至少是截止频率的两倍,即fS=2fc,这就是奈奎斯特准则的要求。然而,对于WD变换,采样频率至少是信号截止频率的四倍,即fS=4fc,才能保证不会混叠,这一点是在多次

　　设计・研究・分析

实验中发现的规律,很重要。这次仿真也同时说明了时频变换后的信噪比比较高,WD变换对噪声不太敏感。

・3・

　　即强弱两信号组成的混合型号的瞬时频率近似为强信号的瞬时频率。

4　应用强线形调频干扰下雷达目标的识别[7]

近年来得到广泛研究的时—频分布是一种分析时变、非平稳信号的有力工具。因此,用时—频分布的方法对目标回波在多普勒域进行分析,有可能达到目标识别和分离的目的。时—频分布中最为典型的两种方法是Wigner分布和短时付氏变换(STFT)。Wigner分布的时—频分辨率最高,且能满足边缘分布、有限时频支撑等许多良好的特性,其最大的缺陷是在多分量信号之间存在着严重的交叉项,而这些交叉项又不代表任何真实的信号。我们称这些交叉项为信号的互分量,这些互分量与表明真实信号的自分量混合在一起,使我们无法确定真实信号。Wigner分布的另一缺点是当信噪比较低时,其分布受噪声的影响很严重,以致于无法区别信号。STFT与Wigner分布相反,它的时、域分辨力受窗长的约束,但它不存在交叉项。在时—频平面表示信号时,一般都采用谱图(spectrogram,即STFT的模平方)。如何减小Wigner分布的交叉项和提高STFT的分辨力是当前时—频分析方法研究的一个热门课题。目前尚没有一种万能的方法。公认较好的解决方法都是自适应的,即根据信号的特点选择相应的核函数来对信号的模糊函数进行加权,进而得到较好的时-频分布。实际上,我们可以将Wigner分布看成是STFT的一种特例,即采用信号自身做为窗函数的STFT。在对信号自分量的分辨方面,由于用信号自身为窗,对于LFM信号使得加窗后的信号为单谐波函数,且窗长由信号的持续时间自动地决定,即充分利用了信号的长度。因此,加窗后的信号既是单频函数又在时间上最长,采用付氏变换就可以得到最高的频率分辨。

假设两信号x1(t),x1(t):频率分别为f1(t)、f2(t),相位

θ分别为θA2,且A2≥A1,组成的混1(t)、2(t),幅度分别为A1、

合信号为x(t):

πθx(t)=x1(t)+x2(t)=A1exp[j21(t)]+←

πθτ)][5]→A2exp[j22(θπθ=A1cos[2π1(t)]+jA1sin[21(t)]+←θπθ(4)→A2cos[2π2(t)]+jA2sin[21(t)]πθπθ={A1cos[21(t)]+A2cos[22(t)]}+←πθπθ→j{A1sin[21(t)]+A2sin[22(t)]}则混合信号x(t)的瞬时频率f(t):

πθπθA1cos[2d1(t)]+A2cos[22(t)]f(t)=tan-1

πθπθdtA1sin[21(t)]+A2sin[22(t)]

=dtan-1dt

A1πθπθcos[21(t)]+cos[22(t)]

A2A1πθπθsin[21(t)]+sin[22(t)]

A2

A1≈0,则:A2

(6)(5)

因为设置信号强度大致相当,所以从信号的频谱中可以观测到信号中存在着两种频率,但不能得出它们的时间持续范围,而从信号的时频平面中不仅可以看出这两种信号的频率,而且还能得出其时间持续范围。4　结论本文的研究表明,利用短时傅立叶变换能够处理非平稳随机信号,但是在窗函数下的一个短的时间间隔是伪平稳的,在时-频窗的形状固定不变时,窗口面积越小,它的时-频局部化描述能力就越强;窗口面积越大,描述能力就越差,并且容易受到各种噪声的影响。Wigner分布是对瞬时自相关运算的傅氏变换,通过对瞬时特性的谱分析,获得时变信号的瞬时频率谱,同时也有利于抑制各种非相关的噪声。同瞬时自相关算法一样,它也是一种非线性运算,不适用于多信号混叠的场合。

参考文献

1李敬勇.对线性调频脉冲压缩雷达干扰的时频分析.1998年,第13

卷,第三期

2纪跃波,秦树人,汤宝平.Wigner分布干扰项抑制及其算法.重庆大

学学报,2001年第4期

3孙晓兵,保铮,罗琳.时频信号分析和雷达的多目标分辨.系统工程

与电子技术,1997年(1)

4张善文,甄蜀春,赵兴录.编队目标架次的识别.航天电子对抗,2001(1)

5张善文,甄蜀春,赵兴录.强线性调频干扰下雷达目标的识别方法.

航天电子对抗2000(4)

6胡国胜,张国红.一种新的信号处理方法———线性调频小波变换.数

学的实践与认识,2003年2月

7孙孔峰,皇甫堪.国防科技大学学报.1997年6月Vol.19No.38丁玉美,阔永红,高新波.数字信号处理.20029张贤达.现代信号处理.清华大学出版社,199410侯慧群.雷达对抗原理.空军第二航空学院,2001

因为A2≥A1,即

πθcos[2d2(t)]f(t)≈tan-1≈f2(t)

πθdtsin[21(t)]