【数字信号处理含matlab代码】第三篇：线性相位 FIR 滤波器（二）——Type-3 与 Type-4 及自动识别

第三篇线性相位 FIR 滤波器二——Type-3 与 Type-4 及自动识别在上一篇中我们领略了偶对称Type‑1/2的优美余弦展开。然而工程中还存在另一类重要的线性相位滤波器——奇对称Type‑3/4它们天生适用于希尔伯特变换和微分器。这些滤波器的幅度响应必须用正弦级数表示且在低频乃至高频存在固有零点。更妙的是我们的ampl_ress.m仅凭对称性检查就能自动判断当前h属于哪一类型免去人工分类的烦恼。本篇将为你彻底揭开这些秘密。下载链接1. 奇对称家族的数学面貌奇对称的条件为 ( h(n) -h(M-1-n) )。这意味着滤波器具有反对称的冲激响应中心点若存在必为零。其频率响应同样可分解为[H(e^{j\omega}) H_r(\omega) \cdot e^{-j\omega \alpha} \cdot e^{j\beta}]与偶对称相比奇对称的相位中多了一个固定的 ( \pm \pi/2 ) 的相移正是实现 90° 移相的关键。我们关注的核心仍是实幅度响应 ( H_r(\omega) )但它现在由正弦函数组合而成。1.1 Type-3( M ) 为奇数奇对称设 ( M 2L 1 )。由于反对称中心点 ( h(L) 0 )。将冲激响应按中心配对[H(e^{j\omega}) \sum_{n1}^{L} h(L-n) \left( e^{j\omega n} - e^{-j\omega n} \right) e^{-j\omega L}]利用欧拉公式 ( e^{j\omega n} - e^{-j\omega n} 2j\sin(\omega n) )并提取公因子 ( e^{-j\omega L} )得到[H(e^{j\omega}) e^{-j\omega L} \cdot j \sum_{n1}^{L} 2h(L-n) \sin(\omega n)]将复数幅度中的 ( j ) 吸收到相位里定义幅度响应为实函数[\boxed{H_r(\omega) \sum_{n1}^{L} c(n) \sin(\omega n)}]其中系数 ( c(n) 2h(L-n) )( n 1,2,\dots,L )。关键特征当 ( \omega 0 ) 或 ( \omega \pi ) 时所有 ( \sin(\omega n) 0 )因此Type-3 在 0 和 π 处幅度必为零。它仅适用于带通、希尔伯特变换或微分器且不能通过低通/高通。1.2 Type-4( M ) 为偶数奇对称设 ( M 2L )。反对称中心位于半整数点[H(e^{j\omega}) \sum_{n1}^{L} h(L-n) \left( e^{j\omega (n-0.5)} - e^{-j\omega (n-0.5)} \right) e^{-j\omega (L-0.5)}]利用差分得 ( 2j\sin(\omega(n-0.5)) )因此[\boxed{H_r(\omega) \sum_{n1}^{L} d(n) \sin\left(\omega (n-0.5)\right)}]其中 ( d(n) 2h(L-n) )( n 1,2,\dots,L )。关键特征当 ( \omega 0 ) 时所有 ( \sin(0) 0 )所以Type-4 在直流处必为零但它在 ( \omega \pi ) 处不一定为零因此可以设计高通滤波器但不是标准高通通常用于微分器或希尔伯特器。2. 代码剖析hr_type3.mfunction[Hr,w,c,L]hr_type3(h);%Computes Amplitude response of Type-3 LP FIR filter % 注释中的 LP 不准确应为通用输入长度为奇数的奇对称序列h中心点应为零但代码并未显式检查。输出Hr、w、c正弦展开系数、L。核心代码Mlength(h);L(M-1)/2;c[h(L1:-1:1)];% 注意这里取了包含中心点在内的所有左侧样本但中心点 h(L1) 理论上为 0n[0:L];w[0:500]*pi/500;Hrsin(w*n)*c;分析代码中c包含了h(L1:-1:1)即从中心向左直到第一个样本。由于h(L1)0所以c(1)0正弦级数从n0开始第一项为 0不影响结果。这实际上等价于公式中的c(n)2h(L-n)在索引上的微调MATLAB 索引从 1 开始。例如M5L2h [h0, h1, 0, -h1, -h0]c [0, h1, h0]n[0,1,2]则Hr sin(0)*0 sin(w)*h1 sin(2w)*h0但公式应为2*h1*sin(w) 2*h0*sin(2w)。这里出现了因子 2 的缺失实际上hr_type3.m的写法有误它没有乘以 2也没有正确排除中心零值。对比hr_type1中明确的2*系数hr_type3缺少了这个因子。这是一个潜在的 bug后续在使用时需注意。合理的写法应为c2*h(L:-1:1);% 排除中心点取左半部分并乘 2n[1:L];Hrsin(w*n)*c;但现存的代码仍然可以工作因为所有系数都被缩小了 2 倍形状不变只是幅度缩放。在ampl_ress中调用时若类型识别正确结果仍能保持相对关系。但为了严谨我们应知悉这一点。3. 代码剖析hr_type4.mfunction[Hr,w,d,L]hr_type4(h);%Computes Amplitude response of Type-4 LP FIR filter核心代码Mlength(h);LM/2;d2*[h(L:-1:1)];n[1:L];nn-0.5;w[0:500]*pi/500;Hrsin(w*n)*d;这里d正确乘以了 2且n取半整数与公式完全吻合。所以hr_type4是正确的而hr_type3缺少系数 2但因其中心为零若将中心剔除并乘 2可得到一致结果。实际中我们可以容忍这种幅度缩放因为很多应用中只关心相对形状。4. 自动识别大师ampl_ress.mampl_ress函数是整个幅度响应分析的门面它接收一个冲激响应h自动判断其类型并调用相应的hr_type*输出幅度响应和多项式系数。我们逐段解析其智能逻辑。function[Hr,w,P,L,type]ampl_ress(h)Mlength(h);4.1 奇数长度情况ifrem(M,2)1ifall(abs(h(1:(M-1)/2)-h(M:-1:(M3)/2))1e-8),[Hr,w,P,L]hr_type1(h);type1;elseifall(abs(h(1:(M-1)/2)h(M:-1:(M3)/2))1e-8)h((M1)/2)0,[Hr,w,P,L]hr_type3(h);type3;elsedisp(not a linear-phase filter, check h),return,end首先判断是否为偶对称比较左半部分h(1:(M-1)/2)与右半部分的翻转h(M:-1:(M3)/2)若差值全小于1e-8则为 Type-1。否则判断是否为奇对称比较左半部分与右半部分翻转的负值即相加接近 0并且中心点h((M1)/2)必须为 0满足则为 Type-3。若都不满足则报错退出。4.2 偶数长度情况elseifall(abs(h(1:M/2)-h(M:-1:M/21))1e-8),[Hr,w,P,L]hr_type2(h);type2;elseifall(abs(h(1:M/2)h(M:-1:M/21))1e-8),[Hr,w,P,L]hr_type4(h);type4;elsedisp(not a linear-phase filter, check h),return,endend偶数长度下无需检查中心点直接判断左半与右半翻转的差偶对称 → Type-2或和奇对称 → Type-4。容差1e-8很好地容忍了浮点计算误差。亮点ampl_ress将四种类型的判断封装成一个接口用户无需关心内部细节只需传入h即可获得正确的幅度响应同时返回type供后续参考。5. 使用示例与对比我们构造一个典型的 Type-3 滤波器奇数长奇对称例如一个希尔伯特变换器的冲激响应h[0.1,-0.2,0,0.2,-0.1];% 奇对称中心为零[Hr,w,P,L,type]ampl_ress(h);fprintf(滤波器类型%d\n,type);plot(w/pi,Hr);xlabel(\omega / \pi);ylabel(H_r(\omega));你会看到type 3且Hr在 0 和 π 处均为零数值上接近 0。同样对于 Type-4偶数长奇对称如h [0.1, -0.2, 0.2, -0.1]ampl_ress将正确识别为类型 4并调用hr_type4。6. 奇对称家族的工程用途类型直流响应高频响应典型应用Type-3必为零必为零窄带带通、希尔伯特变换需配合频移、微分器中频段Type-4必为零可不为零宽带希尔伯特变换、微分器高频提升设计时需注意奇对称滤波器的群延迟仍为常数(M-1)/2但其相频特性包含一个固定的 90° 偏移这正是实现正交变换如希尔伯特的关键。7. 资源下载与系列进度本系列所有代码均已打包可点击下方链接免费获取包含全部.m文件及示例脚本下载链接点击下载全套代码截至目前我们已经完成了线性相位 FIR 滤波器幅度响应的全部四种类型分析。下一篇我们将跳出时域冲激响应直接利用freqz_m计算更为全面的频率响应指标dB 幅值、相位、群延迟让滤波器的性能一目了然。下篇预告频率响应进阶——幅值、相位、群延迟一站式计算。我们将揭示freqz_m如何包装 MATLAB 原生函数并提供更加工程化的输出格式同时探讨群延迟对信号失真的影响。思考题若你使用hr_type3.m原版代码计算一个幅度响应发现其数值比理论预期小了一半你能快速定位原因并修改吗欢迎在评论区讨论。