带扩展系数的zfft
一般来说,由于滤波器的过渡带的问题,ZFFT结果的两端幅值都会减小.经过扩展后的ZFFT可解决这一问题.
function =zoomffta(x,fs,N,fe,D,a)
% 带扩展系数a的zoomfft
% x---输入时间序列,注意x的长度不能小于(N-1)*D+2*M;
% fs--x的采样频率;
% N---做谱点数;
% D---细化倍数;
% a---扩展系数。a取大于0小于等于1的数,越大M越小,计算量越小
M=4*D/a; %滤波器半阶数
k=1:M;
w=0.5+0.5*cos(pi*k/M); %Hanning窗
% 求取带通滤波器上下界;
fl=max(fe-fs/(4*D),-fs/2);
fh=min(fe+fs/(4*D),fs/2);
%求取扩展带通滤波器上下界;
hfl=fl-(fh-fl)*a/2;
hfh=fh+(fh-fl)*a/2;
%构造扩展带通滤波器;
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
%移频,把fl移到0频
yf=D*fl/fs; %移频量
xz=(xrz+j*xiz).*exp(-j*2*pi*(0:N-1)*yf);
xz=fft(xz);
xz=xz(1:N/2)/N; %细化复数谱
fz=(0:N/2-1)*fs/N/D+fl; %细化谱对应的频率
[ 本帖最后由 MVH 于 2006-10-4 15:51 编辑 ] %构造扩展带通滤波器;
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循环里是做卷积吗?能否稍微解释具体点。
谢谢 单独求是因为滤波器在0点是0/0的形式,应该用洛比塔法则来求.
for 里可以说是做卷积,滤波就是一个卷积过程. %选抽滤波
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
这几句是怎样实现抽取和滤波,实在是看不懂了,能否解释的详细一些?
谢谢! 本帖最后由 wdhd 于 2016-9-14 10:31 编辑
原帖由 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 ...
回复 5楼 的帖子
%选抽滤波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
为什么报错啊?
回复 6楼 的帖子
是我自己搞错了,导入数据出的错。 把抽取滤波和变频的顺序倒过来不会有什么问题吗?我看的很多的论文上面都是先变频后滤波的!这样是不是对滤波来说,要简单一些?先变频把关心的频率终点移到零频,这样应该需要的是两个低通滤波器而不是带通滤波器!如果先抽取滤波需要的是带通滤波器~[ 本帖最后由 ruben052013 于 2009-3-6 15:28 编辑 ]
页:
[1]