[原创]zoomfft

引用 shanghai 2006-5-26 10:40

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 2006-5-26 10:47

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

引用 shanghai 2006-5-29 21:52

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

引用 yangzj 2006-5-31 14:04

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

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

引用 shanghai 2006-6-2 10:31

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 2006-6-3 15:01

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

引用 yangzj 2006-6-3 17:31

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



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

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

引用 yangzj 2006-6-3 17:34

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 2006-6-3 19:46

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

引用 shanghai 2006-6-4 20:10

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

引用 yangzj 2006-6-10 12:06

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

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

引用 snort 2006-7-23 09:47

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

引用 liuhuanjun 2006-8-15 21:35

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

引用 liuhuanjun 2006-8-15 23:19

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

引用 miao7mijao 2006-9-22 08:07

现在正在学习这个程序呢，希望已经看明白的人能给我这个正徘徊在此的人一些提示，或说把你们的经验能多少的告诉我；先谢谢了！

引用 miao7mijao 2006-9-22 14:07

程序中的
fl=max(fe-fs/(4*D),-fs/2.2);
fh=min(fe+fs/(4*D),fs/2.2);
是怎么回事？

引用 miao7mijao 2006-9-22 16:16

关于第2个海宁窗公式w=0.5-0.5*cos(2*pi*k/N)
我知道了是单边的，但是第1个k=1：M
w=0.5+0.5*cos(pi*k/M)
而书上为k=-N/2，。。。N/2；
w=0.5+0.5*cos(2pi*k/N)；是不是只要保证k的个数与N的个数相等就可以了
也许我问的问题让大家见笑了，但是实在找不到答案了！
谢谢了！

引用 miao7mijao 2006-9-22 19:08

我这提示：xrz(k)=x(kk+1)*hr(1)+sum(hr(2:M+1).*(x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));这一句超过矩阵界限
那你测试的结果对吗?
我也运行了,也出现类似的问题?该怎么解决?

引用 songzy41 2006-9-23 06:51

原帖由 miao7mijao 于 2006-9-22 16:16 发表
关于第2个海宁窗公式w=0.5-0.5*cos(2*pi*k/N)
我知道了是单边的，但是第1个k=1：M
w=0.5+0.5*cos(pi*k/M)
而书上为k=-N/2，。。。N/2；
w=0.5+0.5*cos(2pi*k/N)；是不是只要保证k的个数与N的个数相等就可以了 ...

对于海宁窗的表达式:
w=0.5+0.5*cos(pi*k/N)
适用于k=-N/2，。。。N/2；而
w=0.5-0.5*cos(2*pi*k/N)
适用于k=0：N-1。这是由笫1式在时间上平移了N/2个点。因为我们实际中信号都是N长，时间上是Δt,2Δt,...NΔt（Δt是采样周期），加海宁窗时用笫2式。

查看全部评论(65)