samp11_1_rankor_samp11_1_rankor_时频分析_源码

clear all; close all fm=50;%主频 dt=0.001;%时间域采样间隔 N=200;%采样点数 t=-N/2+1:N/2; x=(1-2*(pi*fm*t*0.001).^2).*exp(-(pi*fm*t*0.001).^2); figure(1);plot(t,x); title('Ricker-雷克子波道'); xlabel('时间t（ms）'); ylabel('幅值A'); for i=1:200 f(i)=10*(i-1);%计算频回率 end Y=abs(fft(x));%fourier变换，取振幅谱 figure(2) plot(f,Y); title('Ricker子波的振幅谱'); xlabel('频率f（hz）'); ylabel('振幅谱'); Nw=81;%窗口长度，一般应取大于1的奇数 nstep=1; %滑动步长，如果数据过长致使内存不够，可以使用较长的步长 h=window(@hamming,Nw); %给出窗函数,输出为列向量 Ts=[]; %存放时间中心的信息 L=floor(Nw/2); %窗口的半长度 F=[0:(Nw-1)/2]/(Nw*dt);%频率分量 TF=[]; %存放STFT的数组 for ii=1:nstep:N if(ii<L+1) %需要将前面补零以满足窗函数的要求 xw=[zeros(1,L-ii+1),x(1:ii+L)].*h'; %信号和窗函数在时间域内乘积等于在频率域内褶积 elseif(ii>N-L) %需要将后面补零以满足窗函数的要求 xw=[x(ii-L:N),zeros(1,(ii+L)-N)].*h'; %信号和窗函数在时间域内乘积等于在频率域内褶积 else xw=x(ii-L:ii+L).*h'; %信号和窗函数在时间域内乘积等于在频率域内褶积 end Ts=[Ts,ii]; %将时间序号存入时间矢量 temp=fft(xw,Nw); %对信号采用窗长度进行Fourier变换 TF=[TF,[temp(1:(Nw+1)/2)*2/Nw]']; %TF的一维为时间，一维为频率 end figure(3),mesh(Ts,F,abs(TF)); %绘出时频分布立体图 xlabel('时间/s'),ylabel('频率/Hz'),zlabel('振幅') %加上必要的标记 figure(4), pcolor(Ts,F,abs(TF)); %采用时频的平面分布 shading interp; %将图像进行平滑，如果没有此句，则时频平面不光滑 xlabel('时间/s'),ylabel('频率/Hz'); %加上必要的标记 colorbar %加上色标