function Spec,Freq=STFT(Sig,nLevel,WinLen,SampFreq)%计算离散信号的短时傅里叶变换；% Sig 待分析信号；% nLevel 频率轴长度划分（默认值512）；% WinLen 汉宁窗长度（默认值 64）；% SampFreq 信号的采样频率 （默认值1）；if (nargin 1), error(At least one parameter required!);end;Sig=real(Sig);SigLen=length(Sig);if (nargin 4), SampFreq=1;endif (nargin 3), WinLen=64;endif (nargin 2), 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,Under calculation,please wait.);for iLoop=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,20 100 500 430);set(gcf,Color,w);axes(Position,0.1 0.45 0.53 0.5);mesh(t,Freq,Spec);axis(min(t) max(t) 0 max(Freq);colorbarxlabel(t/s);ylabel(f/Hz);title(STFT时频谱图);axes(Position,0.1 0.1 0.55 0.25);plot(t,Sig);axis tightylabel(x(t);title(时域波形);axes(Position,0.73 0.45 0.24 0.5);PSP=abs(fft(Sig);Freq=linspace(0,1,SigLen)*SampFreq;plot(PSP(1:end/2),Freq(1:end/2);title(频谱);