第39章IIR滤波器的Matlab设计

安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第11页页共共2211页页第第3399章章IIIIRR滤滤波波器器的的MMaattllaabb设设计计本章节讲解IIR滤波器的Matlab设计，主要包括巴特沃斯滤波器，切比雪夫I型和II型滤波器以及椭圆滤波器。39.1巴特沃斯滤波器的设计39.2切比雪夫滤波器的设计39.3椭圆滤波器的设计39.4总结3399..11巴巴特特沃沃斯斯滤滤波波器器的的设设计计3399..11..11bbuutttteerr函函数数功能：用于设计Butterworth（巴特沃斯）滤波器语法：[b,a]=butter(n,Wn);说明：butter函数可以设计低通、带通、高通和带阻数字滤波器，其特性可以使通带内的幅度响应最大限度地平坦，但会损失截止频率处的下降斜度，使幅度响应衰减较慢。[b,a]=butter(n,Wn)可以设计截止频率为Wn的n阶低通butterworth滤波器，其中截止频率Wn应满足0≤Wn≤1，Wn=1相当于0.5fs（采样频率）。当Wn=[W1W2]时，butter函数产生一个2n阶的数字带通滤波器，其通带为W1WW2。[b,a]=butter(n,Wn,'ftype')可以设计高通或带阻滤波器。当ftype=high时，可设计截止频率为Wn的高通滤波器；当ftype=stop时，可设计带阻滤波器，此时Wn=[W1W2]，阻带为W1WW2。使用butter函数设计滤波器，可以使通带内的幅度响应最大地平坦，但会损失截止频率处的下降斜度。因此，butter函数主要用于设计通带平坦的数字滤波器。3399..11..22bbuuttttoorrdd函函数数功能：用来选择Butterworth滤波器的阶数。语法：[n,Wn]=buttord（WP,WS,RP,Rs）；说明：buttord函数可以在给定滤波器性能的情况下，选择Butterworth数字滤波器的最小阶数，其中WP和WS分别是通带和阻带的截止频率，其值为0≤Wp（或Ws）≤1，当该值为1时表示0.5fs（采样率）。RP和Rs分别是通带和阻带区的波纹系数和衰减系数。安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第22页页共共2211页页[n,Wn]=buttord（WP,WS,RP,Rs）可以得到高通、带通和带阻滤波器的最小阶数n。当WPWS时，为高通滤波器；当WP,WS为二元矢量时，若WPWS，则为带通和带阻滤波器，此时Wn也为二元矢量。利用buttord函数可得到Butterworth数字滤波器的最小阶数n，并使通带（0，WP）内的波纹系数小于RP,阻带（WS,1）内衰减系数大于Rs。buttord函数还可以得到截止频率Wn，再利用butter函数可产生满足指定性能的滤波器。使用butter函数设计数字滤波器，可以使通带内的幅度响应最大限度地平坦，但在截止频率附件幅度响应衰减慢。如果期望幅度响应下降斜度大，衰减快，可使用Elliptic(椭圆)或Chebyshev（切比雪夫）滤波器。3399..11..33巴巴特特沃沃斯斯低低通通滤滤波波器器设设计计下面我们通过一个实例来讲解巴特沃斯低通滤波器的设计。原始信号是由50Hz正弦波和200Hz的正弦波组成，将200Hz的正弦波当做噪声滤掉，下面通过函数butter设计一组低通滤波器系数，其阶数是2，截止频率为0.25（也就是125Hz），采样率1Kbps。Matlab运行代码如下：fs=1000;%设置采样频率1kN=1024;%采样点数n=0:N-1;t=0:1/fs:1-1/fs;%时间序列f=n*fs/N;%频率序列x1=sin(2*pi*50*t);%信号x2=sin(2*pi*200*t);%噪声x=x1+x2;%信号混合subplot(221);plot(t,x);%绘制原始信号xlabel('时间');ylabel('幅值');title('原始信号');gridon;subplot(222);y=fft(x,N);%绘制原始信号的幅频响应plot(f,abs(y));xlabel('频率/Hz');安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第33页页共共2211页页ylabel('振幅');title('原始信号FFT');gridon;subplot(223);Wc=2*125/fs;%设置截止频率125Hz[b,a]=butter(2,Wc);%获取2阶IIR滤波系数%y2=filter(b,a,x);y2=filtfilt(b,a,x);%计算滤波后的波形y2y3=fft(y2,N);%滤波后波形的幅频响应plot(f,abs(y3));xlabel('频率/Hz');ylabel('振幅');title('滤波后信号FFT');gridon;[H,F]=freqz(b,a,512);subplot(224);plot(F/pi,abs(H));xlabel('归一化频率');%绘制绝对幅频响应ylabel('幅度');Ps=sum(x1.^2);%信号的总功率Pu=sum((y2-x1).^2);%剩余噪声的功率SNR=10*log10(Ps/Pu);%信噪比title(['Order=',int2str(2),'SNR=',num2str(SNR)]);gridon;运算Malab结果如下：安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第44页页共共2211页页从滤波的效果来看，2阶的IIR滤波器能够达到将近20的信噪比，比使用FIR需要更少的阶数。3399..11..44巴巴特特沃沃斯斯高高通通滤滤波波器器设设计计下面我们通过一个实例来讲解巴特沃斯高通滤波器的设计。原始信号是由50Hz正弦波和200Hz的正弦波组成，将50Hz的正弦波当做噪声滤掉，下面通过函数butter设计一组高通滤波器系数，其阶数是2，截止频率为0.25（也就是125Hz），采样率1Kbps。Matlab运行代码如下：fs=1000;%设置采样频率1kN=1024;%采样点数n=0:N-1;t=0:1/fs:1-1/fs;%时间序列f=n*fs/N;%频率序列x1=sin(2*pi*50*t);%噪声x2=sin(2*pi*200*t);%信号x=x1+x2;%信号混合subplot(221);plot(t,x);%绘制原始信号00.51-2-1012时间幅值原始信号050010000200400600频率/Hz振幅原始信号FFT050010000200400600频率/Hz振幅滤波后信号FFT00.5100.51归一化频率幅度Order=2SNR=19.3284安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第55页页共共2211页页xlabel('时间');ylabel('幅值');title('原始信号');gridon;subplot(222);y=fft(x,N);%绘制原始信号的幅频响应plot(f,abs(y));xlabel('频率/Hz');ylabel('振幅');title('原始信号FFT');gridon;subplot(223);Wc=2*125/fs;%设置截止频率125Hz[b,a]=butter(2,Wc,'high');%获取2阶IIR滤波系数%y2=filter(b,a,x);y2=filtfilt(b,a,x);%计算滤波后的波形y2y3=fft(y2,N);%滤波后波形的幅频响应plot(f,abs(y3));xlabel('频率/Hz');ylabel('振幅');title('滤波后信号FFT');gridon;[H,F]=freqz(b,a,512);subplot(224);plot(F/pi,abs(H));xlabel('归一化频率');%绘制绝对幅频响应ylabel('幅度');Ps=sum(x2.^2);%信号的总功率Pu=sum((y2-x2).^2);%剩余噪声的功率SNR=10*log10(Ps/Pu);%信噪比title(['Order=',int2str(2),'SNR=',num2str(SNR)]);安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第66页页共共2211页页gridon;Matlab运行结果如下：从滤波的效果来看，2阶的IIR滤波器效果还是比较好的。3399..11..55巴巴特特沃沃斯斯带带通通滤滤波波器器设设计计下面我们通过一个实例来讲解巴特沃斯带通滤波器的设计。原始信号是由50Hz正弦波和200Hz的正弦波组成，将50Hz的正弦波当做噪声滤掉，下面通过函数butter设计一组带通滤波器系数，其阶数是2，通带为125Hz到300Hz，采样率1Kbps。Matlab运行代码如下：fs=1000;%设置采样频率1kN=1024;%采样点数n=0:N-1;t=0:1/fs:1-1/fs;%时间序列f=n*fs/N;%频率序列x1=sin(2*pi*50*t);%噪声x2=sin(2*pi*200*t);%信号x=x1+x2;%信号混合00.51-2-1012时间幅值原始信号050010000200400600频率/Hz振幅原始信号FFT050010000200400600频率/Hz振幅滤波后信号FFT00.5100.51归一化频率幅度Order=2SNR=19.2927安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--VV55开开发发板板系系统统篇篇手手册册22001155年年0011月月1155日日版版本本：：11..00第第77页页共共2211页页subplot(221);plot(t,x);%绘制原始信号xlabel('时间');ylabel('幅值');title('原始信号');gridon;subplot(222);y=fft(x,N);%绘制原始信号的幅频响应plot(f,abs(y));xlabel('频率/Hz');ylabel('振幅');title('原始信号FFT');gridon;subplot(223);Wn=[125*2300*2]/fs;%设置通带125Hz到300Hz[b,a]=butter(1,Wn);%注意第一个参数虽然是1，但生成的却是2阶IIR滤波器系数%y2=filter(b,a,x);y2=filtfilt(b,a,x);%计算滤波后的波形y2y3=fft(y2,N);%滤波后波形的幅频响应plot(f,abs(y3));xlabel('频率/Hz');ylabel('振幅');title('滤波后信号FFT');gridon;[H,F]=freqz(b,a,512);subplot(224);plot(F/pi,abs(H));xlabel('归一化频率');%绘制绝对幅频响应ylabel('幅度');Ps=sum(x2.^2);%信号的总功率Pu=sum((y2-x2).^2);%剩余噪声的功率安安富富莱莱UUMM440033DDSSPP教教程程SSTTMM3322--