短时DFT问题解惑

网上看到一段用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'; %信号和窗函数在时间域内乘积等于在频域内的卷积。

猜你喜欢