马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
一般来说,由于滤波器的过渡带的问题,ZFFT结果的两端幅值都会减小.
经过扩展后的ZFFT可解决这一问题.
- function [xz fz]=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 编辑 ] |