[原创]zoomfft

yangzj

1. %ZoomFFT谱
   
2. %x-信号序列
   
3. %fs-采样频率
   
4. %N-做谱点数
   
5. %fe-分析中心频率
   
6. %D-细化倍数
   
7. %L-平均段数
   
8. %M-滤波器半阶数
   
9. %f-返回频率向量
   
10. %xz-返回幅值谱
    
11. 
    
12. function [f xz]=ZoomFFT(x,fs,N,fe,D,L,M)
    
13. 
    
14. k=1:M;                          
    
15. w=0.5+0.5*cos(pi*k/M);          %Hanning窗
    
16. 
    
17. fl=max(fe-fs/(4*D),-fs/2.2);
    
18. fh=min(fe+fs/(4*D),fs/2.2);
    
19. 
    
20. yf=D*fl;                       %移频量
    
21. df=fs/D/N;
    
22. f=fl:df:fl+(N/2-1)*df;
    
23. xz=zeros(1,N/2);
    
24. wl=2*pi*fl/fs;
    
25. wh=2*pi*fh/fs;
    
26. hr(1)=(wl-wh)/pi;
    
27. hr(2:M+1)=(sin(wl*k)-sin(wh*k))./(pi*k).*w;
    
28. hi(1)=0;
    
29. hi(2:M+1)=(cos(wl*k)-cos(wh*k))./(pi*k).*w;
    
30. 
    
31. k=0:N-1;
    
32. w=0.5-0.5*cos(2*pi*k/N);
    
33. 
    
34. for i=1:L
    
35.     for k=1:N
    
36.         kk=(k-1)*D+M+(i-1)*N;
    
37.         xrz(k)=x(kk+1)*hr(1)+sum(hr(2:M+1).*(x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));
    
38.         xiz(k)=x(kk+1)*hi(1)+sum(hi(2:M+1).*(x(kk+2:kk+M+1)-x(kk:-1:kk-M+1)));
    
39.     end
    
40.     xzt=(xrz+j*xiz).*exp(-j*2*pi*(0:N-1)*yf/fs);
    
41.     xzt=xzt.*w;
    
42.     xzt=xzt-sum(xzt)/N;
    
43.     xzt=fft(xzt);
    
44.     xz=xz+(abs(xzt(1:N/2))/N*2).^2;
    
45. end
    
46. xz=(xz/L).^0.5;

复制代码

[ 本帖最后由 MVH 于 2006-10-4 15:52 编辑 ]

shanghai

1. for i=1:L
   
2. for k=1:N
   
3. kk=(k-1)*D+M+(i-1)*N;
   
4. xrz(k)=x(kk+1)*hr(1)+sum(hr(2:M+1).*(x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));
   
5. xiz(k)=x(kk+1)*hi(1)+sum(hi(2:M+1).*(x(kk+2:kk+M+1)-x(kk:-1:kk-M+1)));
   
6. end
   
7. xzt=(xrz+j*xiz).*exp(-j*2*pi*(0:N-1)*yf/fs);
   
8. xzt=xzt.*w;
   
9. xzt=xzt-sum(xzt)/N;
   
10. xzt=fft(xzt);
    
11. xz=xz+(abs(xzt(1:N/2))/N*2).^2;
    
12. end

复制代码

上面ZoomFFT里面这段代码没看明白, j指的是什么?

[ 本帖最后由 MVH 于 2006-10-4 15:53 编辑 ]

yangzj

没定义i,j的话,它们都是虚数单位

shanghai

不好意思,还要问一下,语句 xz=(xz/L).^0.5; 中xz是复数,求它的平方根赋给xz ,这样做的意义是什么? 谢谢,谢谢!

yangzj

xz不是复数,而是功率谱,这里开方是求的幅值谱,如果要求其他的谱,可做相应的改动

[ 本帖最后由 MVH 于 2006-10-4 15:53 编辑 ]

shanghai

1. #define N 2048
   
2. #define pi 3.1415926
   
3. #define M 10
   
4. void CAnalysisDlg::ZoomFFT(double *x, double fs, double fe, int D, int L)
   
