脉冲响应不变法将模拟滤波器转换成为数字滤波器的套路

已知一个模拟滤波器系统，要将此模拟滤波器系统转换成为数字滤波器，方法会有很多种，在Matlab中也有专门的函数（impinvar 脉冲响应不变）来实现此功能，如果我们要自己编写算法来做这个过程的话可以考虑使用脉冲响应不变法来实现跟impinvar函数相同的效果。

脉冲响应不变的原理

1】对模拟滤波器进行拉普拉斯逆变换

理论依据是：对传输系统的单位冲击响应求拉普拉斯变换得到的结果就是传输系统的传输函数，因此对传输函数进行拉普拉斯逆变换就是传输系统的单位冲击响应。

2】对单位冲击响应进行采样*时间间隔T，得到等价的脉冲响应序列

理论依据&目的：使用 数字滤波器的响应序列 模仿 模拟滤波器的脉冲响应

3】对得到的数字滤波器脉冲响应序列进行Z变换，得到数字滤波器

理论依据：Z变换就是离散时间信号的拉普拉斯变换，相应目的见第一条

例子

假如我想使用巴特沃兹方法设计一个模拟滤波器，然后使用脉冲响应不变方法将这个模拟滤波器转成数字滤波器，已知通带导带、以及纹波系数、衰减增益等参数，那么按照以下步骤可以完成数字滤波器的设计。

1】先设计满足要求的模拟滤波器

1.1 用butterword指令先计算出滤波器的阶数以及截至频率

去掉下面指令中的‘s’选项，计算出来的结果直接就是数字滤波器的相关结果，可以省去后面模拟滤波器转数字滤波器的过程。

[N,Wn]=buttord(wp,ws,Rp,Rs,'s'); % 求巴特沃斯滤波器阶数

1.2 将计算出来的N\Wn给butter指令计算出s域的多项式传输方程的零极点

[z,p,k]=butter(N,Wn,'s');        % 求巴特沃斯滤波器的零极点以及增益

1.3 用zp2tf计算s域传输方程的多项式系数

[bb,ab] = zp2tf(z,p,k);

或者将1.2跟1.3合并，直接计算s域传输方程的系数

[bb,ab]=butter(N,Wn,'s');        % 求巴特沃斯滤波器系数

2】将模拟滤波器转换成数字滤波器（脉冲响应不变法 或 留数法）

[bd,ad] = impinvar(bb,ab,Fs); % 脉冲响应不变法

或者

[bd,ad]=residuez(T*Ra, pd, ha);	% 用原留数Ra和数字极点pd求得数字滤波器系数

3】 freqs 和 freqz 指令分别查看 s 域 z 域传输函数的频响曲线

freqs(bb,ab) % 接受s域参数，指数横轴坐标单位rad/s，与Hz相差一个2*pi系数

freqz（bd，ad，n，fs） % 接受z域参数，绝对数横轴坐标 hz，横轴刻度（fs/（2*n）:fs/（2*n）:fs/2）
set(gca,'xaxis','log') % 将横轴坐标设置为指数坐标

实现

综上，脉冲时不变响应法设计IIR滤波器的matlab程序见附件

附录

附使用matlab设计滤波器的其他一些demo实例做参考

低通滤波器（可以通过 doc lowpass 在matlab help文档中获取）

clear all;close all;clc;

d = fdesign.lowpass('Fp,Fst,Ap,Ast',0.15,0.25,1,60);
Hd = design(d,'equiripple');
designmethods(d)        % 查看可用的设计方法
fvtool(Hd)

n = 0:159;
x = (0.25*cos((pi/8)*n)+sin((pi/4)*n));
y = filter(Hd,x);
freq = 0:(2*pi)/160:pi;
xdft = fft(x);
ydft = fft(y);
figureplot(freq/pi,abs(xdft(1:length(x)/2+1)))
hold on
plot(freq/pi,abs(ydft(1:length(y)/2+1)))
hold off
legend('Original Signal','Filtered Signal')
ylabel('Magnitude')
xlabel('Normalized Frequency (\times\pi rad/sample)')

其中

Ap — amount of ripple allowed in the pass band in decibels (the default units). Also called Apass.
Ast — attenuation in the stop band in decibels (the default units). Also called Astop.
F3db — cutoff frequency for the point 3 dB point below the passband value. Specified in normalized frequency units.
Fc — cutoff frequency for the point 6 dB point below the passband value. Specified in normalized frequency units.
Fp — frequency at the start of the pass band. Specified in normalized frequency units. Also called Fpass.
Fst — frequency at the end of the stop band. Specified in normalized frequency units. Also called Fstop.
N — filter order.
Na and Nb are the order of the denominator and numerator.

 

高通滤波器、带通滤波器都可以通过同样的方式从matlab help中获取demo代码。

 

clear all;close all;clc;

d = fdesign.highpass('Fst,Fp,Ast,Ap',20,30,60,1,16000);
Hd = design(d,'ifir');
designmethods(d)        % 查看可用的设计方法
fvtool(Hd)

% 生成待滤波信号
n = 0:0.00001:1;
x = (3*cos(2000*(2*pi)*n)+sin(10*(2*pi)*n));
y = filter(Hd,x);

figure
plot(x)
hold on
plot(y)

根据模拟滤波器的特性参数，设计数字滤波器（巴特沃斯滤波器 切比雪夫I型滤波器 切比雪夫II型滤波器 椭圆型滤波器 ）

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) % 显示滤波器阶数
[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')
axis([0 max(we/pi) -30 2]); %grid;
legend('巴特沃斯滤波器','切比雪夫I型滤波器','切比雪夫II型滤波器','椭圆型滤波器')
xlabel('角频率{\omega}/{\pi}'); ylabel('幅值/dB')
set(gcf,'color','w'); 

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','--');

impinvar函数的使用方法