已知一個模擬濾波器系統，要將此模擬濾波器系統轉換成為數字濾波器，方法會有很多種，在Matlab中也有專門的函式（impinvar）來實現此功能，如果我們要自己編寫演算法來做這個過程的話可以考慮使用脈衝響應不變法來實現跟impinvar函式相同的效果。

1】對模擬濾波器進行拉普拉斯逆變換

理論依據是：對傳輸系統的單位衝擊響應求拉普拉斯變換得到的結果就是傳輸系統的傳輸函式，因此對傳輸函式進行拉普拉斯逆變換就是傳輸系統的單位衝擊響應。

2】對單位衝擊響應進行取樣*時間間隔T，得到等價的脈衝響應序列

理論依據&目的：使用 數字濾波器的響應序列 模仿 模擬濾波器的脈衝響應

3】對得到的數字濾波器脈衝響應序列進行Z變換，得到數字濾波器

理論依據：Z變換就是離散時間訊號的拉普拉斯變換，相應目的見第一條

參考網路上脈衝時不變響應法設計IIR濾波器的matlab程式如下：

clear all; clc; close all;
 
bs=[1,1];as=[1,5,6];            % 系統分子分母系數向量
Fs=10; T=1/Fs;                  % 取樣頻率和取樣間隔
[Ra,pa,ha]=residue(bs, as);        % 將模擬濾波器係數向量變為模擬極點和留數
pd=exp(pa*T);                    % 將模擬極點變為數字（z平面）極點pd
[bd,ad]=residuez(T*Ra, pd, ha);    % 用原留數Ra和數字極點pd求得數字濾波器係數
t=0:0.1:3;                      % 時間序列
ha=impulse(bs,as,t);            % 計算模擬系統的脈衝響應
hd=impz(bd,ad,31);              % 數字系統的脈衝響應
% 呼叫impinvar函式計算數字濾波器係數
[Bd,Ad]=impinvar(bs,as,Fs);
fprintf('bd=%5.4f   %5.4f   ad=%5.4f   %5.4f   %5.4f\n\n',bd,ad);
fprintf('Bd=%5.4f   %5.4f   Ad=%5.4f   %5.4f   %5.4f\n',Bd,Ad);
% 作圖
plot(t,ha*T,'r','linewidth',3); hold on; grid on;
plot(t,hd,'k');
legend('模擬濾波器脈衝響應','數字濾波器脈衝響應');
xlabel('時間/s'); ylabel('幅值/dB');
title('原模擬濾波器的脈衝響應與數字濾波器的脈衝響應比較')
set(gcf,'color','w')   

附使用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函式的使用方法

clear all; clc; close all;
 
K=1.74802;                           % K(r)中的引數值
w1=2*pi*9.15494;
w2=2*pi*2.27979;
w3=2*pi*1.22535;
w4=2*pi*21.9;
lemda=2*pi*4.05981;
Fs=400;                              % 取樣頻率
% 把K(r)中的各引數值轉為2個子系統的係數
b(1)=K*w1; b(2)=0;                   % 第1子系統分子
a(1)=1/w2; a(2)=1;                   % 第2子系統分子
c(1)=1; c(2)=2*lemda; c(3)=w1*w1;    % 第1子系統分母
d(1)=1/w3/w4; d(2)=1/w3+1/w4; d(3)=1;% 第2子系統分母
 
B=conv(b,a);                         % 求出模擬系統的分子的係數
A=conv(c,d);                         % 求出模擬系統的分母的係數
[Hs,whs]=freqs(B,A);                 % 求模擬系統響應曲線
 
[num,den]=impinvar(B,A,Fs);          % 脈衝不變法求出數字系統的分子和分母系數
[Hz,wz]=freqz(num,den);              % 求數字系統響應曲線
% 作圖
line(whs/2/pi,abs(Hs),'color',[.6 .6 .6],'linewidth',3); 
axis([0 30  0 1.2]); box on; hold on
plot(wz/pi*Fs/2,abs(Hz),'k');
title('K(s)模擬響應曲線和數字響應曲線比較');
xlabel('頻率/Hz'); ylabel('幅值')
legend('模擬系統','數字系統');
set(gcf,'color','w');