DSP实验5 IIR滤波器的设计与信号滤波

实验五IIR滤波器的设计与信号滤波一、实验目的设计IIR滤波器，实现对存在加性干扰的时域离散信号进行滤波。二、实验内容已知带加性干扰的信号用x(n)表示，x(n)=xs(n)+η(n)，式中xs(n)是有用的信号，是一个0～0.2πrad的带限信号。η(n)是一个干扰信号，其频谱分布在0.3πrad以上。要求设计一个巴特沃斯IIR数字滤波器对信号x(n)进行滤波，将干扰η(n)滤除。要求在xs(n)所在的通带内滤波器幅度平坦，在0.2πrad处幅度衰减不大于1dB，在噪声所在的0.3πrad以上的频带内滤波器幅度衰减大于等于40dB。三、实验步骤1.根据题目要求确定要设计的数字滤波器的技术指标(低通滤波器指标：wp=0.2πrad,ws=0.3πrad,αp=1dB,αs=40dB)；2.用双线性变换法频率转换公式，将上面的技术指标转换为巴特沃斯模拟滤波器的技术指标；3.调用MATLAB函数buttord和butter，设计该模拟滤波器；4.用双线性变换法将设计的模拟滤波器转换为数字滤波器，求出该滤波器的系统函数并画出它的直接型网络结构图；5.画出该滤波器的幅度特性曲线；6.根据1所确定的技术指标，调用MATLAB函数buttord和butter，直接设计数字滤波器，观察设计结果与上面用双线性变换法的设计结果是否相同。8.滤波仿真：调用MATLAB工具箱函数filter对下面给出的带加性干扰的信号x(n)进行滤波，观察滤波效果(绘制滤波前后信号的时域和频域波形)。实验信号为心电图信号采样序列，人体心电图信号在测量过程中往往受到工业高频干扰，所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。下面给出一实际心电图信号采样序列样本x(n)，其中存在高频干扰。在实验中，以x(n)作为输入序列，滤除其中的干扰成分。x(n)：={-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0}四、实验用MATLAB函数简介filter功能：一维数字滤波器直接Ⅱ型实现函数。调用格式：yn=filter(B,A,xn)：按直接Ⅱ型实现结构对输入信号向量xn滤波，返回滤波器输出信号向量yn，调用参数B和A分别为滤波器系统函数的分子和分母多项式系数向量。其实质是求解差分方程：A(1)y(n)=B(1)x(n)+B(2)x(n-1)+…+B(M+1)x(n-M)-A(2)y(n-1)-…-A(N+1)y(n-N如果A(1)不等于1时，则对系数关于A(1)归一化后计算输出信号y(n)。当A=1时，对应IIR-DF的直接Ⅱ型实现。五、实验报告要求1.写出所设计的数字滤波器的主要技术指标、系统函数、网络结构；2.写出两种IIR数字滤波器的设计原理、步骤；3.写出仿真滤波试验的结果。附录1：用双线性变换法设计IIR数字滤波器的程序：%输入信号及其波形x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,...0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,...4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];n=0:55;subplot(2,2,1);stem(n,x,'.');axis([060-10050]);holdon;n=0:60;m=zeros(61);plot(n,m);xlabel('n');ylabel('x(n)');title('心电图信号采样序列x(n)');X=fft(x,256);subplot(2,2,2)k=0:255;plot(2*pi*k/256,abs(X))xlabel('k');ylabel('X(k)');title('心电图信号采样序列的频谱');%用双线性变换法滤波器设计IIR-DFT=1;Wp=0.2*pi;Ws=0.3*pi;Rp=1;Rs=40;Omegap=2*tan(Wp/2)/T;%模拟波器通带截止频率Omegas=2*tan(Ws/2)/T;%模拟波器通带截止频率Fs=1/T;%采样频率[N,Wc]=buttord(Omegap,Omegas,Rp,Rs,'s');%计算巴特沃斯模拟滤波器阶数N及%3dB截止频率Wc[Z,P,K]=buttap(N);%求butterworth归一化原型H(p)，零极点增益形式[b,a]=zp2tf(Z,P,K);%转换成归化系统函数形式，H(p)为分子分母多项式形式[bt,at]=lp2lp(b,a,Wc);%去归一化转换为实际低通Ha(s)%[b1,a1]=butter(N,Wc,'s');%可以代替以上三个步骤[bb,ab]=bilinear(bt,at,Fs);%双线性变换得H(z)，Fs为采样频率%求数字滤波器的幅频特性[H,w]=freqz(bb,ab,100);%对应采样间隔为0.02*pimag=abs(H);db=20*log10((mag+eps)/max(mag));subplot(2,2,3);plot(w/pi,db);axis([0,0.6,-50,10]);gridontitle('滤波器的幅频响应曲线');%用设计滤波器对输入信号进行滤波y=filter(bb,ab,x);subplot(2,2,4);n=0:55;stem(n,y,'.');axis([060-10050]);holdon;n=0:60;m=zeros(61);plot(n,m);xlabel('n');ylabel('x(n)');title('进行滤波后的心电图信号');附录2：直接设计IIR巴特沃斯数字滤波器的程序：%直接设计巴特沃斯数字滤波器Wp=0.2*pi;Ws=0.3*pi;Rp=1;Rs=40;wp=Wp/pi;ws=Ws/pi;%求数字滤波器的归一化频率[N,wc]=buttord(wp,ws,Rp,Rs);%计算巴特沃斯数字滤波器阶数N及3dB截止频率[b1,a1]=butter(N,wc);%设计巴特沃斯数字滤波器%求数字滤波器的幅频特性[H,w]=freqz(b1,a1,100);%对应采样间隔为0.02*pimag=abs(H);db=20*log10((mag+eps)/max(mag));subplot(3,1,1);plot(w/pi,db);axis([0,0.6,-50,10]);gridontitle('滤波器的幅频响应曲线');%输入信号及其波形x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,...0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,...4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];n=0:55;subplot(3,1,2);stem(n,x,'.');axis([060-10050]);holdon;n=0:60;m=zeros(61);plot(n,m);xlabel('n');ylabel('x(n)');title('心电图信号采样序列x(n)');%用设计滤波器对输入信号进行滤波y=filter(b1,a1,x);subplot(3,1,3);n=0:55;stem(n,y,'.');axis([060-10050]);holdon;n=0:60;m=zeros(61);plot(n,m);xlabel('n');ylabel('x(n)');title('进行滤波后的心电图信号');比较直接设计IIR数字滤波器和用双线性变换法设计IIR滤波器的结果是否相同？脉冲响应不变法和双线性变换法在IIR数字滤波器的设计中起到什么作用？比较他们的优缺点脉冲响应不变法的优点：

