function[Spec,Freq]=STFT(Sig,nLevel,WinLen,SampFreq)%计算离散信号的短时傅里叶变换;%Sig待分析信号;%nLevel频率轴长度划分(默认值512);%WinLen汉宁窗长度(默认值64);%SampFreq信号的采样频率(默认值1);if(nargin1),error('Atleastoneparameterrequired!');end;Sig=real(Sig);SigLen=length(Sig);if(nargin4),SampFreq=1;endif(nargin3),WinLen=64;endif(nargin2),nLevel=513;endnLevel=ceil(nLevel/2)*2+1;WinLen=ceil(WinLen/2)*2+1;WinFun=exp(-6*linspace(-1,1,WinLen).^2);WinFun=WinFun/norm(WinFun);Lh=(WinLen-1)/2;Ln=(nLevel-1)/2;Spec=zeros(nLevel,SigLen);wait=waitbar(0,'Undercalculation,pleasewait...');foriLoop=1:SigLen,waitbar(iLoop/SigLen,wait);iLeft=min([iLoop-1,Lh,Ln]);iRight=min([SigLen-iLoop,Lh,Ln]);iIndex=-iLeft:iRight;iIndex1=iIndex+iLoop;iIndex2=iIndex+Lh+1;Index=iIndex+Ln+1;Spec(Index,iLoop)=Sig(iIndex1).*conj(WinFun(iIndex2));end;close(wait);Spec=fft(Spec);Spec=abs(Spec(1:(end-1)/2,:));Freq=linspace(0,0.5,(nLevel-1)/2)*SampFreq;t=(0:(SigLen-1))/SampFreq;clfset(gcf,'Position',[20100500430]);set(gcf,'Color','w');axes('Position',[0.10.450.530.5]);mesh(t,Freq,Spec);axis([min(t)max(t)0max(Freq)]);colorbarxlabel('t/s');ylabel('f/Hz');title('STFT时频谱图');axes('Position',[0.10.10.550.25]);plot(t,Sig);axistightylabel('x(t)');title('时域波形');axes('Position',[0.730.450.240.5]);PSP=abs(fft(Sig));Freq=linspace(0,1,SigLen)*SampFreq;plot(PSP(1:end/2),Freq(1:end/2));title('频谱');