STFT短时傅立叶变换程序有问题说我mesh数据维度有问题，求指导

盛志鹏

1. 
   
2. function [Spec,Freq] = STFT(Sig,nLevel,WinLen,SampFreq)
   
3. %计算离散信号的短时Fourier变换
   
4. %   Sig：待分析信号
   
5. %nLevel:频率轴长度划分（默认值512）
   
6. %WinLen：汉宁窗长度（默认值64）
   
7. %SamFreq：信号的采样频率（默认值1）
   
8. 
   
9. 
   
10. if(nargin<1),
    
11.     error ('At least one parameter required!');
    
12. end;
    
13. 
    
14. Sig = real(Sig);
    
15. SigLen=length(Sig);
    
16. if(nargin<4),
    
17.     SampFreq = 1;
    
18. end
    
19. if(nargin<3),
    
20.     WinLen=64;
    
21. end
    
22. if (nargin<2),
    
23.     nLevel=513;
    
24. end
    
25. 
    
26. 
    
27. nLevel=ceil(nLevel/2)*2+1;
    
28. WinLen=ceil(WinLen/2)*2+1;
    
29. WinFun=exp(-6*linspace(-1,1,WinLen).^2);
    
30. WinFun=WinFun/norm(WinFun);
    
31. Lh=(WinLen-1)/2;
    
32. Ln=(nLevel-1)/2;
    
33. Spec=zeros(nLevel,SigLen);
    
34. 
    
35. wait=waitbar(0,'Under calculation,please wait...');
    
36. for iLoop =1:SigLen,
    
37.     waitbar(iLoop/SigLen,wait);
    
38.     iLeft=min([iLoop-1,Lh,Ln]);
    
39.     iRight=min([SigLen-iLoop,Lh,Ln]);
    
40.     iIndex=iLeft:iRight;
    
41.     iIndex1=iIndex+iLoop;
    
42.     iIndex2=iIndex+Lh+1;
    
43.     Index=iIndex+Ln+1;
    
44.     Spec(Index,iLoop)=Sig(iIndex1).*conj(WinFun(iIndex2));
    
45. end;
    
46. 
    
47. 
    
48. close(wait);
    
49. Spec=fft(Spec);
    
50. Spec=abs(Spec(1:(end-1)/2,:));
    
51. Freq=linspace(0,0.5,(nLevel-1)/2)*SampFreq;
    
52. t=(0:(SigLen))/SampFreq;
    
53. clf
    
54. set(gcf,'Position',[20 100 500 430]);
    
55. set(gcf,'Color','w');
    
56. axes('Position',[0.1 0.45 0.53 0.5]);
    
57. mesh(t,Freq,Spec);
    
58. axis([min(t) max(t) 0 max(Freq)]);
    
59. colorbar
    
60. xlabel('t/s');
    
61. ylabel('f/Hz');
    
62. title('STFT时频谱图');
    
63. axes('Position',[0.1 0.1 0.55 0.25]);
    
64. polt(t,Sig);
    
65. axis tight
    
66. ylabel('x(t)');
    
67. title('时域波形');
    
68. axes('Position',[0.73 0.45 0.24 0.5]);
    
69. PSP=abs(fft(Sig));
    
70. Freq=linspace(0,1,SigLen)*SampFreq;
    
71. plot(PSP(1:end/2),Freq(1:end/2));
    
72. tiltle('频谱');
    
73. 
    
74. 
    
75. 
    
76.    
    
77.         

复制代码

盛志鹏

1. SampFreq = 100;
   
2. t = 0:1/SampFreq:5;
   
3. Sig=sin(2*pi*(5*t+4*t.^2));
   
4. Sig=Sig+sin(2*pi*(45*t-4*t.^2));
   
5. STFT(Sig,512,128,100);
   

复制代码

上面的这个是要画的函数..

jiangwj

数据没有一一对应吧

梦时无意

学习了~~~

zghitboy

提示你索引超出维度了吗 注意数据矩阵的对应关系