1. 程式人生 > >【 MATLAB 】使用 MATLAB 求由差分方程表示的濾波器的響應的兩種方法

【 MATLAB 】使用 MATLAB 求由差分方程表示的濾波器的響應的兩種方法

例題:

一個3階低通濾波器由下面差分方程描述:

y(n) = 0.0181 x(n) + 0.0543 x(n-1) + 0.0543 x(n-2) + 0.0181 x(n-3) + 1.76 y(n-1) - 1.1829 y(n-2) + 0.2781 y(n-3)

畫出這個濾波器的幅度和相位響應,並驗證它是一個低通濾波器。

第一種方法是博文裡給出的:

第二個案例中,類比,如果我知道一個LTI系統的脈衝響應h(n),那麼我也能求出它的頻率響應:

k = [0:M];
n = [n1:n2];
X = x * (exp(-j * pi/M)).^(n'*k);

博文裡面由具體的推薦,上述程式中X就是頻率響應,也就是x的DTFT,如果x換成h,則X可以換成H。

第二種方法是,通過差分方程直接求出系統的頻率響應,求解的方法是通過向量化的方法。

設某一LTI系統的差分方程表示為:

y(n)+ \sum_{l=1}^{N}a_ly(n-l)= \sum_{m = 0}^{M}b_mx(n-m)

可以用一種簡單的矩陣向量乘法來完成。如果在[0,\pi]k=0,1,...,K個等分頻率上求H(e^{jw}),那麼

注意,上述的b以及a向量都是行向量。

先給出第一種方法的指令碼:

clc
clear
close all

b = [0.0181,0.0543,0.0543,0.0181];
a = [1.0000,-1.7600,1.1829,-0.2781];
[h,t]=impz(b,a);
k = [0:500];
w = (pi/500)*k;
t = t';
h = h';
H = h * ( exp(-j*pi/500) ).^(t'*k);

magH = abs(H);
angH = angle(H);

subplot(2,1,1);
plot(w/pi,magH);
title('Magnitude part');


subplot(2,1,2);
plot(w/pi,angH);
title('Angle part');

這種方法的思路是通過差分方程可以得到有理傳遞函式或者頻率響應的分子和分母系數,通過impz函式得到脈衝響應,之後由脈衝響應h(n)得到頻率響應。

clc
clear
close all

b = [0.0181,0.0543,0.0543,0.0181];
a = [1.0000,-1.7600,1.1829,-0.2781];
m = 0:length(b)-1;
l = 0:length(a)-1;
k = 0:500;
w = (pi/500)*k;
nume = b * exp(-j * m' * w);
den = a * exp(-j * l' * w);
H = nume ./ den;
magH = abs(H);
angH = angle(H);

subplot(2,1,1);
plot(w/pi,magH);
title('Magnitude Response');

subplot(2,1,2);
plot(w/pi,angH);
title('Phase Response');

從圖可以看出這確實是一個低通濾波器。