华工【数字信号处理】(DSP)实验三四 FIR,IIR滤波器设计

一、 用窗函数法设计如下三种类型的FIR滤波器，滤波器的性能指标：
低通滤波器：通带截止频率1000Hz，阻带截止频率1200Hz，采样率4000Hz，通带波纹1dB，最小阻带衰减40dB；
高通滤波器：通带截止频率1400Hz，阻带截止频率1000Hz，采样率4000Hz，通带波纹1dB，最小阻带衰减40dB；
带通滤波器：通带截止频率fp1=900Hz, fp2=1300Hz，阻带截止频率fs1=600Hz， fs2=1500Hz, 采样率4000Hz，通带波纹1dB，最小阻带衰减40dB。
二、 再用双线性变换法设计上面的三种IIR滤波器（可以选择butteworth or chebyshey）
三、要画出各种滤波器的频率响应特性（用freqz函数）

FIR代码：
1.FIR低通

%FIR lowpass
close all;
clear all;
ft=4000;    %采样频率
%课本p356
wp=(2*pi*1000)/ft;
ws=(2*pi*1200)/ft;
bowen=1;    %通带波纹
%40dB,采用hann窗
M=3.11*pi/(ws-wp);
M=ceil(M);
N=2*M+1;    %滤波器长度
W=hann(N);
wc=(ws+wp)/2/pi;  %截止频率
h=fir1(N-1,wc,'low',W);    %FIR低通滤波器冲激响应
[H,w]=freqz(h);     %低通滤波器的频率响应
H=abs(H);
f=ft*w/2/pi;
plot(f,H)
xlabel('f/MHz');
ylabel('幅度');
title('FIR低通频谱图');

2.FIR高通

ws=(2*pi*1000)/ft;
bowen=1;    %通带波纹
%40dB,采用hann窗
M=3.11*pi/(wp-ws);
M=ceil(M);
N=2*M+1;    %滤波器长度
W=hann(N);
wc=(ws+wp)/2/pi;  %截止频率
h=fir1(N-1,wc,'high',W);    %FIR低通滤波器冲激响应
[H,w]=freqz(h);     %高通滤波器的频率响应
H=abs(H);
f=ft*w/2/pi;
plot(f,H)
xlabel('f/MHz');
ylabel('幅度');
title('FIR高通频谱图');

3.FIR带通

%FIR bandpass
%通带截止频率fp1=900Hz, fp2=1300Hz，阻带截止频率fs1=600Hz， fs2=1500Hz,
%采样率4000Hz，通带波纹1dB，最小阻带衰减40dB。
close all;
clear all;
ft=4000;    %采样频率
wp1=(2*pi*900)/ft;
wp2=(2*pi*1300)/ft;
ws1=(2*pi*600)/ft;
ws2=(2*pi*1500)/ft;
bowen=1;    %通带波纹
%40dB,采用hann窗
 %滤波器的阶数估计
if((wp1-ws1)<(ws2-wp2))
    M=3.11*pi/(wp1-ws1);
else M=3.11*pi/(ws2-wp2);
end
M=ceil(M);
N=2*M+1;    %滤波器长度
W=hann(N);
wc1=(ws1+wp1)/2/pi;  %截止频率
wc2=(ws2+wp2)/2/pi;  
h=fir1(N-1,[wc1,wc2],W);    %FIR低通滤波器冲激响应
[H,w]=freqz(h);     %带通滤波器的频率响应
H=abs(H);
f=ft*w/2/pi;
plot(f,H)
xlabel('f/MHz');
ylabel('幅度');
title('FIR带通频谱图');

IIR代码
1.IIR低通

%IIR lowpass
close all;
clear all;
ft=4000;    %采样频率
rp_db=1;
rs_db=40;
%rp:通带波纹;rs:最小阻带衰减
rp=10^(-rp_db/20);
rs=10^(-rs_db/20);
%课本p356
wp=(2*pi*1000)/ft;
ws=(2*pi*1200)/ft;
%预畸变
wp=tan(wp/2);
ws=tan(ws/2);
 
