(完整版)matlab信号的时频分析

信号的时频分析实验一、实验目的1、掌握Matlab对信号时频分析方法。2、掌握能量信号、周期性功率信号和非周期性功率信号的概念。3、掌握能量信号和功率信号的截断信号的时频域特性。4、掌握相关函数的概念及与功率谱的关系。二、实验原理1、能量信号的时频分析（1）能量信号：能量有限的信号，满足20()Estdt。如时间受限信号。（2）能量信号的频谱密度能量的频谱密度为：2()()jtSfstedt()Sf的逆变换为原信号：2()()jftstSfedf也可以表示为：()()jtSstedt，1()()2jtstSed（3）能量谱密度根据帕塞瓦定理（能量守恒），可以知道22|()||()|Estdtsfdf因此，可以将2()|()|Gfsf看成是信号的能量谱密度，表示能量随频率的分布。（4）能量信号的相关函数能量信号的自相关函数定义：()()()Rststdt能量信号的互相关函数定义：1212()()(),Rststdt（5）能量信号的相关函数和能量谱密度的关系222()222[()]()*()()*()[()]*()|()|ififtiftvtziftifvfRststeddtstesteddtstedtsvedvSf由此可看到，能量信号的自相关与其能谱密度是一对傅立叶变换对。2、功率信号的时频分析（1）功率信号：如果信号的能量无穷大，但其功率存在，则称该信号为功率信号，（2）功率信号的频谱函数设()st为周期性功率信号，0T为周期，则有000/220/201()()TjnftnTCCnfstedtT式中，0()Cnf为复数，表示为0()||njnCnfCe式中，||nC是频率为0nf的分量的振幅；n是频率为0nf的分量的相位。由傅立叶级数知识可知，上式为周期信号s(t)的傅里叶级数的系数，即周期信号s(t)可以表示为：0022/()jnftjntTnnnnstCeCe（3）功率谱密度令()st的截短信号为(),/2/2TstTtT，则()Tst为能量信号。有/222/2()()TTTTEstdtSfdf定义功率谱密度为：21()lim()TTPfSfT得到信号功率：/22/21lim()()TTTTPSfdfPfdfT若功率信号为周期信号，则：00/22220/201()|||()|()TnTnnPstdtCCffnfdfT，其中0,()0,nCfnfCff其他，由于()PPfdf，则20()|()|()nPfCffnf（4）功率信号的相关函数功率信号的自相关函数定义：/2/21()lim()()TTTRststdtT功率信号的互相关函数定义：/21212/21()lim()(),TTTRststdtT（5）功率信号的相关函数和功率谱密度的关系2()()jfPfRed2()()jfRPfedf相关函数和功率谱密度互为傅立叶变换。3、信号的频域分析方法（1）连续非周期信号的傅里叶变换（FT）对于连续非周期信号（能量信号），即满足绝对可积条件|()|stdt时。可以采用傅里叶变换来分析信号的频域特性。2()()jftSfstedt2()()jftstSfedf频域为连续非周期信号。（2）连续周期信号的傅立叶级数（FS）对于连续周期信号，如果满足Dirichlet条件，可以采用傅立叶级数来分析信号的频域特性。002001()()tTjnfttSnfstedtT0022/00()()()jnftjntTnnstSnfeSnfe当引入冲激函数时，可以用傅里叶变换来分析连续周期信号的频域。频域为离散非周期信号。（3）离散非周期信号的傅里叶变换（DTFT）2()()jfnnSfsne2()()jfnsnSfedf频域为连续周期信号。（4）离散周期信号的傅立叶级数（DFS）12/0()()NjnkNnSksne12/01()()NjnkNnsnSkeN频域为离散周期信号。时域的离散对应于频域的周期；时域的连续对应于频域的非周期；时域的周期对应于频域的离散；时域的非周期对应于频域的连续。（5）离散傅里叶变换（DFT）将离散周期信号的傅立叶级数各取一个周期（N个采样点），可以得到如下离散傅里叶变换。即将离散傅里叶变换作周期扩展，可以得到傅立叶级数。12/0()(),0,1,....,1NjnkNnSksnekN12/01()(),0,1,....,1NjnkNnsnSkenNN（6）快速傅里叶变换（FFT）对于离散傅里叶变换的规律，可以进行快速傅里叶变换FFT。它包含两类：一类为Cooley和Tukey提出的N等于2的整数次幂的算法；另一类为Winograd提出的N不等于2的整数次幂的算法。这两类方法都大大降低了DFT的乘法运算的次数。4、连续信号的离散化为了便于对信号通过计算机进行处理，如FFT。需要对信号进行采样。根据采样定理，当带限信号带宽为W时，采样频率2sFW，即采样间隔1/(2)sTW时。采样信号表示为：()()()ssnxtxnTtnT，其傅立叶变换为：1(/),()1(),||snssXfnTfTXfXffWT对全部，()Xf为()xt的傅立叶变换因此通过一个带宽为W，增益为sT的低通滤波器可以将原信号恢复。重建公式为：()()sin(2())ssnxtxnTcWtnT。由于离散序列的离散傅立叶变换(DFT)表示为2()()sjnTdnXfxne则()()/dsXfXfT。Matlab中只是对离散序列信号进行了DFT，因此，若要得到该离散信号的原连续信号的频谱，还需要将得到的信号乘以sT或除以sf，才能得到原连续信号的频谱。5、matlab对信号时频分析的方法（1）相关：Matlab中通过xcorr函数计算互相关。自相关：调用格式为c=xcorr(x)或c=xcorr(x,maxlags)，其中x为输入序列，c为它的自相关函数估计。maxlags为x的最大时延，若缺省，函数返回值c长度是2N-1；若不缺省，返回值c长度是2*maxlags+1。互相关：调用格式有四种c=xcorr(x,y)；c=xcorr(x,y,’option’)；c=xcorr(x,y,maxlags,’option’)；[c,lags]=xcorr(x,y,maxlags,’option’)。其中x,y为两个独立的信号序列，长度相同。maxlags意义同上。option为选择项：’bised’计算有偏互相关函数估计，'unbiased'计算无偏互相关函数估计，‘coeff’序列归一化使零延迟的自相关函数为1，‘none’为缺省情况。（2）傅立叶变换：Matlab中通过fft函数来计算离散傅立叶变换。（3）直接法求功率谱密度由公式21()lim()TTPfSfT可直接求得功率谱密度。（4）利用自相关函数求功率谱密度也可以通过自相关函数和功率谱密度之间的关系来得到信号的功率谱密度图。M=512;2()()jfPfRed2()()jfRPfedf自相关函数可由下面语句求得：M=N;Ryy=zeros(1,M);form=1:Mforn=1:N-m+1Ryy(m)=Ryy(m)+y(n)*y(n+m-1);endRyy(m)=Ryy(m)/(N-m+1);endsubplot(3,2,4)plot(t,Ryy(1:M));title('信号自相关');xlabel('时间(s)');axis([01-1.21.2]);三、实验步骤1、输入余弦信号()cos(2*80)stt，采样间隔取dt=0.001s，采样点数取N=512。将其变换到频域观察其频谱，再将其反变换到时域与原信号进行比较，并绘制该信号的时域和频域波形。2、通过直接求解法求该余弦信号的功率谱密度，并绘制该信号的功率谱密度波形。3、分别用xcorr函数和自编制自相关函数的方式，求该余弦信号的自相关函数，并绘制自相关函数波形。4、利用自相关函数和功率谱密度的关系且该余弦信号的功率谱密度，并绘制该信号的功率谱密度波形。5、改变抽样间隔dt=0.002s，重复步骤1、2、3、4。6、改变采样点数N=1024，重复1、2、3、4、5。7、输入信号()cos(2*80)2cos(2*50)sttt，重复1、2、3、4、5、6。四、实验结果1、写出完成实验步骤的程序。实验程序bt=0;%开始时间dt=0.001;%时间间隔N=1024;%傅立叶变换点数et=bt+N*dt-dt;%结束时间t=bt:dt:et;%时间域TT=et-bt;%总的时间y=cos(2*pi*80*t);%待分析的信号第二次为cos(2*pi*80*t)+2*cos(2*pi*50*t)figure(1);subplot(3,2,1)plot(t,y);title('信号');xlabel('时间(s)');axis([01-1.21.2]);Y=fftshift(fft(y,N));df=1/TT;%频率间隔Tf=N*df;%分析的频宽f=-Tf/2+df:df:Tf/2;%频率域subplot(3,2,2)plot(f,abs(Y));title('频谱');xlabel('频率(Hz)');axis([-1001000max(abs(Y))]);Pyy=abs(Y/Tf).^2/TT;%计算功率谱密度subplot(3,2,3)plot(f,Pyy);title('根据定义求出的功率谱');xlabel('频率(Hz)');axis([-1001000max(abs(Pyy))]);M=N;Ryy=zeros(1,M);form=1:Mforn=1:N-m+1Ryy(m)=Ryy(m)+y(n)*y(n+m-1);endRyy(m)=Ryy(m)/(N-m+1);endsubplot(3,2,4)plot(t,Ryy(1:M));title('信号自相关');xlabel('时间(s)');axis([01-1.21.2]);Pyy2=fftshift(fft(Ryy,M))/Tf;subplot(3,2,5)plot(f,abs(Pyy2));title('自编自相关函数求出的功率谱');xlabel('频率(Hz)');axis([-1001000max(abs(Pyy2))]);Ryy2=xcorr(y,y,'unbiased');L=length(Ryy2);Pyy3=fftshift(fft(Ryy2(L-N+1:L),N))/Tf;subplot(3,2,6)plot(f,abs(Pyy3));title('系统自相关函数求出的功率谱');xlabel('频率(Hz)');axis([-1001000max(abs(Pyy3))]);2、绘制实验步骤中要求的图形。1,n=512点dt=0.0012,n=512dt=0.0023,n=1024dt=0.0014,n=1024dt=0.0025,n=512dt=0.0016n=,512,dt=0,0027,n=1024dt=0.0018,n=1024dt=0.002分析：由以上实验可以验证采样点数越多，采样间隔越小越能放映真是的波形包络和频谱特征。