网上看到一段用matlab写的代码对短时DFT来处理信号:数组
%某一数字信号序列x(n),n=0,1,2,...,N-1 , N=500, 在(50,150)和在(250,300)范围内分别有两个不一样的正弦信号。 %采用STFT分析其时频分布,采用长度为81的hanmming窗,时间的滑动步长为1 clear all; close all N=500; %信号的长度 dt=1; %信号的采样间隔 t=0:dt:N-1; %时间序列 x=zeros(size(t)); x(50:150)=cos(2*pi*1/20*(t(50:150)-50)); %频率为0.05Hz x(250:350)=cos(2*pi*1/40*(t(250:350)-250)); %频率为0.025Hz subplot(2,2,1),plot(x),xlabel('时间/s');grid on; X=fft(x); subplot(2,2,2),plot([0:N-1/2]/N/dt,abs(X*2/N));grid on;%给出频率轴的表示 xlabel('频率/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 subplot(2,2,3),mesh(Ts,F,abs(TF)); %绘出时频分布立体图 xlabel('时间/s'),ylabel('频率/Hz'),zlabel('振幅') %加上必要的标记 subplot(2,2,4), pcolor(Ts,F,abs(TF)); %采用时频的平面分布 shading interp; %将图像进行平滑,若是没有此句,则时频平面不光滑 xlabel('时间/s'),ylabel('频率/Hz'); %加上必要的标记 colorbar %加上色标
其中有一句:xw=x(ii-L:ii+L).*h'; %信号和窗函数在时间域内乘积等于在频域内的卷积。函数