如何在 MATLAB 中绘制数字滤波器的频率响应？

Question

我正在尝试绘制频率响应图。我被要求在 MATLAB 中使用 filter，但由于我已经阅读了手册，我仍然不明白它是如何执行 Z 转换的。

我有下面写的数字滤波器的脉冲响应：

for i=1:22;
   y(i)= 0;
end
x(22) = 1;
for k=23:2523
    x(k) = 0;
end
for n = 22:2522;
    y(n) = ((1/21)*x(n))+((20/21)*y(n-21));
end
plot(y);

只是y[n] = 1/21*x[n] + 20/21*y[n-21]的反馈系统

以下是我计算上述系统的Z 变换的计算，它最终决定了脉冲响应：

Z(y) = Z((1/21)*x(n)+(20/21)*y(n-21))

Y(Z) = (1/21)X(Z)+(20/21)*Z.^-21Y(Z)

Z(Z)-(20/21)*Z.^-21Y(Z) = (1/21)X(Z)

Y(Z)(1-(20/21)*Z.^-21) = (1/21)X(Z)  // divide by X(Z)*(1-(20/21)*Z.^-21)

Y(Z)/X(Z) = (1/21)/(1-(20/21)*Z.^-21)

H(Z) = (1/21)/(1-(20/21)*Z.^-21) // B = 1/21, A = 20/21

H(Z) = (B*Z.^21)/(Z.^21-A)

如何绘制H(Z) 的频率响应？我应该使用filter吗？

Answer 1

如果您只需要绘制脉冲响应，这很容易。脉冲响应是数字滤波器对狄拉克脉冲的响应。你已经有了差分方程，所以你已经在'z'并且你不关心's'，你不必执行's'到'z'的转换（这是一个主题本身！）。

所以只需生成一个由零组成的信号 x(n)，除了第一个样本 x(1) 为 1 之外。将它通过滤波器（是的，你的差分方程）。你得到的 y(n) 就是你的脉冲响应 h(n)。这基本上就是你所做的。

当然，如果你对这个 h(n) 进行 FFT，你会得到相位和幅度响应。

Answer 2

使用信号处理工具箱中的freqz（希望你拥有它）。您首先需要找到 Z 变换，您已经在这里完成了：

Y(Z)/X(Z) = (1/21)/(1-(20/21)*Z.^-21)

freqz 接受与传递函数的分子和分母相对应的系数向量。它是这样称呼的：

freqz(b, a);

b 和 a 是传递函数的分子和分母系数。这将产生一个显示上述系统幅度和相位响应（因此是频率响应）的图形。

因此，您只需要这样做：

b = 1/21;
a = [1 zeros(1,20) -(20/21)];
freqz(b, a)

特别注意a 向量。它有 1，后跟 20 个零，然后是 -(20/21)。因为你有一个 -21 次方的系数，除了它之前的 1 之外别无其他，这意味着 -1 到 -20 之间的那些系数是零，并且这些系数中有 20 个总共为零，这就是为什么我们需要在 1 和 -(20/21) 术语之间用零填充向量。

我们得到：

---

如果您想绘制滤波器的极点和零点，请使用tf 和pzmap 的组合：

sys = tf(b, a, -1);
pzmap(sys);

tf 通过指定滤波器的分子和分母系数来创建传递函数，-1 暗示它是一个离散时间滤波器，但我们不知道采样时间是多少。 pzmap 绘制 z 域中的极点和零点，单位圆重叠。

我们得到这个：

这是有道理的，因为您的系统中没有零点和 21 个极点，因为您在检查示例中的离散时间序列时有 21 个延迟元素。

Answer 3

来自 Matlab 的 filter 文档：

使用分别由分子和分母系数 b 和 a 定义的 rational transfer function 过滤输入数据 x。

您可能知道，过滤脉冲输入会得到脉冲响应。然后仍然需要获得那些b 和a 系数。 您可以直接从差分方程中获得那些

y[n] = 1/21*x[n] + 20/21*y[n-21];

（如rational transfer function link above 中所示）或等效地从您导出的有理传递函数：

%H(Z) = (B*Z.^21)/(Z.^21-A)

在任何一种情况下，您都应该得到以下a 和b 系数：

a = [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -20/21];
% or equivalently: a=zeros(22,1); a(1)=1; a(22)=-20/21;
b = [1/21];

因此，

% setup the impulse input
x = zeros(2500,1);
x(1) = 1;

% compute the impulse response
y = filter(b, a, x);

这个脉冲响应可以像你所做的那样绘制出来：

plot(y);

就这与变换H(Z) 的关系而言，您获得的封闭形式表达式可以根据Laurent series 展开进行评估，然后具有时域y 的系数（脉冲响应) 系列。

但是，H(z) 是在复平面收敛的|z| > R（在您的情况下为R=power(20/21,1/21)）区域中定义的解析函数。更典型的是绘制频率响应，它对应于在单位圆上评估的H(z)（即，对于满足|z|=1 或等效的z = exp(j * theta) 的复数，θ 在 [0-2pi] 范围内）。计算单位圆上规则间隔点处H(z) 的值的一种有效方法是对脉冲响应进行 FFT：

FrequencyResponse = fft(y);
figure(1);
plot(abs(FrequencyResponse));
figure(2);
plot(phase(FrequencyResponse));

P.S.：如果您有信号处理工具箱，filter 和 fft 的计算可以在单个调用 freqz 中完成（尽管特别要求您使用 filter )。