第四篇频率响应进阶——幅值、相位、群延迟一站式计算在前三篇中我们深入剖析了线性相位 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归一化幅度dBmag线性幅度绝对值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对于 FIRa1。输出五个变量。第一行调用标准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));第三行计算线性幅度绝对值。magabs(H);第四行计算归一化 dB 幅度。注意这里使用了(mageps)/max(mag)的技巧——加eps防止mag0时出现log10(0)负无穷同时将最大值归一化为 1即 0 dB。所有幅度都相对于峰值。db20*log10((mageps)/max(mag));第五行直接用angle计算相位主值区间 (-\pi\sim\pi)。phaangle(H);第六行调用 MATLAB 的grpdelay计算群延迟传入b,a和频率向量w注意grpdelay接受的频率单位是弧度/样本w正是如此。返回的grd就是对应每个频率点的群延迟单位样本。grdgrpdelay(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 阶M31低通 FIR截止频率 0.4*piwc0.4*pi;M31;hfir1(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 椭圆滤波器群延迟 —— 明显波动);你会看到 IIR 的群延迟在通带内剧烈变化这就是相位失真的来源。5. 与前几篇的衔接我们已经有了ampl_ress计算幅度响应 (H_r)有符号而freqz_m给出的是线性幅度mag |H(e^{jω})|。两者关系为[\text{mag} |H_r(\omega)|]当 (H_r) 为负值时相位会出现 (-\pi) 的跳变这称为“相位卷绕”。而freqz_m的相位输出是angle(H)会自动处理这些跳变。我们可以在同一幅图上比较Hr和mag加深对符号的理解。6. 工程技巧与注意事项归一化基准db输出中最大值始终为 0 dB便于观察阻带衰减如 -60 dB 表示衰减 1000 倍。防除零eps的加入非常关键否则在阻带深度零点处log10(0)会产生-Inf。频率点数固定 1000 点实际 501 点已经足够光滑若需更高分辨率可修改内部数字但建议保持默认。群延迟计算grpdelay对 FIR 计算准确但对 IIR 可能因相位非线性导致数值波动较大但结果可信。7. 总结与下篇预告本篇我们学会了使用freqz_m一站式获取五种频率响应指标理解群延迟的定义及其对保真度的影响通过实例对比了线性相位 FIR 与非线性相位 IIR 的差异。至此我们完成了滤波器分析工具的全部讲解ampl_ressfreqz_m。接下来我们将进入滤波器结构转换的世界先从直接型→级联型开始探索如何将高阶滤波器拆解为稳定的二阶节。所有代码均已打包点击下方链接免费获取下载链接下篇预告滤波器结构转换一——直接型与级联型的互转。我们将深入dir2cas.m剖析零极点配对算法并揭示为什么级联型在高阶滤波器中更受青睐。敬请期待