[N,wn]=buttord(wp,ws,rp_db,rs_db,'s')   ;  %巴特沃斯低通滤波器
[num,den]=butter(N,wn,'low','s');
[numz,denz]=bilinear(num,den,0.5);      %双线性变换，采样频率Fs取0.5
[H,w]=freqz(numz,denz);     %低通滤波器的频率响应
H=abs(H);
f=ft*w/2/pi;
plot(f,H)
xlabel('f/MHz');
ylabel('幅度');
title('IIR低通频谱图');

2.IIR高通

%IIR highpass
close all;
clear all;
ft=4000;    %采样频率
rp_db=1;
rs_db=40;
%rp:通带波纹;rs:最小阻带衰减
rp=10^(-rp_db/20);
rs=10^(-rs_db/20);
%课本p356
wp=(2*pi*1400)/ft;
ws=(2*pi*1000)/ft;
%预畸变
wp=tan(wp/2);
ws=tan(ws/2);
 
[N,wn]=cheb1ord(ws,wp,rp_db,rs_db,'s');
[num,den]=cheby1(N,1,wn,'s');
[numT,denT]=lp2hp(num,den,wp);    %lp2hp低通转高通
[numz,denz]=bilinear(numT,denT,0.5);    %双线性变换，此处采样频率Fs取0.5
[H,w]=freqz(numz,denz);     %高通滤波器的频率响应
H=abs(H);
f=ft*w/2/pi;
plot(f,H)
xlabel('f/MHz');
ylabel('幅度');
title('IIR高通频谱图');

3.IIR带通

%IIR bandpass
close all;
clear all;
%通带截止频率fp1=900Hz, fp2=1300Hz，阻带截止频率fs1=600Hz， fs2=1500Hz,
%采样率4000Hz，通带波纹1dB，最小阻带衰减40dB。
ft=4000;    %采样频率
rp_db=1;
rs_db=40;
%rp:通带波纹;rs:最小阻带衰减
rp=10^(-rp_db/20);
rs=10^(-rs_db/20);
%课本p356
%归一化角频率
wp1=(2*pi*900)/ft;
wp2=(2*pi*1300)/ft;
ws1=(2*pi*600)/ft;
ws2=(2*pi*1500)/ft;
%预畸变
wp1=tan(wp1/2);
wp2=tan(wp2/2);
ws1=tan(ws1/2);
ws2=tan(ws2/2);
%通带带宽
BW=wp2-wp1;
%通带边界乘积
wp=wp1*wp2;
%阻带边界乘积
ws=ws1*ws2;
%通带几何平均
wo=wp^(1/2);
 
if (wp==ws)           %判断wp==ws==wo^2
    ws11=ws1;           %相等则不需要修改ws1
else
    ws11=ws/ws2;        %不相等则求出新的ws1
end
 %低通滤波器的阻带边界频率
 %课本p629
 %归一化wp=1;
wstop=(wp-ws11^2)/(BW*ws11);   
[N,wn]=buttord(1,wstop,rp_db,rs_db,'s');           %巴特沃斯低通滤波器
[B,A]=butter(N,wn,'s');
[BT,AT]=lp2bp(B, A,wo,BW);          %低通滤波器转换成带通滤波器
[numz,denz]=bilinear(BT,AT,0.5);    %双线性变换，此处采样频率Fs取0.5
[H,w]=freqz(numz,denz);     %带通滤波器的频率响应
H=abs(H);
f=ft*w/2/pi;
plot(f,H)
xlabel('f/MHz');
ylabel('幅度');
title('IIR带通频谱图');

分析:
FIR低通滤波器的阶数N=65
FIR高通滤波器的阶数N=33
FIR带通滤波器的阶数N=65
可见,FIR滤波器的阶数都比较大,一般大于10
IIR低通滤波器的阶数N=17
IIR高通滤波器的阶数N=5
IIR带通滤波器的阶数N=5
可见,IIR滤波器的阶数会比较低,除个别外一般小于10

此外,从频域波形上看,FIR和IIR滤波器没有明显的区别