如何使用 Matlab 的滤波器命令绘制系统的脉冲响应?

How can I use Matlab's filter command to plot impulse response of system?

我试图通过使用 Matlab 的滤波器命令而不是其他预先存在的 Matlab 函数来绘制 IIR 和 FIR 系统的脉冲响应图。过滤器必须能够处理样本输入,例如 plot([1 2 2], [0 1 .8])。问题看起来像这样:

function iplot(b, a)

% IPLOT Impulse response of system.
% IPLOT(B,A,N) is the N-point impulse response
% of the filter B/A:

%                          -1                -nb 
%     
%           B(z) b(1) + b(2)z + .... + b(nb+1)z
%
%    H(z) = ---- = --------------------------------- 
% 
%                          -1                -n         
%           A(z) a(1) + a(2)z + .... + a(na+1)z


% If h[n] FIR, then N = length(h);

% given numerator and denominator coefficients in vectors B and A.

% N is specified according to the following rule:

% If If h[n] is IIR and increasing (i.e. |h[n]|-->inf as n-->inf),

% then N=20;

% h[n] is IIR and decreasing (i.e. |h[n]|-->0 as n-->inf), 

% then the maximum N is determined such that

% rms value of h(1:N) = 0.999 * rms value of h(1:1000).

% However, in this case, N must also be chosen such that 10 <=

% N <= 100

这就是我目前所拥有的。我知道它不正确,现在它没有绘制任何东西,但据我所知。任何帮助深表感谢!

h = filter(a,b,[1,zeros(1,999)]); 
N = length(h); 

plot(h, zeros(length(h),999), 'b'); 

poles = roots(a); 

if poles <= 1

    N = 20; 
    plot(h, N, 'b'); 

end


end

第一行代码正确。但是调用filter时需要反转ba的顺序。 b代表分子系数,a代表分母系数。这也在你的文档字符串中。

您正确地将输入指定为单个输入。但是,您的绘图代码不正确。 plot 的工作原理是它接受一组 x 值和一组 y 值,我们将这些点绘制成对。因为您正在绘制离散输出信号,所以我建议您使用 stem 而不是 plot。这样,每个点都用一条从水平轴到感兴趣点的线绘制。

无论如何,您需要提供一组 x 个值。具体来说,您需要一个从 0 到 999 的向量。y 值由 filter 引起的输出脉冲响应组成。目前,您的 x 值包含输出脉冲响应,而您的 y 值是一堆零。这将被可视化为聚集在 y = 0 处的一大堆点,这可能不是您想要的。

因此您只需要这样做:

h = filter(b,a,[1,zeros(1,999)]); 
plot(0:999, h, 'b');  %// Change

关于我的建议,我会使用 stem 代替:

stem(0:999, h, 'b');

代码的最后一部分是绘制滤波器的极点。您还应该制作一个 separate 图形并绘制这些结果。你没有这样做,所以当你尝试绘制极点时,包含你的脉冲响应的图形 将被覆盖 。在当前 window 打开的情况下多次调用 plot 将覆盖最后一次调用 plot 的情节。因此,在继续之前生成一个新的 figure window:

figure; %// New window
h = filter(b,a,[1,zeros(1,999)]); 
stem(0:999, h, 'b');  %// Change

figure; %// New window
poles = roots(a); %// Find roots
%// Plot the roots as single dots
plot(real(poles), imag(poles), 'b.');

但是,这并没有回答您的问题陈述。问题陈述说:

If h[n] is FIR, then N = length(h); given numerator and denominator coefficients in vectors B and A.

N is specified according to the following rule: If h[n] is IIR and increasing (i.e. |h[n]|-->inf as n-->inf), then N=20;

h[n] is IIR and decreasing (i.e. |h[n]|-->0 as n-->inf), then the maximum N is determined such that the rms value of h(1:N) = 0.999 * rms value of h(1:1000). However, in this case, N must also be chosen such that 10 <= N <= 100

我们知道 FIR 滤波器是指当您指定 A 为 1 或只是一个值时。当检测到A的长度为1时,则N = length(h);。但是,如果 a 而不是 只是一个值,我们必须检查两种情况。第一种情况是脉冲响应增加时。您可以通过使用 diff 并确保您的脉冲响应的所有相邻差异都是正数来检查这一点。下一种情况是 h 正在减少。你可以通过确保你的脉冲响应的所有相邻差异都是负数来检查这一点。

如果是递增的大小写,则选择N = 20。如果是递减的情况,则选择 N 值为您在上面看到的值,您可以使用 MATLAB 中的 rms 函数来帮助您做到这一点。我们还需要确保这个值在 10 到 100 之间,并且我们还希望 舍入 这个值,因为我们试图选择一个整数点和 rms 可能会给你一个浮点数。

因此,使用默认的 1000 点执行 filter,然后您需要 截断 响应,具体取决于函数的输入。首先我们需要检查 FIR,然后是 IIR 情况来决定显示 N 的哪些值。

因此:

figure; %// New window
h = filter(b,a,[1,zeros(1,999)]); 

%// Decide on the value of N
if length(a) == 1
    N = length(h);
else
    d = diff(h);
    if all(d >= 0) %// Check for increasing
        N = 20;
    else %// Decreasing case
        %// Find RMS value
        N = round(0.999*rms(h));

        %// Ensure that 10 <= N <= 100
        if N < 10
            N = 10;
        end
        if N > 100
            N = 100;
        end
    end
end

%// Plot just up to the first N values
stem(0:N-1, h(1:N), 'b');  

figure; %// New window
poles = roots(a); %// Find roots
%// Plot the roots as single dots
plot(real(poles), imag(poles), 'b.');