$方法一：利用\text{freqz()}$

[h, w] = freqz(b, a);

%%  利用freqz()求解系统的频率响应
b = [0.1321 0 -0.3963 0 0.3963 0 -0.1321];
a = [1 0 0.34319 0 0.60439 0 0.20407];
figure('color','w');

subplot(2,1,1)
[h, w] = freqz(b, a);
plot(w/pi, abs(h), w/pi, angle(h));
title('系统幅度和角度响应');
xlabel('归一化频率\omega /\pi');
ylabel('幅度/角度');
legend('幅度', '角度'); 

N = 13;n = 0:N-1;
h2 = impz(b, a, n);
subplot(2,1,2)
stem(n, h2,'filled');
title('系统冲激响应');
xlabel('n');
ylabel('幅度');

$方法二：利用自定义\text{freqz\_m()}$

有自定义函数：

function [db,mag,pha,grd,w] = freqz_m(b,a)
%% db  记录一组对应[0,pi]频率区域的相对幅值响应值
%% mag 记录一组对应[0,pi]频率区域的绝对幅值响应值
%% pha 记录一组对应[0,pi]频率区域的相位响应值
%% grd 记录一组对应[0,pi]频率区域的群迟延响应值
    [H, w] = freqz(b,a,1000,"whole");
    H = (H(1:501))';
    w = (w(1:501))';
    mag = abs(H);
    db = 20*log10((mag+eps)/max(mag));
    pha = angle(H);
    grd = grpdelay(b,a,w);
end

[db,mag,pha,grd,w] = freqz_m(b,a)

%% 利用freqz_m求解系统相对幅频响应，绝对幅频响应，相位频率响应和群延迟

b = [0.1321 0 0.3963 0 0.3963 0 0.1321];
a = [1 0 -0.34319 0 0.60439 0 -0.20407];

[db,mag,pha,grd,w] = freqz_m(b,a);

figure('color','w');
subplot(221)
plot(w/pi,mag,'color','#FFBB00','LineWidth',1);grid on;grid minor;

subplot(222)
plot(w/pi,pha,'color','#C9007A','LineWidth',1);grid on;grid minor;

subplot(223)
plot(w/pi,db,'color','#A5EF00','LineWidth',1);grid on;grid minor;

subplot(224)
plot(w/pi,grd,'color','#1924B1','LineWidth',1);grid on;grid minor;