武汉理工大学 基于矩形窗、三角窗、海明窗、汉宁窗、布拉克曼窗的FIR数字滤波器设计 下载本文

子程序:

function hd=ideal_hp1(Wc,N) alpha=(N-1)/2; n=0:1:N-1; m=n-alpha+eps;

%hd=[sin(pi*m)-sin(Wc*m)]/(pi*m); hd=[sin(pi*m)-sin(Wc*m)]./(pi*m);

function[db,mag,pha,w]=freqz_m2(b,a) [H,w]=freqz(b,a,1000,'whole'); H=(H(1:1:501))'; w=(w(1:1:501))'; mag=abs(H);

db=20*log10((mag+eps)/max(mag)); pha=angle(H);

主程序:

clear all; Ws=6.5/15*pi Wp=3.5/15*pi

tr_width=Ws-Wp;%过渡带宽度

N=ceil(8*pi/tr_width);%滤波器长度 n=0:1:N-1;

Wc=(Ws+Wp)/2;%理想滤波器的截止频率

hd=ideal_hp1(Wc,N);%理想滤波器的单位冲击响应

w_ham=(triang(N))'; string=['三角窗','N=',num2str(N)];

h=hd.*w_ham;%截取得到的实际的单位脉冲响应

[db,mag,pha,w]=freqz_m2(h,[1]);%计算实际滤波器的幅度响应 delta_w=2*pi/1000;

Ap=-(min(db(Ws/delta_w+1:1:501)));%实际通带波纹 As=-round(max(db(1:1:Wp/delta_w+1)));%实际阻带波纹 subplot(3,2,1);stem(n,hd);title('理想脉冲响应hd(n)'); axis([0,N-1,-0.5,0.5]);xlabel('n');ylabel('hd(n)'); subplot(3,2,2);stem(n,w_ham);axis([0,N-1,0,1.1]); xlabel('n');ylabel('w(n)');text(1.5,1.3,string);

subplot(3,2,3);stem(n,h);title('实际脉冲响应hd(n)'); axis([0,N,-1.4,1.4]);xlabel('n');ylabel('hd(n)'); subplot(3,2,4);plot(w,pha);title('想频特性');

axis([0,3.15,-4,4]);xlabel('频率(pi)');ylabel('相位(o)'); subplot(3,2,5);plot(w/pi,db);title('幅度特性(dB)');

axis([0,1,-100,10]);xlabel('频率(pi)');ylabel('(分贝数)'); subplot(3,2,6);plot(w,mag);title('频率特性');

axis([0,2.95,0,1.5]);xlabel('频率(rad)');ylabel('幅值'); fs=12000; t=(0:100)/fs;

x=sin(2*pi*t*750)+sin(2*pi*t*1500)+sin(2*pi*t*3000); q=filter(h,1,x); [a,f1]=freqz(x); f1=f1/pi*fs/2; [b,f2]=freqz(q); f2=f2/pi*fs/2; figure(2);

subplot(2,1,1); plot(f1,abs(a));

title('输入波形频谱图');

xlabel('频率');ylabel('幅度') subplot(2,1,2); plot(f2,abs(b));

title('输出波形频谱图');

xlabel('频率');ylabel('幅度') 检验程序:

用含多种频率成份离散时间信号作输入信号,得到通过滤波器后与输出信号的对比结果

fs=20000;

t=(0:100)/fs;

x=sin(2*pi*t*1500)+sin(2*pi*t*4000)+sin(2*pi*t*8000);%输入信号 q=filter(h,1,x);%滤波处理,得到输出信号 [a,f1]=freqz(x); f1=f1/pi*fs/2; [b,f2]=freqz(q); f2=f2/pi*fs/2; figure(2);

subplot(2,1,1); plot(f1,abs(a));

title('输入波形频谱图'); xlabel('频率'); ylabel('幅度'); subplot(2,1,2); plot(f2,abs(b));

title('输出波形频谱图'); xlabel('频率'); ylabel('幅度')

图形示例:低通部分

检验图形:低通部分