带扩展系数的zfft

yangzj

一般来说,由于滤波器的过渡带的问题,ZFFT结果的两端幅值都会减小.
经过扩展后的ZFFT可解决这一问题.

1. function [xz fz]=zoomffta(x,fs,N,fe,D,a)
   
2. % 带扩展系数a的zoomfft
   
3. % x---输入时间序列,注意x的长度不能小于(N-1)*D+2*M；
   
4. % fs--x的采样频率；
   
5. % N---做谱点数；
   
6. % D---细化倍数；
   
7. % a---扩展系数。a取大于0小于等于1的数，越大M越小，计算量越小
   
8. 
   
9. M=4*D/a;                        %滤波器半阶数
   
10. k=1:M;                          
    
11. w=0.5+0.5*cos(pi*k/M);          %Hanning窗
    
12. 
    
13. % 求取带通滤波器上下界；
    
14. fl=max(fe-fs/(4*D),-fs/2);
    
15. fh=min(fe+fs/(4*D),fs/2);
    
16. 
    
17. %求取扩展带通滤波器上下界；
    
18. hfl=fl-(fh-fl)*a/2;
    
19. hfh=fh+(fh-fl)*a/2;
    
20. 
    
21. %构造扩展带通滤波器；
    
22. wl=2*pi*hfl/fs;
    
23. wh=2*pi*hfh/fs;
    
24. hr(1)=(wh-wl)/pi;
    
25. hr(2:M+1)=(sin(wh*k)-sin(wl*k))./(pi*k).*w;
    
26. hi(1)=0;
    
27. hi(2:M+1)=(cos(wl*k)-cos(wh*k))./(pi*k).*w;
    
28. 
    
29. %选抽滤波
    
30. for k=1:N
    
31.     kk=(k-1)*D+M;
    
32.     xrz(k)=x(kk+1)*hr(1)+sum(hr(2:M+1).*(x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));
    
33.     xiz(k)=x(kk+1)*hi(1)+sum(hi(2:M+1).*(-x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));
    
34. end
    
35. 
    
36. %移频，把fl移到0频
    
37. yf=D*fl/fs;                       %移频量
    
38. xz=(xrz+j*xiz).*exp(-j*2*pi*(0:N-1)*yf);
    
39. 
    
40. xz=fft(xz);
    
41. xz=xz(1:N/2)/N;                  %细化复数谱
    
42. fz=(0:N/2-1)*fs/N/D+fl;          %细化谱对应的频率

复制代码

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

yixiaofan

%构造扩展带通滤波器；
wl=2*pi*hfl/fs;
wh=2*pi*hfh/fs;
hr(1)=(wh-wl)/pi;
hr(2:M+1)=(sin(wh*k)-sin(wl*k))./(pi*k).*w;
hi(1)=0;
hi(2:M+1)=(cos(wl*k)-cos(wh*k))./(pi*k).*w;

%选抽滤波
for k=1:N
    kk=(k-1)*D+M;
    xrz(k)=x(kk+1)*hr(1)+sum(hr(2:M+1).*(x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));
    xiz(k)=x(kk+1)*hi(1)+sum(hi(2:M+1).*(-x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));
end


请问hr(1)和hi(1)为什么要单独求呢？
for循环里是做卷积吗？能否稍微解释具体点。
谢谢

yangzj

单独求是因为滤波器在0点是0/0的形式,应该用洛比塔法则来求.
for 里可以说是做卷积,滤波就是一个卷积过程.

yixiaofan

%选抽滤波
for k=1:N
    kk=(k-1)*D+M;
    xrz(k)=x(kk+1)*hr(1)+sum(hr(2:M+1).*(x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));
    xiz(k)=x(kk+1)*hi(1)+sum(hi(2:M+1).*(-x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));
end

这几句是怎样实现抽取和滤波，实在是看不懂了，能否解释的详细一些？
谢谢！

yangzj

　　原帖由 yixiaofan 于 2007-5-15 21:13 发表
　　%选抽滤波
　　for k=1:N
　　kk=(k-1)*D+M; %从第M点开始隔D点进行滤波
　　xrz(k)=x(kk+1)*hr(1)+sum(hr(2:M+1).*(x(kk+2:kk+M+1)+x(kk:-1:kk-M+1))); %滤波即做一个卷积,卷积对每一个点来说又是对滤波器进行反转再做内积运算,具体参考卷积的步骤;这里考滤到滤波器实部的偶对称,只用了一半来表示滤波器的实部;
　　xiz(k)=x(kk+1)*hi(1)+sum(hi(2:M+1).*(-x(kk+2:kk+M+1)+x(kk:-1:kk-M+1))); %虚部同样,只不过滤波器虚部是奇对称的
　　e ...

happyliu0712

%选抽滤波
for k=1:N
    kk=(k-1)*D+M;
    xrz(k)=x(kk+1)*hr(1)+sum(hr(2:M+1).*(x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));
    xiz(k)=x(kk+1)*hi(1)+sum(hi(2:M+1).*(-x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));
end
为什么报错啊？

happyliu0712

是我自己搞错了，导入数据出的错。

ruben052013

把抽取滤波和变频的顺序倒过来不会有什么问题吗？我看的很多的论文上面都是先变频后滤波的！这样是不是对滤波来说，要简单一些？先变频把关心的频率终点移到零频，这样应该需要的是两个低通滤波器而不是带通滤波器！如果先抽取滤波需要的是带通滤波器~

[ 本帖最后由 ruben052013 于 2009-3-6 15:28 编辑 ]