【数字信号处理】频率响应进阶——幅值、相位、群延迟一站式计算【含matlab代码】-拓冰建站

第四篇：频率响应进阶——幅值、相位、群延迟一站式计算

在前三篇中，我们深入剖析了线性相位 FIR 的幅度响应 (H_r(\omega))，但它只是频率响应的“骨架”。在实际工程中，我们更常需要观察幅频响应（dB 单位）、相频响应以及群延迟，以便全面评估滤波器的选频特性、相位失真和实时性。MATLAB 自带的freqz功能强大，但输出格式不够工程化。为此，我们的freqz_m.m对freqz进行了封装，将频率限制在 ([0,\pi])，归一化幅值，同时一次性输出四种常用曲线。本篇将逐行解析该函数，并带你掌握群延迟的物理意义和计算方法。

1. 为什么要自定义`freqz_m`？

标准freqz(b,a,N,'whole')会返回 (N) 个频率点（覆盖 (0\sim 2\pi)），且幅值未归一化，相位为弧度，群延迟需要单独调用grpdelay。在工程绘图时，我们通常只需要 ([0,\pi]) 区间，幅值用 dB 表示并归一化到 0 dB 峰值，相位展开更方便观察。freqz_m统一了这些需求，一次调用即可获得：

db：归一化幅度（dB）
mag：线性幅度（绝对值）
pha：相位（弧度，未展开）
grd：群延迟（样本数）
w：频率向量（(0\sim\pi)）

此外，它支持任意 IIR 或 FIR 滤波器（系数b,a），通用性极强。

2. 代码逐行解读

function[db,mag,pha,grd,w]=freqz_m(b,a);%Modified version of freqz subroutine

输入：分子系数b，分母系数a（对于 FIR，a=1）。
输出：五个变量。

第一行：调用标准freqz，使用 1000 个点覆盖整个单位圆（‘whole’），得到复数频率响应H和频率向量w。

[H,w]=freqz(b,a,1000,'whole');

第二行：只取前 501 个点（对应 (0\sim\pi)），并转置为行向量（便于后续矩阵运算）。这里固定 1000 点，等效于 501 个频率点，分辨率约为 (\pi/500)。

H=(H(1:501))'; w = (w(1:501))';

第三行：计算线性幅度（绝对值）。

mag=abs(H);

第四行：计算归一化 dB 幅度。注意，这里使用了(mag+eps)/max(mag)的技巧——加eps防止mag=0时出现log10(0)负无穷，同时将最大值归一化为 1（即 0 dB）。所有幅度都相对于峰值。

db=20*log10((mag+eps)/max(mag));

第五行：直接用angle计算相位（主值区间 (-\pi\sim\pi)）。

pha=angle(H);

第六行：调用 MATLAB 的grpdelay计算群延迟，传入b,a和频率向量w（注意grpdelay接受的频率单位是弧度/样本，w正是如此）。返回的grd就是对应每个频率点的群延迟（单位：样本）。

grd=grpdelay(b,a,w);

3. 群延迟的本质与工程意义

群延迟定义为相位对频率的负导数：

[
\tau_g(\omega) = -\frac{d\phi(\omega)}{d\omega}
]

它表示某一频率分量通过滤波器时，其包络发生的时延。对于线性相位滤波器，群延迟为常数（等于(M-1)/2），意味着所有频率分量的包络延迟相同，信号波形不会发生相位畸变。而非线性相位滤波器（如普通 IIR）的群延迟随频率变化，会导致输出波形失真（如脉冲展宽）。

grpdelay函数采用数值差分法估计导数，对于 FIR 滤波器，其理论值应为(M-1)/2，但实际计算会因频率离散而略有波动。我们可通过freqz_m绘图验证。

4. 实战演练：对比 FIR 与 IIR 的频率响应

4.1 设计一个低通 FIR（窗函数法）

% 设计一个 31 阶（M=31）低通 FIR，截止频率 0.4*piwc=0.4*pi;M=31;h=fir1(M-1,wc/pi,'low',hamming(M));% fir1 自动生成偶对称 Type-1[db,mag,pha,grd,w]=freqz_m(h,1);% 绘制四合一图figure;subplot(2,2,1);plot(w/pi,db);grid;xlabel('\omega/\pi');ylabel('dB');title('归一化幅度 (dB)');subplot(2,2,2);plot(w/pi,mag);grid;xlabel('\omega/\pi');ylabel('|H|');title('线性幅度');subplot(2,2,3);plot(w/pi,pha);grid;xlabel('\omega/\pi');ylabel('Radians');title('相位');subplot(2,2,4);plot(w/pi,grd);grid;xlabel('\omega/\pi');ylabel('Samples');title('群延迟');

观察群延迟图，在通带内grd稳定在(M-1)/2 = 15样本，验证了线性相位特性。

4.2 对比一个 IIR 椭圆滤波器

[b,a]=ellip(6,1,40,0.4);% 6阶椭圆低通[db_iir,mag_iir,pha_iir,grd_iir,w]=freqz_m(b,a);figure;plot(w/pi,grd_iir);grid;xlabel('\omega/\pi');ylabel('群延迟 (样本)');title('IIR 椭圆滤波器群延迟 —— 明显波动');