網上看到一段用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'; %信號和窗函數在時間域內乘積等於在頻域內的卷積。函數