MATLAB 提供了函數filter，可以實現差分方程的遞規求解。
設差分方程的形式為\(a_0y(n) + a_1y(n-1) + \cdots + a_my(n-m)=b_0x(n)+b_1x(n-1)+\cdots+b_sx(n-s)\)
基本的調用方法：
yn = filter(B, A, xn)
- \(B = [b_0, \cdots, b_s], A = [a_0, \cdots, a_m]\)；
xn是輸入信號；yn是輸入信號通過系統的零狀態響應。
如果輸入是單位脈沖函數，則輸出就是系統單位脈沖響應。
yn = filter(B, A, xn, xi)
- 前3個參數的用法同上。
xi是等效初始條件序列，一般是通過調用函數filtic得到的。yi是由系統初始狀態引起的零輸入相應和輸入信號引起的零狀態相應之和，即全響應。
xi = filtic(B, A, ys, xs)
- 參數
B,A的用法如上 ys和xs表示的輸入和輸出在\(n_0\)時刻之前的初始條件：\(ys =[y(-1),y(-2),\cdots, y(-N)]\)，\(xs=[x(-1),x(-2),\cdots, x(-N)]\)。如果系統是因果系統，顯然輸入為0，即xs = 0，此時可以省略xs。
例子
- 求系統的單位脈沖響應和零輸入相應
% 滑動平均濾波器的差分方程如下
% y(n) = 1/5 * (x(n) + x(n-1) + x(n-2) + x(n-3) + x(n-4))
windowSize = 5;
B = ones(1, windowSize) / windowSize;
A = 1;
% 求該系統的單位脈沖響應
subplot(411)
delta = [1, zeros(1,33)]; % 單位脈沖信號
y0 = filter(B, A, delta);
stem(0:length(y0)-1, y0)
title('單位脈沖響應')
% 求系統在某個輸入信號下的零狀態響應
% 輸入信號
subplot(412)
xn = [ones(1, 32), zeros(1,4)];
xn(13)=2; xn(16)=0.5;xn(20)=1.5;
stem(0:length(xn)-1, xn);
title('輸入信號')
% 零狀態相應
subplot(413)
yn = filter(B, A, xn, xi);
stem(0:length(yn)-1, yn);
title('零狀態相應輸出')
% 卷積驗證
subplot(414)
yc = conv(xn, y0)
stem(0:length(yc)-1, yc)
title('卷積驗證')
xlim([0, 35])
運行結果：
- 求系統的全響應
% 設因果系統 y(n) = 0.8*y(n-1) + x(n)，輸入為x(n)=\delta(n)
% 在不同初始條件下的輸出：（1）y(-1)=0; (2) y(-1)=1
A = [1, -0.8];
B = 1;
xn = [1, zeros(1, 30)];
% y(-1)=0, 輸出即是零狀態輸入
subplot(211)
y0 = filter(B, A, xn);
stem(0:length(y0)-1, y0)
title('初始狀態為y(-1)=0')
ylim([0, 2])
% y(-1)=1
subplot(212)
ys = 1; % 初始狀態
xi = filtic(B, A, ys);
yn = filter(B, A, xn, xi);
stem(0:length(yn)-1, yn);
title('初始狀態為y(-1)=1')
ylim([0, 2])
運行結果