5. {
   
6. double wl,wh,fl,fh,yf,df,c;
   
7. int kk;
   
8. complex N1,L1,w2,x1,x2,sumxzt;
   
9. 
   
10. int i,j;
    
11. 
    
12. for(i=0;i<M;i++)
    
13. {
    
14. w0=0.5+0.5*cos(pi*i/M);
    
15. }
    
16. fl=(fe-fs/(4*D))>(-fs/2.2) ? (fe-fs/(4*D)):(-fs/2.2);
    
17. fh=(fe+fs/(4*D))<(fs/2.2) ? (fe+fs/(4*D)):(fs/2.2);
    
18. yf=D*fl;
    
19. df=fs/D/N;
    
20. 
    
21. for(i=0;i<=(N/2-1);i++)
    
22. {
    
23. f=fl+i*df;
    
24. }
    
25. 
    
26. wl=2*pi*fl/fs;
    
27. wh=2*pi*fh/fs;
    
28. 
    
29. hr[0]=(wl-wh)/pi;
    
30. for(i=1;i<=M;i++)
    
31. {
    
32. hr=(sin(wl*i)-sin(wh*i))/(pi*i)*w0[i-1];
    
33. }
    
34. hi[0]=0;
    
35. for(i=1;i<=M;i++)
    
36. {
    
37. hi=(cos(wl*i)-cos(wh*i))/(pi*i)*w0[i-1];
    
38. }
    
39. 
    
40. for(i=0;i<N;i++)
    
41. {
    
42. w1=0.5-0.5*cos(2*pi*i/N);
    
43. }
    
44. 
    
45. for(i=0;i<L;i++)
    
46. {
    
47. for(j=0;j<N;j++)
    
48. {
    
49. kk=j*D+M+i*N;
    
50. for(int m=1;m<=M;m++)
    
51. {
    
52. sumhr=sumhr+hr[m]*(x[kk+1+m]+x[kk-m+1]);
    
53. sumhi=sumhi+hi[m]*(x[kk+1+m]-x[kk-m+1]);
    
54. }
    
55. xrz[j]=x[kk+1]*hr[0]+sumhr;
    
56. xiz[j]=x[kk+1]*hi[0]+sumhi;
    
57. }
    
58. }
    
59. N1=complex(N,0);
    
60. for(int k=0;k<N;k++)
    
61. {
    
62. x1=complex(xrz[k],xiz[k]);
    
63. x2=complex(cos(2*pi*k*yf/fs),-sin(2*pi*k*yf/fs));
    
64. xzt[k]=x1*x2;
    
65. w2=complex(w1[k],0);
    
66. xzt[k]=xzt[k]*w2;
    
67. sumxzt=sumxzt+xzt[k];
    
68. xzt[k]=xzt[k]-sumxzt/N1;
    
69. }
    
70. fft(xzt,11);
    
71. for(i=0;i<N/2;i++)
    
72. {
    
73. xz=(xzt.Real()*xzt.Real()+xzt.Imag()*xzt.Imag())/(1024*1024);
    
74. }
    
75. 
    
76. for(i=0;i<N/2;i++)
    
77. {
    
78. xz=sqrt(xz/L);
    
79. }
    
80. }

复制代码

该程序是我把上面的Matlab版 ZoomFFT改写成Vc++版,该程序通过编译没有问题,我调用函数ZoomFFT(dc,100000.0,3125.0,8,5),但是其中的xiz,xrz始终为0,所以后面的xz当然也是0了,请问这个函数错误在哪里,请指点指点,谢谢!!!

[ 本帖最后由 MVH 于 2006-10-4 15:54 编辑 ]

shanghai

下标问题我已经注意到了,Vc++中 Hanning窗的一半,是0到M-1,但还是不行.     请问有没有ZoomFFT算法方面的资料,仔细研究一下,谢谢!

yangzj

我采用的是丁康教授提出的基于复解析带通滤波器的复调制细化选带频谱分析，你可以参考这两篇论文:
[1] 丁康,谢明等. 基于复解析带通滤波器的复调制细化频谱分析方法研究[J]. 振动工程学报. 2001, 14(1):29~35.



[2] 谢明,丁康. 基于复解析带通滤波器的复调制细化谱分析的算法研究[J]. 振动工程学报. 2002, 15(4):479~483.

[ 本帖最后由 MVH 于 2006-10-4 15:55 编辑 ]

yangzj

1. for(i=0;i<M;i++)
   
2. {
   
3. w0=0.5+0.5*cos(pi*i/M);
   
4. }
   
5. 应该改成
   
6. for(i=0;i<M;i++)
   
7. {
   
8. w0=0.5+0.5*cos(pi*(i+1)/M);
   
9. }  

复制代码

[ 本帖最后由 MVH 于 2006-10-4 15:55 编辑 ]

shanghai

谢谢,我再看看这方面的资料!

shanghai

请问你的ZoomFFT 中平均段数L是什么意思,丁康教授资料中并没有提到这个参数呀!

yangzj

哦,这个是为了减少实际工程信号中随机信号的影响而取了多段信号,然后在功率谱上做平均

[ 本帖最后由 MVH 于 2006-10-4 15:55 编辑 ]

snort

丁康教授教授的论文我也下了，好好看看，向大家学习。

liuhuanjun

滤波器的半阶数是什么意思啊

liuhuanjun

kk=(k-1)*D+M+(i-1)*N;是什么意思？
这个程序运行有问题啊？