如何在 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 之间的那些系数是 0，并且有 20 个这些系数总共为零，这就是为什么我们需要在 1 和 -(20/21) 项之间用零填充向量。

我们得到：

如果要绘制滤波器的极点和零点，请组合使用 tf and pzmap:

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

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

我们得到这个：

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

Answer 3

来自 Matlab 的 filter 文档：

filters the input data, x, using a rational transfer function defined by the numerator and denominator coefficients b and a, respectively.

您可能知道，过滤脉冲输入会给您脉冲响应。然后剩下的就是获得那些 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）。

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

How do I plot the frequency response of a digital filter in MATLAB?

matlab

signal-processing