带通带阻高通低通滤波器的设计(IIR)

一：带通滤波器的设计 %周期方波时域 f1=72; fs=10000; t=0:1/fs:0.1;

x1=square(2*pi*f1*t); subplot(321); plot(t,x1);

axis([0,0.1,-1.3,1.3]); xlabel('t/s'); ylabel('x(t)'); title('时域谱'); %周期方波频域 L=1250;

X1=fft(x1,L); ws=2*pi*fs; w=(0:L-1)*ws/L; subplot(322);

stem((0:L-1),abs(X1)); xlabel('w/rad'); ylabel('X'); title('幅度谱');

%采用BW和双线性变换法 Wp1=0.1*pi,Wp2=0.5*pi; Ws1=0.01*pi,Ws2=0.6*pi; Ap=1;As=50;

wp1=2*fs*tan(Wp1/2),wp2=2*fs*tan(Wp2/2); ws1=2*fs*tan(Ws1/2),ws2=2*fs*tan(Ws2/2); wp=[wp1 wp2],ws=[ws1 ws2]; [N,wc]=buttord(wp,ws,Ap,As,'s'); [num,den]=butter(N,wc,'s');

[numd,dend]=bilinear(num,den,fs); w=linspace(0,pi,1024); H=freqz(numd,dend,w); subplot(324); plot(w/pi,abs(H)); xlabel('W'); ylabel('H');

title('幅度响应');

[h,t1]=impz(numd,dend); subplot(323); plot(t1,h);

axis([0,60,min(h),max(h)]); xlabel('t/s'); ylabel('h');

title('单位冲击响应'); %输出信号

y=filter(numd,dend,x1); subplot(325); plot(t,y);

axis([0,0.1,-1.3,1.3]); Y=fft(y,L); ws=2*pi*fs; w=(0:L-1)*ws/L; subplot(326);

stem((0:L-1),abs(Y));

二：带阻滤波器的设计 %周期方波时域 f1=64; fs=10000; t=0:1/fs:0.1;

x1=square(2*pi*f1*t); subplot(321); plot(t,x1);

axis([0,0.1,-1.3,1.3]); xlabel('t/s'); ylabel('x(t)'); title('时域谱') %周期方波频域 L=1250;

X1=fft(x1,L); ws=2*pi*fs; w=(0:L-1)*ws/L; subplot(322);

stem((0:L-1),abs(X1)); xlabel('w/rad'); ylabel('X'); title('幅度谱');

%采用BW和双线性变换法 Wp1=0.1*pi,Wp2=0.4*pi; Ws1=0.2*pi,Ws2=0.3*pi; Ap=1;As=50;

wp1=2*fs*tan(Wp1/2),wp2=2*fs*tan(Wp2/2); ws1=2*fs*tan(Ws1/2),ws2=2*fs*tan(Ws2/2); B=ws2-ws1,w0=sqrt(ws1*ws2);

ws=1;

wp1=(B*wp1)/(-wp1^2+w0^2),wp2=(B*wp2)/(-wp2^2+w0^2); wp=max(abs(wp1),abs(wp2)); %wp=[wp1 wp2],ws=[ws1 ws2]; [N,wc]=buttord(wp,ws,Ap,As,'s'); [num,den]=butter(N,wc,'s');

[numt,dent]=lp2bs(num,den,w0,B); [numd,dend]=bilinear(numt,dent,fs); w=linspace(0,pi,1024); H=freqz(numd,dend,w); subplot(324);

plot(w/pi,20*log10(abs(H))); grid;

xlabel('W'); ylabel('H/db'); title('幅度响应');

[h,t1]=impz(numd,dend); subplot(323); plot(t1,h);

axis([0,60,min(h),max(h)]); xlabel('t/s'); ylabel('h');

title('单位冲击响应'); %输出信号

y=filter(numd,dend,x1); subplot(325); plot(t,y);

axis([0,0.1,-1.3,1.3]); Y=fft(y,L); ws=2*pi*fs; w=(0:L-1)*ws/L; subplot(326);

stem((0:L-1),abs(Y));

三．高通滤波器的设计 %周期方波时域 f1=64; fs=10000; t=0:1/fs:0.1;

x1=square(2*pi*f1*t); subplot(321); plot(t,x1);

axis([0,0.1,-1.3,1.3]); xlabel('t/s');

ylabel('x(t)'); title('时域谱') %周期方波频域 L=1250;

X1=fft(x1,L); ws=2*pi*fs; w=(0:L-1)*ws/L; subplot(322);

stem((0:L-1),abs(X1)); xlabel('w/rad'); ylabel('X'); title('幅度谱');

%采用BW和双线性变换法 Wp=0.5*pi,Ws=0.1*pi; Ap=1,As=50;

wp=2*fs*tan(Wp/2);ws=2*fs*tan(Ws/2); wp1=1/wp,ws1=1/ws;

[N,wc]=buttord(wp1,ws1,Ap,As,'s'); [num,den]=butter(N,wc,'s'); [numt,dent]=lp2hp(num,den,1); [numd,dend]=bilinear(numt,dent,fs); w=linspace(0,pi,1024); H=freqz(numd,dend,w); subplot(324); plot(w/pi,abs(H)); xlabel('W'); ylabel('H');

title('幅度响应');

[h,t1]=impz(numd,dend); subplot(323); plot(t1,h);

axis([0,max(t1),min(h),max(h)]); xlabel('t/s'); ylabel('h');

title('单位冲击响应'); %输出信号

y=filter(numd,dend,x1); subplot(325); plot(t,y);

axis([0,0.1,-1.3,1.3]); Y=fft(y,L);

w=(0:L-1)*ws/L; subplot(326);

stem((0:L-1),abs(Y));

四：低通滤波器的设计 %周期方波时域 f1=64; fs=10000; t=0:1/fs:0.1;

x1=square(2*pi*f1*t); subplot(321); plot(t,x1);

axis([0,0.1,-1.3,1.3]); xlabel('t/s'); ylabel('x(t)'); title('时域谱') %周期方波频域 L=1250;

X1=fft(x1,L); ws=2*pi*fs;

w1=(0:L-1)*ws/pi; subplot(322);

stem((0:L-1),abs(X1)); xlabel('w/rad'); ylabel('X'); title('幅度谱');

%采用BW和双线性变换法 Wp=0.1*pi;Ws=0.5*pi; Ap=1;As=50;

wp=2*fs*tan(Wp/2);ws=2*fs*tan(Ws/2); [N,wc]=buttord(wp,ws,Ap,As,'s'); [num,den]=butter(N,wc,'s');

[numd,dend]=bilinear(num,den,fs); w=linspace(0,pi,1024); H=freqz(numd,dend,w); subplot(324); plot(w/pi,abs(H)); xlabel('W'); ylabel('H');

title('幅度响应');

[h,t1]=impz(numd,dend); subplot(323); plot(t1,h);

axis([0,max(t1),min(h),max(h)]); xlabel('t/s'); ylabel('h');

title('单位冲击响应'); %输出信号

y=filter(numd,dend,x1); subplot(325); plot(t,y);

axis([0,0.1,-1.3,1.3]); Y=fft(y,L); ws=2*pi*fs;

w=(0:L-1)*ws/pi; subplot(326);

stem((0:L-1),abs(Y));