1，模拟频率到数字频率的转换时线性的；

2，数字滤波器单位脉冲响应的数字表示近似原型的模拟滤波器单位脉冲响应，因此时域特性逼近好。缺点：

会产生频谱混叠现象，只适合带限滤波器双线性变换法优点：

克服多值映射得关系，可以消除频率的混叠缺点：

时域到频域的变换是非线性的，在高频处有较大的失真。实验要求：自己编写用脉冲响应不变法设计IIR数字滤波器的程序(调用impinvar，掌握方法即可)。若把上述滤波器改为高通滤波器，试分别用脉冲响应不变法和双线性变换法设计对应的数字滤波器，观察设计结果，并说明原因。答：用脉冲响应不变法设计IIR数字滤波器附1双线性变换法设计对应的数字高通滤波器附2脉冲响应不变法不能用于设计数字高通滤波器和带阻滤波器；附1、用脉冲响应不变法设计IIR数字滤波器T=2;

%设置采样周期为2

fs=1/T;

%采样频率为周期倒数

Wp=0.25*pi/T;

Ws=0.35*pi/T;

%设置归一化通带和阻带截止频率

Ap=20*log10(1/0.9);

As=20*log10(1/0.18);

%设置通带最大和最小衰减[N,Wc]=buttord(Wp,Ws,Ap,As,'s');

%调用butter函数确定巴特沃斯滤波器阶数

[B,A]=butter(N,Wc,'s');

%调用butter函数设计巴特沃斯滤波器

W=linspace(0,pi,400*pi);

%指定一段频率值

hf=freqs(B,A,W);

%计算模拟滤波器的幅频响应

subplot(2,1,1);

plot(W/pi,abs(hf)/abs(hf(1)));

%绘出巴特沃斯模拟滤波器的幅频特性曲线grid

on;

title('巴特沃斯模拟滤波器');

xlabel('Frequency/Hz');

ylabel('Magnitude');

[D,C]=impinvar(B,A,fs);

%调用脉冲响应不变法

Hz=freqz(D,C,W);

%返回频率响应

subplot(2,1,2);

plot(W/pi,abs(Hz)/abs(Hz(1)));

%绘出巴特沃斯数字低通滤波器的幅频特性曲线grid

on;

title('巴特沃斯数字滤波器');

xlabel('Frequency/Hz');

ylabel('Magnitude');

附2、双线性变换法设计对应的数字高通滤波器Fs=10000;

Ts=1/Fs;

%采样频率和采样周期wp=3000*2*pi;ws=2000*2*pi;

%低通、高通数字滤波器截止频率Rp=3;Rs=14;

%通带波纹和阻带衰减

wpl=Fs*2*tan(ws*Ts/2);

%按频率转换公式进行转换wsl=Fs*2*tan(wp*Ts/2);

[N,Wn]=buttord(wpl,wsl,Rp,Rs,'s');

%计算模拟滤波器的最小阶数[z,p,k]=buttap(N);

%设计模拟原型低通滤波器

[b,a]=zp2tf(z,p,k);

%变零极点增益形式为系统传递函数形式[H,w]=freqs(b,a);

%求模拟低通滤波器的频率响应

figure;

plot(w,abs(H));grid;

%绘制低通频率特性title('模拟低通原型');

xlabel('w');ylabel('幅度');[e1,f1]=lp2lp(b,a,2000*2*pi)

%原型低通转换为低通模拟滤波器

[h1,w1]=freqs(e1,f1);

%求模拟低通滤波器的频率响应figu