【MATLAB】数字滤波器的设计
一、引言
在信号处理过程中，所处理的信号往往混有噪声，从接收到的信号中消除或减弱噪声是信号传输和处理中十分重要的问题。根据有用信号和噪声的不同特性，提取有用信号的过程称为滤波,实现滤波功能的系统称为滤波器。在以往的模拟电路中用的都是模拟滤波器，在近代电信设备和各类控制系统中，由于数字化的普及，数字滤波器已得到了广泛的应用。数字滤波器与传统模拟滤波器在实现方式上存在很大的差异。传统的模拟滤波器主要是硬件实现,它的硬件部分主要包括电容,电感和电阻等元件,而数字滤波器在硬件实现上主要涉及A/D转换器、D/A转换器,、寄存器、存储器及微处理器等。数字滤波器的另一特点是可以用软件实现,即通过编程用算法来实现。数字滤波器与模拟滤波器相比,有其独特的优点，比如体积小、成本低、参数容易调整、有较高的精度、工作效率高等，但它们也有共同之处,例如滤波器的选频特性,即都用频率响应作为滤波器的主要技术指标。
二、数字滤波器
1.数字滤波器的分类
数字滤波器有多种分类方法，每一种方法都从不同的侧面揭示了数字滤波器的特性，主要有按频率分布特性的分类方法、按实现方式的分类方法和按脉冲响应特性的分类方法。
（1）按频率分布特性分类
按频率分布特性分类方法与传统的模拟滤波器分类方法一致,将数字滤波器分为低通、高通、带通和带阻四类,如图所示。这种分类方法从频域的角度阐述了滤波器对频率的选择特性。在实际应用中,对于频率分布可分的信号,可选择性通过，以达到滤波的目的。
（2）按实现方式分类
按实现方式的不同，数字滤波器可以分为非递归型(Nonrecursive)数字滤波器和递归型(Recursive)数字滤波器两大类。
（3）按脉冲响应特性分类
按数字滤波器对脉冲响应的特性来分类,数字滤波器可以分为有限脉冲响应FIR(FiniteImpulse Response)数字滤波器和无限脉冲响应 IIR(Infinite Impulse Response)数字滤波器FIR滤波器的传递函数只有零点,不含极点,它的单位脉冲响应h(k)只包含有限个非零值,即这种数字滤波器的脉冲响应是时间有限的。
IIR滤波器既有零点又有极点,它的脉冲响应h(k)中含有无限多个非零值,即这种滤波器的脉冲响应是无限长时间序列,在一定的时间后可能变小,但不会为零。
一般情况下,FIR滤波器比较适合用非递归型方式来实现,而IIR滤波器比较适合用递归型方式来实现。但无论是FIR滤波器还是IIR滤波器,都可用两种方法中的任一种来实现,只是按上述方法更简便而已。
2.典型模拟低通滤波器
介绍一些典型的模拟滤波器,主要说明数字滤波器的设计和应用。为什么要介绍模拟滤波器呢?原因是有些数字滤波器是以模拟滤波器为原型,把它们转换成数字滤波器的。在本节中介绍模拟滤波器并不是要说明模拟滤波器如何设计和实现,而是要说明模拟滤波器的一些特性,为以后采取某些方法转换成数字滤波器打下基础。典型模拟低通滤波器有巴特沃斯(Butterworth)、切比雪夫(Chebyshev)I型、切比雪夫Ⅱ型和椭圆型(Elliptic)滤波器等。
3.模拟原型低通滤波器的频率变换
设计模拟滤波器时,已知滤波器的具体指标后常常先设计模拟原型低通滤波器,通过适当的频率变换可转换成满足滤波器具体指标的各类滤波器(低通、高通、带通和带阻)。
4.模拟滤波器设计的 MATLAB 函数
在这里讨论模拟滤波器并不是要构成一个模拟滤波器，而是在模拟滤波器的基础上通过某种变换，转换成数字滤波器，使相应的数字滤波器也有模拟滤波器的某些特性。
在MATLAB中常用的模拟滤波器设计数如下：
模拟滤波器阶数选择函数
(1)巴特沃斯模拟滤波器阶数的选择
函数:buttord;功能:计算巴特沃斯模拟滤波器的阶数;
调用格式:[n,Wn】= buttord(Wp,Ws,Rp,Rs ,'s')(2)切比雪夫工型模拟滤波器阶数的选择
函数:cheblord;功能:计算切比雪夫工型模拟滤波器的阶数
调用格式;[n,Wn]= cheblord( p,Ws,Rp, Rs ,'s')(3)切比雪夫Ⅱ型模拟滤波器阶数的选择
函数:cheb2ord;功能:计算切比雪夫Ⅱ型模拟滤波器的阶数
调用格式:[n,Wn]= cheb2ord(wp,Ws ,Rp,Rs , 's')(4)椭圆型模拟滤波器阶数的选择
函数:ellipord;功能:计算椭圆型模拟滤波器的阶数;
调用格式:[n,Wn]= ellipord(wp,Ws,Rp,Rs,'s')在MATLAB中设计模拟滤波器的函数
(1)巴特沃斯模拟滤波器的设计
函数:butter功能:巴特沃斯模拟滤波器的设计
调用格式:[b,a]= butter(n,Wn,'s');[b,a]= butter(n, Wn,'ftype' ,'s');(2)切比雪夫工型模拟滤波器的设计
函数:chebyl功能:切比雪夫工型模拟滤波器的设计(通带等波纹)
调用格式:[b,a]= cheby1(n,Rp, Wn, 's');[b,a]= cheby1(n,Rp,Wn,'ftype' ,'s');(3)切比雪夫Ⅱ型模拟滤波器的设计函数:cheby2功能:切比雪夫Ⅱ型模拟滤波器的设计(阻带等波纹)
调用格式:[b,a]= cheby2(n,Rs , Wn,'s')[b,a]= cheby2(n,Rs,Wn, 'ftype', 's')(4)椭圆型模拟滤波器的设计
函数:ellip功能:椭圆型模拟滤波器的设计
调用格式:[b,a]= ellip(n,Rp,Rs, Wn,'s');[b,a]= ellip(n,Rp,Rs, Wn,'ftype','s')以上4个函数中都带有's',表示设计模拟滤波器,输人参数中n表示为低通模拟滤波器的阶次,Wn表示为截止频率,无论高通、带通和带阻滤波器,在设计中最终都等效于一个低通滤波器,而等效低通滤波器的n和Wn可由以上求阶次的函数计算得到。Rp表示通带内的波纹(单位为dB),Rs表示阻带的衰减(单位为dB)。以上的函数并不是都用到这4个参数,有的用2个参数,有的用3个参数或4个参数。当Wn=[W1 W2](WI<W2)时,表示设计一个带通滤波器,函数将产生一个2n阶的数字带通滤波器,其通带频率为W1<w<W2。
当带有参数'ftype'时,表示可设计出高通或带阻滤波器:
>当ftype=high时,设计出截止频率为Wn的高通滤波器。
>当ftype=stop时,设计出带阻滤波器,这时Wn=[W1 W2],且阻带频率为 W1<w<W2.
模拟低通滤波器原型设计函数
(1)巴特沃斯滤波器原型
函数:buttap功能:巴特沃斯模拟低通滤波器原型
调用格式:[z,pk]= buttap(n)(2)切比雪夫工型滤波器原型函数:cheblap功能:切比雪夫工型模拟低通滤波器原型
调用格式:[z,p,k]= cheblap(n,Rp)(3)切比雪夫Ⅱ型滤波器原型
函数:cheb2ap功能:切比雪夫Ⅱ型模拟低通滤波器原型
调用格式:[z,p,k= cheb2ap(n,Rs)(4)椭圆型滤波器原型
函数:ellipap功能:椭圆型模拟低通滤波器原型
调用格式:[z,p,k=ellipap(n,Rp,Rs)说明:以上4个函数中,n是由模拟滤波器阶次函数求出的阶次;Rp、RS分别是通带内的波纹系数和阻带内的波纹系数,单位都是dB。
频率变换函数
(1)低通到带通变换
函数:lp2bp功能:低通到带通模拟滤波器变换
调用格式:[bt ,at]= lp2bp(b,a,wo,Bw)[At,Bt,Ct,Dt]= lp2bp(A,B,C,D,Wo,Bw)(2)低通到高通变换
函数:1p2hp功能:低通到高通模拟滤波器变换
调用格式:[bt,at]= lp2hp(b,a,wo)At,Bt,Ct,Dt]=lp2hp(A,B,C,D,Wo)(3)低通到带阻变换函数:1p2bs功能:低通到带阻模拟滤波器变换
调用格式:[bt,at]=lp2bs(b,a,wo,Bw)At,Bt,Ct,Dt]=lp2bs(A,B,C,D, Wo,Bw)(4)低通到低通变换
函数:1p2lp功能:低通到低通模拟滤波器变换
调用格式:[bt,at]=lp2lp(b,a,wo)At,Bt,Ct,Dt]=1p21p(A,B,C,D,Wo)滤波器频率分析函数
函数:freqs功能:已知模拟滤波器系数b和a,求出模拟滤波器的幅值响应和相角响应
调用格式:[H,w]= freqs(b,a);[H,w]= freqs(b,a,n);H=freqs(b,a,w)说明:在调用格式[H,w]=freqs(b,a)中,freqs函数自定义计算200个模拟频率,响应曲线在H中,对应的模拟角频率在矢量w中。参数b和a是模拟滤波器的系数,对应的传递函数如式(3-2-17)所示。当要修改模拟频率的个数时,可以设置参数n,调用格式[H,w]=fregs(b,a,n)。也可以设置模拟角频率矢量w,调用格式 H=freqs(b,a,w)。
巴特沃斯、切比雪夫Ⅰ型、Ⅱ型和椭圆型滤波器的相同和不同之处
clear ; close all; clc;
wp=[0.2*pi 0.3*pi]; % 设置通带频率
ws=[0.1*pi 0.4*pi]; % 设置阻带频率
Rp=1; Rs=20; % 设置波纹系数
% 巴特沃斯滤波器设计
[N,Wn]=buttord(wp,ws,Rp,Rs,'s'); % 求巴特沃斯滤波器阶数
fprintf('巴特沃斯滤波器 N=%4d\n',N) % 显示滤波器阶数
[bb,ab]=butter(N,Wn,'s'); % 求巴特沃斯滤波器系数
W=0:0.01:2; % 设置模拟频率
[Hb,wb]=freqs(bb,ab,W); % 求巴特沃斯滤波器频率响应
plot(wb/pi,20*log10(abs(Hb)),'b')% 作图
hold on
% 切比雪夫I型滤波器设计
[N,Wn]=cheb1ord(wp,ws,Rp,Rs,'s'); % 求切比雪夫I型滤波器阶数
fprintf('切比雪夫I型滤波器 N=%4d\n',N) % 显示滤波器阶数
[bc1,ac1]=cheby1(N,Rp,Wn,'s'); % 求切比雪夫I型滤波器系数
[Hc1,wc1]=freqs(bc1,ac1,W); % 求切比雪夫I型滤波器频率响应
plot(wc1/pi,20*log10(abs(Hc1)),'k')% 作图
% 切比雪夫II型滤波器设计
[N,Wn]=cheb2ord(wp,ws,Rp,Rs,'s'); % 求切比雪夫II型滤波器阶数
fprintf('切比雪夫II型滤波器 N=%4d\n',N) % 显示滤波器阶数
[bc2,ac2]=cheby2(N,Rs,Wn,'s'); % 求切比雪夫II型滤波器系数
[Hc2,wc2]=freqs(bc2,ac2,W); % 求切比雪夫II型滤波器频率响应
plot(wc2/pi,20*log10(abs(Hc2)),'r')% 作图
% 椭圆型滤波器设计
[N,Wn]=ellipord(wp,ws,Rp,Rs,'s'); % 求椭圆型滤波器阶数
fprintf('椭圆型滤波器 N=%4d\n',N) % 显示滤波器阶数
[be,ae]=ellip(N,Rp,Rs,Wn,'s'); % 求椭圆型滤波器系数
[He,we]=freqs(be,ae,W); % 求椭圆型滤波器频率响应
% 作图
plot(we/pi,20*log10(abs(He)),'g')
line([0 max(we/pi)],[-20 -20],'color','k','linestyle','--');
line([0 max(we/pi)],[-1 -1],'color','k','linestyle','--');
line([0.2 0.2],[-30 2],'color','k','linestyle','--');
line([0.3 0.3],[-30 2],'color','k','linestyle','--');
axis([0 max(we/pi) -30 2]); %grid;
legend('巴特沃斯滤波器','切比雪夫I型滤波器','切比雪夫II型滤波器','椭圆型滤波器')
xlabel('角频率{\omega}/{\pi}'); ylabel('幅值/dB')
set(gcf,'color','w'
在频带变换的模拟滤波器设计中,怎样计算Wn和Bs
设计了低通滤波器的原型后,按设计要求转换到带通(或带阻滤波器),频带转换函数有lp2bp(ba,as,Wo,Bs)或lp2bs(ba,as,Wo,Bs),其中Wo、Bs为Wo=sqrt(Wl* W2),Bs=W2-W1。那么W1是W2通带频率,还是阻带频率,还是其他频率?我们还是以一个实例来说明
设计一个巴特沃斯模拟带通滤波器,用频带转换方法编程。带通值为 wp1=0.2π,wp2=0.3π,带阻值为wsl=0.1π,ws2=0.4π,Rp=1,Rs=20。程序清单如下:
clear all; close all; clc;
wp=[0.2*pi 0.3*pi]; % 设置通带频率
ws=[0.1*pi 0.4*pi]; % 设置阻带频率
Rp=1; Rs=20; % 设置波纹系数
% 巴特沃斯滤波器设计
[N,Wn]=buttord(wp,ws,Rp,Rs,'s'); % 求巴特沃斯滤波器阶数
fprintf('巴特沃斯滤波器 N=%4d\n',N) % 显示滤波器阶数
[bn,an]=butter(N,Wn,'s'); % 求巴特沃斯滤波器系数
W1=Wn(1); W2=Wn(2); % 设置W1,W2
Wo=sqrt(W1*W2);
Bw=W2-W1;
[z,p,k]=buttap(N); % 设计低通原型数字滤波器
[Bap,Aap]=zp2tf(z,p,k); % 零点极点增益形式转换为传递函数形式
[bb,ab]=lp2bp(Bap,Aap,Wo,Bw); % 低通滤波器频率转换
W=0:0.01:2; % 设置模拟频率
[Hn,wn]=freqs(bn,an,W); % 求巴特沃斯滤波器频率响应
[Hb,wb]=freqs(bb,ab,W); % 求巴特沃斯滤波器频率响应
% 作图
plot(wn/pi,20*log10(abs(Hn)),'r','linewidth',2)
hold on
plot(wb/pi,20*log10(abs(Hb)),'k')% 作图
line([0 max(wb/pi)],[-20 -20],'color','k','linestyle','--');
line([0 max(wb/pi)],[-1 -1],'color','k','linestyle','--');
line([0.2 0.2],[-30 2],'color','k','linestyle','--');
line([0.3 0.3],[-30 2],'color','k','linestyle','--');
title('两种编程方法设计巴特沃斯带通滤波器');
xlabel('圆频率{\omega}/{\pi}'); ylabel('幅值/dB')
set(gcf,'color','w');
axis([0 max(wb/pi) -30 2]);
legend('第1种编程方法','第2种编程方法')
4.IIR 数字滤波器设计的MATLAB 函数
在IIR数字滤波器设计中有直接的设计函数,即在已知数字滤波器的指标后调用函数直接设计得到滤波器的系数,或者先设计模拟滤波器,再通过转换得到数字滤波器的系数。以下分别给出它们的MATLAB函数。不论用哪一种方法设计的数字滤波器都以巴特沃斯、切比雪夫I型、切比雪夫Ⅱ型和椭圆型模拟滤波器为原型,并把它们映射为数字滤波器。
IIR 数字滤波器阶数选择函数
(1)巴特沃斯数字滤波器阶数的选择函数:buttord功能:计算巴特沃斯数字滤波器的阶数
调用格式:[n,Wn]=buttord(Wp,Ws,Rp,Rs)(2)切比雪夫工型数字滤波器阶数的选择函数:cheblord功能:计算切比雪夫工型数字滤波器的阶数
调用格式:[n,Wn]= cheblord(Wp,Ws,Rp,Rs)(3)切比雪夫Ⅱ型数字滤波器阶数的选择函数:cheb2ord功能:计算切比雪夫Ⅱ型数字滤波器的阶数
调用格式:[n,Wn]=cheb2ord(wp,Ws,Rp,Rs)(4)椭圆型数字滤波器阶数的选择函数:ellipord功能:计算椭圆型数字滤波器的阶数
调用格式:[n,Wn]=ellipord(Wp,Ws,Rp,Rs)说明:以上4个函数与在模拟滤波器设计中一样,只是在输入参数中不再带有's。输出参数中的n是求出滤波器最小的阶数,Wn是等效低通滤波器的截止频率;Wp和Ws分别是通带和阻带的频率(截止频率),是归一化的角频率,单位为rad/pi,数值在0~1之间。当Wp>Ws时,为高通滤波器;当Wp和Ws为二元矢量时,为带通或带阻滤波器,这时求出的 Wn也为二元矢量。
模拟滤波器数字化
(1)脉冲响应不变法
函数:impinvar功能:脉冲响应不变法实现模拟到数字的滤波器变换
调用格式:[bz ,az]=impinvar(b,a,Fs);[bz,az]= impinvar(b,a)说明:b和a是模拟滤波器系数,bz和az为数字滤波器的系数,两者通过脉冲不变法的变换,把模拟滤波器的脉冲响应按Fs采样后等同于数字滤波器的脉冲响应。而在使用impinvar(b,a)时其中 Fs 默认为 1 Hz。
(2)双线性Z变换法
函数:bilinear功能:双线性Z变换法实现模拟到数字的滤波器变换
调用格式:[bz,az]=bilinear(b,a,Fs);[zd,pd,kd]=bilinear(z,p,k,Fs);说明:b和a是模拟滤波器系数,bz和az为数字滤波器的系数,两者通过双线性Z变换法把模拟滤波器的系数变换为数字滤波器的系数。z、p、k为模拟滤波器的零极点和增益,zd、pd、kd为数字滤波器的零极点和增益,Fs为采样频率.
直接设计数字滤波器的函数
(1)巴特沃斯模拟滤波器的设计
函数:butter功能:巴特沃斯模拟滤波器设计
调用格式:[b,a]= butter(n,Wn);[b,a]= butter(n,Wn,'ftype')(2)切比雪夫工型模拟滤波器的设计
函数:cheby1功能:切比雪夫工型模拟滤波器设计(通带等波纹)
调用格式:[b,a]= cheby1(n,Rp,Wn);[b,a]= cheby1(n, Rp, Wn, 'ftype')(3)切比雪夫Ⅱ型模拟滤波器的设计
函数:cheby2功能:切比雪夫Ⅱ型模拟滤波器设计(阻带等波纹)
调用格式:[b,a]= cheby2(n,Rs,Wn);[b,a]= cheby2(n,Rs,Wn, 'ftype');(4)椭圆型模拟滤波器的设计
函数:ellip功能:椭圆型模拟滤波器设计
调用格式:[b,a]= ellip(n,Rp,Rs ,Wn);[b,a]= ellip(n,Rp,Rs,Wn,'ftype');说明:以上4个函数与在模拟滤波器设计中一样,只是在输人参数中不再带有's'。输人参数中,n表示低通模拟滤波器的阶次,Wn表示截止频率,无论高通、带通或带阻滤波器,在设计中最终都等效于一个低通滤波器,而等效低通滤波器的n和Wn可由“IIR数字滤波器阶次选择的函数”中的有关函数计算得到。Rp表示通带内的波纹(单位为dB),Rs表示为阻带的衰减(单位为dB)。以上函数并不是都用到这4个参数,有的用2个参数,有的用3个参数或4个参数。当Wn=[W1 W2](WI<W2)时,表示设计一个带通滤波器,函数将产生一个2n阶的数字带通滤波器,W1 W2
滤波器特性分析函数
(1)数字滤波器频率特性分析1
函数:freqz功能:已知数滤波器系数b和a,求出数字滤波器的幅值响应和相角响应
1️⃣格式一：调用freqz函数，返回幅度响应H和角频率w。
调用格式：[H,w]=freqz(b,a)
2️⃣格式二：调用freqz函数，返回幅度响应H、角频率w和采样点数n。
调用格式：[H,w]=freqz(b,a,n)
3️⃣格式三：调用freqz函数，返回幅度响应H、频率f和采样点数n。其中采样频率为Fs。
调用格式：[H,f]=freqz(b,a,n,Fs)
4️⃣格式四：调用freqz函数，返回幅度响应H和整个角频率范围的数据w。
调用格式：[H,w]=freqz(b,a,n,'whole')
5️⃣格式五：调用freqz函数，返回幅度响应H、整个频率范围f和采样点数n。其中采样频率为Fs。
调用格式：[H,f ]=freqz(b,a,n,'whole',Fs)
👉Bonus：你也可以使用freqs函数来计算幅度响应H，采用的是角频率w。调用格式：H=freqs(b,a,w)
👉别忘了，如果没有输入参数，调用format函数会返回默认的幅度响应说明:在调用格式[H,w]=freqz(b,a)中,freqz函数自定义计算512个数字频率,响应曲线在H中,对应的数字角频率在矢量w中,w在0~x区间内,表示归一化角频率。参数b和a是数字滤波器的系数。设置参数n后,即表示频率矢量和H都有n个元素;当带有采样频率Fs后,输出频率矢量1,将在0~Fs/2的区间内。设置参数'whole后频率将在0~2-区间内或在0~Fs的区间内。当不带任何输出参数时,freqz(b,a)将直接显示出该滤波器的幅频响应曲线与相频响应曲线。
(2)数字滤波器频率特性分析2
函数:fvtool功能:显示出数字滤波器的响应曲线
调用格式:fvtool(b,a)
fvtool(Hd)
fvtool(b,a,bz,az,...b.,an)
fvtool(Hd,,Hd...)说明:输入参数b和a是数字滤波器系数;b1,a1,b2,a2,…,bn,an是多个滤波器的系数要求多个滤波器的响应曲线在一张图上画出);Hd是滤波器系数集合;Hd1,Hd2,…是多个滤波器的系数集合。从调用格式可以看出,多个滤波器的响应曲线通过本函数一次就能同时画在一张图上。
(3)群延迟的频率分析
函数:grpdelay功能:已知数字滤波器系数b和a,求出数字滤波器群延迟随频率变化的曲线
调用格式:[gd,w]= grpdelay(b,a)
[gd,w]=grpdelay(b,a,n)
[gd,f]= grpdelay(b,a,n,Fs)
[gd,w]= grpdelay(b,a,n,'whole')
[gd,f]= grpdelay(b,a,n,'whole',Fs)
gd= grpdelay(b,a,w)
grpdelay(b.a)说明:grpdelay函数类似于freqz函数,自定义计算512个数字频率,设置n后就有n个数字频率。延迟量曲线在gd中,对应的数字角频率在矢量w中,w在0~π区间内,表示归一化角频率。w和gd的长度或为512个(默认值)或为n个。
(4)产生数字滤波器的脉冲响应
函数:impz功能:已知数字滤波器系数b和a,求出数字滤波器的脉冲响应
调用格式:[h,t]= impz(b,a);[h,t]=impz(b,a,n);[h,t]=impz(b,a,n,Fs);impz(b,a,n,Fs)说明:b和a为数字滤波器系数;n为脉冲响应样点总数;Fs为采样频率,默认值为1;h为滤波器单位脉冲响应向量;t为[0:n-1]’,是h对应的时间向量。当函数输出默认时,绘制滤波器脉冲响应图;当n为默认时,函数自动选择n值。
(5)绘制离散系统的零极点图
函数:zplane功能:已知数字滤波器零极点z和p,或系数b和a,绘出系统的零极点图
调用格式:zplane(z,p);zplane(b,a);说明:z和p为零极点向量(为复数),b和a为滤波器的系数(为实数)。函数在Z平面绘出零点和极点。极点用x表示,零点用O表示。
滤波使用函数
(1)数字滤波
函数:filter功能:实现数字滤波器对数据的滤波
调用格式:y=filter(b,a,x);[y,zf]= filter(b,a,x,zi)说明:b和a为数字滤波器系数;x为滤波器的离散输人;y为滤波器的输出;y为与x具有相同大小的向量;zi和zf表示x在分段滤波时的初始状态和最终状态(在3.7.8小节中有更详细的介绍)。
(2)零相位数字滤波
函数:filtfilt;功能:实现数字滤波器对数据进行零相位滤波
调用格式:Y=filtfilt(b,a,x);说明:b和a为数字滤波器系数;x为滤波器的离散输人;y为滤波器的输出
数字陷波器和全相位数字滤波器的相关函数
(1)数字陷波器设计
函数:iirnotch功能:数字陷波器设计
调用格式:[b,a=iirnotch(Wo,Bw)说明:Wo是陷波器的中心频率,Bw是陷波器的带宽,参数b和a是数字滤波器的系数
(2)全相位数字滤波器
函数:iirgrpdelay;功能:设计一个在已知频率范围内接近规定群延迟的全相位数字滤波器
调用格式:[b,a] = iirgrpdelay(N,F,Edges , Gd)说明:N是全通滤波器的阶数(N必须是偶数):F和Gd是指在某一频率区间内群延迟的指标,F和Gd都是矢量,Gd是通带的边缘值;b和a是数字滤波器的系数。
7.数字系统留数的计算函数
函数:residuez;功能:计算数字系统的留数和极点或计算数字系统系数。
调用格式:[r,p,k]=residuez(B,A);[B,A]=residuesz(r,p,k);说明:调用格式为[r,p,k]=residuez(B,A)它和模拟系统的residue函数类似。这里是当已知数字系统的极点向量p,留数向量r和剩余多项式k时,计算出数字系统的系数B和A。
三、案例分析
1.设计一个椭圆型滤波器,通带截止频率为3400Hz,阻带截止频率为3700Hz,通带波纹为0.8dB,阻带衰减为50dB。先设计成模拟滤波器,再经双线性Z变换转换成数字滤波器。用该滤波器处理San2.wav文件,把语音中的啸叫声过滤掉。程序清单如下:
clear all; close all; clc;
[y,Fs]=audioread('San2.wav'); % 读入数据
fc=3400;fb=3700; % 设置通带和阻带截止频率
Rp=3;Rs=60; % 设置通带波纹和阻带衰减
wp=2*pi*fc/Fs;ws=2*pi*fb/Fs; % 计算归一化频率
Ts=1/Fs;
Wp=2/Ts*tan(wp/2.);Ws=2/Ts*tan(ws/2.);% 模拟频率进行预畸
[M,Wn]=ellipord(Wp,Ws,Rp,Rs,'s'); % 得到模拟滤波器原型阶数和带宽
[bs,as]=ellip(M,Rp,Rs,Wn,'s'); % 得到模拟滤波器系数
[b,a]=bilinear(bs,as,Fs); % 双线性Z变换得数字滤波器系数
x=filter(b,a,y); % 对数据进行滤波
Y=fft(y); % 求输入和输出信号的谱图
X=fft(x);
N=length(x);
n2=1:N/2;
freq=(n2-1)*Fs/N;
% 作图
[H,ff]=freqz(b,a,1000,Fs); % 观察数字滤波器的响应曲线
plot(ff,20*log10(abs(H)),'k');
xlabel('频率/Hz'); ylabel('幅值/dB')
title('椭圆滤波器幅频响应曲线');
ylim([-80 10]); grid;
set(gcf,'color','w');
figure
subplot 211; plot(freq,abs(Y(n2)),'k'); % 输入信号谱图
xlabel('频率/Hz'); title('输入信号谱图')
subplot 212; plot(freq,abs(X(n2)),'k'); % 输出信号谱图
xlabel('频率/Hz'); title('输出信号谱图')
set(gcf,'color','w');
2.若低通滤波器的截止频率或带通滤波器的中心频率相对采样频率低(fc/fs),或者滤波器的过渡带很窄,则在设计滤波器后往往可发现滤波器的阶数很大,或是滤波器的响应曲线很差。在这种情况,通过降低信号的采样频率,再以降低的采样频率设计数字滤波器,会取得较好的效果。
例子：信号从bzsdata.mat文件读入,它由数个正弦信号与随机噪声叠加而成,采样频率为250Hz。信号中有一个5Hz的正弦信号,要求设计一个巴特沃斯滤波器,通带频率为 1.5~10Hz,阻带频率为1Hz和12Hz,而Ap=3,As=15。程序清单如下:
clear all; clc; close all;
load bzsdata.mat % 读入数据
N=length(bzs); % 原始数据长
t=(0:N-1)/Fs; % 设置时间
x=resample(bzs,1,5); % 降采样
N1=length(x); % 降采样后的长度
fs=Fs/5; % 降采样后的采样频率
fs2=fs/2; % 降采样后采样频率的一半
t1=(0:N1-1)/fs; % 降采样后的时间刻度
fp1=[1.5 10]; % 通带频率
fs1=[1 12]; % 阻带频率
wp1=fp1/fs2; % 归一化通带频率
ws1=fs1/fs2; % 归一化阻带频率
Ap=3; As=15; % 通带波纹和阻带衰减
[n,Wn]=buttord(wp1,ws1,Ap,As); % 求滤波器原型阶数和带宽
[bn1,an1]=butter(n,Wn); % 求数字滤波器系数
[H,f]=freqz(bn1,an1,1000,fs); % 求数字滤波器幅频曲线
y1=filter(bn1,an1,x); % 对降采样后的数据进行滤波
y=resample(y1,5,1); % 对滤波器输出恢复原采样频率
% 作图
figure(1)
subplot 311; plot(t,bzs,'k');
xlabel('时间/秒'); title('原始数据波形')
subplot 312; plot(t1,x,'k');
xlabel('时间/秒'); title('降采样后数据波形')
subplot 313; plot(t,y,'k')
xlabel('时间/秒'); title('滤波后数据波形')
set(gcf,'color','w');
figure(2)
plot(f,abs(H),'k');
grid; axis([0 25 0 1.1]);
xlabel('频率/Hz'); ylabel('幅值')
title('巴特沃斯滤波器的幅值响应')
set(gcf,'color','w');
这一个问题实际上并不是太难解决的,但:①采样频率不是用Hz,而是用1/min,这与我们平时用的单位稍有不同,会不太习惯;②在要求数据滤波时,只给出一个通带的要求,没有给出其他指标,如何来设计滤波器呢?
在工程应用中,一些具体的测量其采样率往往不是以Hz(1/s)为单位,而是以1/min、1/h、1/d(天)或一些空间长度1/cm、1/m,或速度、加速度等为单位的。但在处理中,不论采样率是什么样的单位,都可以当作单位为Hz进行处理。
同样是在工程上往往说不出滤波器的具体指标,只给出一个通带的要求。我们在进行滤波器设计时,可以先假定某一种滤波器,将该滤波器的阶数设为4阶,然后进行滤波器的设计并进行数据滤波处理。观察处理后是否满足工程上的要求,如果不满足,则可以增减滤波器阶数或改变滤波器类型等来满足工程上的要求。
在本例中,设采样率是1个样点/min,则信号周期为10~20h对应的频率为1/1200~1/600(1/min),这样fs/f0之比大约在1000的量级,显然太大了。所以按例3-7-4所述把信号降采样,降到1个样点/h(即要降采样率60倍),选用切比雪夫Ⅱ型滤波器。
滤波器的通带频率在采样频率为1/h的条件下,fp=[1/201/10]=[0.050.1],我们把阻带频率设为fs=[0.0250.15],通带波纹设为Ap=1,阻带衰减为As=50。有了这些条件就能设计滤波器了。程序清单如下:
clear all; clc; close all;
load jandatas.mat % 导入数据
N1=length(z); % 原始数据长
Fsm=1; % 原始采样频率1分钟1样点
y=resample(z,Fsm,60); % 降采样60倍
N=length(y); % 降采样后的长度
hour=0:N-1;
Fsh=1; % 降采样后采样频率1小时1样点
fp=[0.05 0.1]; % 通带频率
fs=[0.025 0.15]; % 阻带频率
Ap=1; As=50; % 通带波纹和阻带衰减
Wp=fp*2/Fsh; Ws=fs*2/Fsh; % 归一化通带和阻带频率
[M,Wn]=cheb2ord(Wp,Ws,Ap,As); % 求滤波器原型阶数和带宽
[bn,an]=cheby2(M,As,Wn); % 求数字滤波器系数
[H,f]=freqz(bn,an,1000,1); % 求数字滤波器幅频曲线
x=filter(bn,an,y); % 对降采样后的数据进行滤波
xx=resample(x,60,1); % 对滤波器输出恢复原采样频率
xx=xx(1:N1); % 求取与输入序列相同长度和单位
% 作图
figure(1)
plot(f,20*log10(abs(H)),'k');
axis([0 0.2 -70 10]); grid;
title('椭圆型滤波器幅频响应曲线');
xlabel('时间/小时'); ylabel('幅值/dB');
set(gcf,'color','w');
figure(2)
subplot 211; plot(hour,y,'k');
xlim([0 max(hour)]);
title('降采样后的数据'); xlabel('时间/小时');
subplot 212; plot(minute/10000,xx,'k');
title('滤波后周期为10-20小时的数据');xlabel('时间/万分钟');
set(gcf,'color','w');
3.如何使用数字陷波器滤除工频信号
在实际测量时经常会受到工频信号(交流50Hz)的干扰,有时干扰还很大,有用信号完全被淹没了。可以应用数字陷波器来消除工频信号的干扰。有一个心电信号,但测量时由于设备的问题受到工频信号的干扰,并完全被干扰所淹没请设计一个数字陷波器来恢复心电信号。心电信号的数据在noisyecg.mat中,设计数字陷波器滤除干扰噪声。程序清单如下:
clear all; clc; close all;
load('noisyecg.mat'); % 读入信号数据和采样频率
x=noisyecg; % 信号为x
N=length(x); % 信号长N
t =(0:N-1)./fs; % 时间刻度
fs2=fs/2; % 设置奈奎斯特频率
W0=50/fs2; % 陷波器中心频率
BW=0.1; % 陷波器带宽
[b,a]=iirnotch(W0,BW); % 设计IIR数字陷波器
[H,w]=freqz(b,a); % 求出滤波器的频域响应
y=filter(b,a,x); % 对信号滤波
% 作图
plot(w*fs/2/pi,20*log10(abs(H)),'k');
xlabel('频率/Hs'); ylabel('幅值/dB');
title('陷波器幅值响应图')
axis([0 125 -45 5]); grid;
set(gca,'XTickMode','manual','XTick',[0,40,50,60,100])
set(gca,'YTickMode','manual','YTick',[-40,-30,-20,-10,0])
set(gcf,'color','w');
figure(2)
subplot 211; plot(t,x,'k','linewidth',0.5);
xlabel('时间/s'); ylabel('幅值');
title('带噪心电波形图')
subplot 212; plot(t,y,'k');
xlabel('时间/s'); ylabel('幅值');
title('消噪后心电波形图')
set(gcf,'color','w');
从图中可以看到,在滤波前,带噪信号中的噪声完全淹没了心电信号,除了几个峰值外什么也看不出来;但在滤波后恢复了心电信号,但在初始部分还是有瞬态现象存在。有时,在工频干扰中不限于只有基频50Hz,往往还带有它的谐波。所以在处理之前先做谱分析,观察除基波外还有几个谐波。除基波的陷波器外,还可设计谐波的陷波器,把它们级联在一起,以滤除基波和谐波。有时处理的不是50Hz工频信号,而是在信号中混有其他的正弦信号。这时先要做谱分析,甚至对干扰的正弦信号要用校正法(将在第6章中介绍)求出该信号的频率,才能进一步进行陷波处理。
获取代码请关注MATLAB科研小白的个人公众号（即文章下方二维码），并回复数字滤波器的设计本公众号致力于解决找代码难，写代码怵。各位有什么急需的代码，欢迎后台留言~不定时更新科研技巧类推文，可以一起探讨科研，写作，文献，代码等诸多学术问题，我们一起进步。