特别说明:完整45期数字信号处理教程,原创高性能示波器代码全开源地址:链接
第39章 IIR滤波器的Matlab设计
本章节讲解IIR滤波器的Matlab设计,主要包括巴特沃斯滤波器,切比雪夫I型和II型滤波器以及椭圆滤波器。 39.1 巴特沃斯滤波器的设计 39.2 切比雪夫滤波器的设计 39.3 椭圆滤波器的设计 39.4 总结
39.1 巴特沃斯滤波器的设计
39.1.1 butter函数功能:用于设计Butterworth(巴特沃斯)滤波器 语法:[b, a] = butter(n, Wn); 说明:butter函数可以设计低通、带通、高通和带阻数字滤波器,其特性可以使通带内的幅度响应最大限度地平坦,但会损失截止频率处的下降斜度,使幅度响应衰减较慢。 l [b,a] = butter(n,Wn)可以设计截止频率为Wn的n阶低通butterworth滤波器,其中截止频率Wn应满足0Wn1,Wn=1相当于0.5fs(采样频率)。当Wn = [W1 W2]时,butter函数产生一个2n阶的数字带通滤波器,其通带为W1 < W < W2。 l [b,a] = butter(n,Wn,'ftype')可以设计高通或带阻滤波器。当ftype=high时,可设计截止频率为Wn的高通滤波器;当ftype=stop时,可设计带阻滤波器,此时Wn = [W1 W2],阻带为W1 < W < W2。 使用butter函数设计滤波器,可以使通带内的幅度响应最大地平坦,但会损失截止频率处的下降斜度。因此,butter函数主要用于设计通带平坦的数字滤波器。
39.1.2 buttord函数功能:用来选择Butterworth滤波器的阶数。 语法:[n, Wn] = buttord(WP, WS, RP, Rs); 说明:buttord函数可以在给定滤波器性能的情况下,选择Butterworth数字滤波器的最小阶数,其中WP和WS分别是通带和阻带的截止频率,其值为0Wp(或Ws)1,当该值为1时表示0.5fs(采样率)。RP和 Rs分别是通带和阻带区的波纹系数和衰减系数。 [n, Wn] = buttord(WP, WS, RP, Rs)可以得到高通、带通和带阻滤波器的最小阶数n。 当WP>WS时,为高通滤波器;当WP, WS为二元矢量时,若WP<WS,则为带通和带阻滤波器,此时Wn也为二元矢量。 利用buttord函数可得到Butterworth数字滤波器的最小阶数n,并使通带(0,WP)内的波纹系数小于RP,阻带(WS, 1)内衰减系数大于Rs。buttord函数还可以得到截止频率Wn,再利用butter函数可产生满足指定性能的滤波器。 使用butter函数设计数字滤波器,可以使通带内的幅度响应最大限度地平坦,但在截止频率附件幅度响应衰减慢。如果期望幅度响应下降斜度大,衰减快,可使用Elliptic(椭圆)或Chebyshev(切比雪夫)滤波器。
39.1.3 巴特沃斯低通滤波器设计 下面我们通过一个实例来讲解巴特沃斯低通滤波器的设计。原始信号是由50Hz正弦波和200Hz的正弦波组成,将200Hz的正弦波当做噪声滤掉,下面通过函数butter设计一组低通滤波器系数,其阶数是2,截止频率为0.25(也就是125Hz),采样率1Kbps。Matlab运行代码如下: - fs=1000; %设置采样频率 1k
- N=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('原始信号');
- grid on;
-
- subplot(222);
- y=fft(x,N); %绘制原始信号的幅频响应
- plot(f,abs(y));
- xlabel('频率/Hz');
- ylabel('振幅');
- title('原始信号FFT');
- grid on;
-
- subplot(223);
- Wc=2*125/fs; %设置截止频率125Hz
- [b,a]=butter(2,Wc); %获取2阶IIR滤波系数
- % y2=filter(b,a,x);
- y2=filtfilt(b,a,x); %计算滤波后的波形y2
- y3=fft(y2,N); %滤波后波形的幅频响应
- plot(f,abs(y3));
- xlabel('频率/Hz');
- ylabel('振幅');
- title('滤波后信号FFT');
- grid on;
-
- [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)]);
- grid on;
复制代码运算Malab结果如下: 从滤波的效果来看,2阶的IIR滤波器能够达到将近20的信噪比,比使用FIR需要更少的阶数。
39.1.4 巴特沃斯高通滤波器设计 下面我们通过一个实例来讲解巴特沃斯高通滤波器的设计。原始信号是由50Hz正弦波和200Hz的正弦波组成,将50Hz的正弦波当做噪声滤掉,下面通过函数butter设计一组高通滤波器系数,其阶数是2,截止频率为0.25(也就是125Hz),采样率1Kbps。Matlab运行代码如下: - fs=1000; %设置采样频率 1k
- N=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('原始信号');
- grid on;
-
- subplot(222);
- y=fft(x,N); %绘制原始信号的幅频响应
- plot(f,abs(y));
- xlabel('频率/Hz');
- ylabel('振幅');
- title('原始信号FFT');
- grid on;
-
- 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); %计算滤波后的波形y2
- y3=fft(y2,N); %滤波后波形的幅频响应
- plot(f,abs(y3));
- xlabel('频率/Hz');
- ylabel('振幅');
- title('滤波后信号FFT');
- grid on;
-
- [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)]);
- grid on;
复制代码Matlab运行结果如下: 从滤波的效果来看,2阶的IIR滤波器效果还是比较好的。
39.1.5 巴特沃斯带通滤波器设计 下面我们通过一个实例来讲解巴特沃斯带通滤波器的设计。原始信号是由50Hz正弦波和200Hz的正弦波组成,将50Hz的正弦波当做噪声滤掉,下面通过函数butter设计一组带通滤波器系数,其阶数是2,通带为125Hz到300Hz,采样率1Kbps。Matlab运行代码如下: - fs=1000; %设置采样频率 1k
- N=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('原始信号');
- grid on;
-
- subplot(222);
- y=fft(x,N); %绘制原始信号的幅频响应
- plot(f,abs(y));
- xlabel('频率/Hz');
- ylabel('振幅');
- title('原始信号FFT');
- grid on;
-
- subplot(223);
- Wn=[125*2 300*2]/fs; %设置通带125Hz到300Hz
- [b,a]=butter(1,Wn); %注意第一个参数虽然是1,但生成的却是2阶IIR滤波器系数
- % y2=filter(b,a,x);
- y2=filtfilt(b,a,x); %计算滤波后的波形y2
- y3=fft(y2,N); %滤波后波形的幅频响应
- plot(f,abs(y3));
- xlabel('频率/Hz');
- ylabel('振幅');
- title('滤波后信号FFT');
- grid on;
-
- [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)]);
- grid on;
复制代码Matlab的计算结果如下: 从滤波的效果来看,2阶的IIR滤波器效果还是比较好的。
39.1.6 巴特沃斯带阻滤波器设计 下面我们通过一个实例来讲解巴特沃斯带阻滤波器的设计。原始信号是由50Hz正弦波和200Hz的正弦波组成,将200Hz的正弦波当做噪声滤掉,下面通过函数butter设计一组带阻滤波器系数,其阶数是2,阻带为125Hz到300Hz,采样率1Kbps。Matlab运行代码如下: - fs=1000; %设置采样频率 1k
- N=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('原始信号');
- grid on;
-
- subplot(222);
- y=fft(x,N); %绘制原始信号的幅频响应
- plot(f,abs(y));
- xlabel('频率/Hz');
- ylabel('振幅');
- title('原始信号FFT');
- grid on;
-
- subplot(223);
- Wn=[125*2 300*2]/fs; %设置阻带125Hz到300Hz
- [b,a]=butter(1,Wn, 'stop'); %注意第一个参数虽然是1,但生成的却是2阶IIR滤波器系数
- % y2=filter(b,a,x);
- y2=filtfilt(b,a,x); %计算滤波后的波形y2
- y3=fft(y2,N); %滤波后波形的幅频响应
- plot(f,abs(y3));
- xlabel('频率/Hz');
- ylabel('振幅');
- title('滤波后信号FFT');
- grid on;
-
- [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)]);
- grid on;
复制代码Maltab运行结果如下: 从滤波的效果来看,2阶的IIR滤波器效果还是比较好的。
|
39.2.2 cheby1ord函数
39.2.3 cheby2函数
39.2.4 cheby2ord函数
39.2.5 切比雪夫I型低通滤波器设计
39.2.6 切比雪夫I型高通滤波器设计
39.2.7 切比雪夫I型带通滤波器设计
39.2.8 切比雪夫I型带阻滤波器设计
39.3.1 ellip函数
39.3.2 ellipord函数
39.3.3 椭圆滤波器设计