IIR数字滤波器设计报告

实验五IIR数字滤波器设计报告实验目的1.掌握IIR数字滤波器设计方法；2.掌握利用数字滤波器处理连续信号的方法。实验内容1.熟悉用双线性变换法设计IIR数字滤波器的原理与方法。2.人体心电信号的主要频率范围为0．05～100Hz，设计带通滤波器，滤除含噪的心电信号。3.通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。具体实验步骤及实验结果此次试验，我运用了两种办法设计IIR数字滤波器。一种是直接运用MATLAB中已有的函数直接带入参数设计。一种是运用比较复杂的双线性法设计。并且分别进行了滤波，和对心电信号的滤波。函数法：运用MATLAB中已有的函数，buttord，butter，freqz，filter等，直接代入滤波器特性的参数Wp，Ws，Rp，As就可以了。程序如下：clearcloseallclcwp=100;%通带频率ws=300;%阻带频率rp=1;as=50;n=0.0005;E=0.1;f1=20;f2=400;%1000开始混叠T=0.0005;fs=1/T;[N,wn]=buttord(wp/1000,ws/1000,rp,as);[b,a]=butter(N,wn);freqz(b,a,fs,2000);[hz,w]=freqz(b,a,fs,2000);%加入模拟信号t=0:n:E-n;x=cos(2*pi*t*f1)+cos(2*pi*t*f2);y=filter(b,a,x);h=fft(x)/(E/n);i=fft(y)/(E/n);figure(2);subplot(232);stem(t/E*(1/n),abs(h));axis([0,0.5/n,-1,1]);xlabel('输入频谱(Hz)');subplot(235);stem(t/E*(1/n),abs(i),'r');xlabel('输出频谱(Hz)');axis([0,0.5/n,-1,1]);subplot(231);plot(t,x);xlabel('原始输入(s)');subplot(234);plot(t,y);xlabel('滤波之后(s)');subplot(233);plot(w,abs(hz));grid;xlabel('滤波器特性(Hz)');01002003004005006007008009001000-600-400-2000Frequency(Hz)Phase(degrees)01002003004005006007008009001000-600-400-2000200Frequency(Hz)Magnitude(dB)0100200300400500600700800900100000.20.40.60.811.21.4滤波器特性(Hz)05001000-1-0.500.51输入频谱(Hz)05001000-1-0.500.51输出频谱(Hz)00.050.1-2-1012原始输入(s)00.050.1-2-1012滤波之后(s)我的信号x=cos(2*pi*t*f1)+cos(2*pi*t*f2);为f1=20Hz，其中带有f2=400Hz的噪声，所以我设计了一个低通滤波器滤除信号中的高频噪声。截止频率为122.73Hz。滤波器阶数为6。通带频率wp=100;阻带频率ws=300;rp=1;as=50;（为了计算旁边，我选用2000Hz的采样频率，这样就可以比较容易计算pi和频率的换算。）在图一可以看到滤波器的dB表示形式，图二为滤波器的幅频特性。可以看到在噪声频率400Hz处的响应已经几乎为零了。所以足够滤除噪声。图三中第一个是信号在时域的图像，可以看到噪声比较明显。第二个图是对带噪声信号的频谱分析。可以看到两条谱线，一条在20Hz位置，噪声谱线在400Hz位置。而第三个是经过滤波的信号在时域的图像，已经把噪声滤除。第四个是经过滤波后的频谱。噪声谱线已经消失，不过在信号谱线位置发现有轻微的泄漏现象。滤波认为基本成功。双线性法：先对滤波器进行模拟设计，再用双线性的办法转换成数字形式。程序如下：clear;closeall;clc;wp=0.3*pi;%通带频率ws=0.5*pi;%阻带频率rp=1;as=30;n=0.0005;E=0.1;f1=20;f2=600;%1000开始混叠T=0.0005;fs=1/T;P=(2/T)*tan(wp/2);S=(2/T)*tan(ws/2);[cs,ds]=buttord(P,S,rp,as,'s');[z,p,k]=buttap(cs);[bt,at]=zp2tf(z,p,k);[b,a]=lp2lp(bt,at,ds);[b,a]=bilinear(b,a,fs);[hz,w]=freqz(b,a,fs,2);%加入模拟信号t=0:n:E-n;x=cos(2*pi*t*f1)+cos(2*pi*t*f2);y=filter(b,a,x);h=fft(x)/(E/n);i=fft(y)/(E/n);figure(2);subplot(222);stem(t/E*(1/n),abs(h));axis([0,0.5/n,-1,1]);xlabel('输入频谱(Hz)');subplot(224);stem(t/E*(1/n),abs(i));xlabel('输出频谱(Hz)');axis([0,0.5/n,-1,1]);subplot(221);plot(t,x);xlabel('原始输入(s)');subplot(223);plot(t,y);xlabel('滤波之后(s)');figure(1);plot(w,abs(hz));grid;xlabel('滤波器特性(pi)');00.10.20.30.40.50.60.70.80.9100.20.40.60.811.21.4滤波器特性(pi)05001000-1-0.500.51输入频谱(Hz)05001000-1-0.500.51输出频谱(Hz)00.050.1-2-1012原始输入(s)00.050.1-2-1012滤波之后(s)输入信号x=cos(2*pi*t*f1)+cos(2*pi*t*f2)为20Hz，噪声为600Hz。所以设计了通带频率wp=0.3*pi阻带频率ws=0.5*pirp=1;as=30的低通滤波器。他的幅频响应可从一图看出。二图中第一个为带噪声信号的时域图形，噪声明显。第二个是它的频谱分析，可以看到信号频谱线和噪声频谱线。第三个是经过滤波的信号图形，已经没有了噪声，比较圆滑。第四个是滤波后的频谱，噪声频谱线已经消失。滤除心电信号：把四个周期的心电信号加噪后进行滤波。程序如下：clearcloseallclcwp=450;%通带频率ws=600;%阻带频率rp=1;as=50;E=0.2;T=0.0005;fs=1/T;[N,wn]=buttord(wp/1000,ws/1000,rp,as);[b,a]=butter(N,wn);freqz(b,a,fs,2000);[hz,w]=freqz(b,a,fs,2000);%加入心电信号t=0:T:E-T;xdt=[00.0014124290.0028248590.0042372880.0056497180.0070621470.0084745760.0556900730.1029055690.1501210650.1973365620.2445520580.2917675540.3389830510.2976694920.2563559320.2150423730.1737288140.1324152540.0911016950.0497881360.0084745760.0075329570.0065913370.0056497180.0047080980.0037664780.0028248590.0018832390.000941620-0.144067797-0.288135593-0.0305084750.2271186440.4847457630.74237288110.632203390.26440678-0.103389831-0.471186441-0.838983051-0.41949152500.0008920610.0017841210.0026761820.0035682430.0044603030.0053523640.0062444250.0071364850.0080285460.0089206070.0098126670.0107047280.0115967890.0124888490.013380910.0142729710.0151650310.0160570920.0169491530.0399515740.0629539950.0859564160.1089588380.1319612590.154963680.1779661020.154963680.1319612590.1089588380.0859564160.0629539950.0399515740.0169491530.0127118640.0084745760.0042372880000000000000000000];x1=[xdt,xdt,xdt,xdt];x=x1+cos(1400*pi*t);y=filter(b,a,x);h=fft(x)/(E/T);i=fft(y)/(E/T);figure(2);subplot(223);stem(t/E*(1/T),abs(h));xlabel('输入频谱(Hz)');axis([0,0.5/T,0,0.06]);subplot(224);stem(t/E*(1/T),abs(i),'r');xlabel('输出频谱(Hz)');axis([0,0.5/T,0,0.06]);subplot(222);plot(t,x);xlabel('原始输入(s)')subplot(221);plot(w,abs(hz));grid;xlabel('滤波器特性(Hz)');figure(1);plot(t,x1,'r');holdon;plot(t,y);xlabel('滤波之后(s)');0500100000.020.040.06输入频谱(Hz)0500100000.020.040.06输出频谱(Hz)00.050.10.150.2-2-1012原始输入(s)0500100000.511.5滤波器特性(Hz)00.020.040.060.080.10.120.140.160.180.2-1-0.8-0.6-0.4-0.200.20.40.60.81滤波之后(s)首先，对新电信后进行频谱分析，发现心电信号主要频率成分都为低频。噪声频谱在信号的主要频谱之外。所以设计了一个低通滤波器，足以滤除高频噪声。由于，我设计的滤波器要求质量较高，所以为40阶，所以发生相异。不过总体来说还是符合要求的。结论：经过此次试验，对于IIR滤波器的了解可谓是深入了一大步。也了解了他和FIR滤波器本质的区别。甚至可以应用它去做出设计，真正的滤波。比较两种方法，函数方法无疑简单不少，但是双线性法又让我们更清楚的了解了双线性法的本质，可以很好的对其市值进行研究。巴特沃斯滤波器可以较好的完成对于心电信号的过滤。