免费获取学习方案
ARTICLE DETAIL

资讯详情

深耕编程基础知识与建站技术分享的一线实战洞察。

MATLAB FIR带阻滤波器设计:Kaiser窗参数与实现

MATLAB FIR带阻滤波器设计:Kaiser窗参数与实现 简介这份文档面向学习数字信号处理、需要在MATLAB中实现滤波器的学生与工程人员围绕长度为N45、阻带衰减AS60dB的FIR带阻滤波器展开设计讲解。内容以凯塞-贝塞尔窗函数法为主线说明参数beta如何影响主瓣宽度、旁瓣大小与过渡带宽度并给出Beta0.1102*(As-8.7)的计算依据同时提供freqz.m、ideal_lp、bsfilter.m等源程序演示频率响应、相位响应与群延时的求解过程以及理想冲激响应、窗函数和实际冲激响应的绘图对比。资源包为单个doc文件约60KB结构紧凑便于直接查阅与复现。目前已有1376人学习下载适合希望掌握窗函数选型、滤波器性能分析与MATLAB实现思路的读者参考。1. 从一段能跑的 FIR 带阻代码说起很多人在 MATLAB 里做数字滤波器设计第一步就卡在“窗函数参数到底怎么定”。这份matlab设计FIR带阻滤波器.doc给了一个很具体的答案N45、阻带衰减 As60dB、用凯塞—贝塞尔窗Kaiser 窗来逼近。它没有停在公式推导而是直接给了bsfilter.m、ideal_lp.m、freqz.m三个可运行片段把理想冲激响应、窗函数、实际冲激响应和幅度响应四张图一次性画出来。带阻滤波器要干的事很明确把某一段频率这里是 π/3 到 2π/3压下去其余频段尽量原样保留。FIR 的好处是线性相位、天然稳定代价是要用足够长的阶数去换过渡带陡度。这份文档的价值在于它把“As 决定 beta、beta 决定窗形、窗形决定过渡带”这条链路用代码串了起来适合正在做课程设计、信号处理大作业或者第一次用 MATLAB 手写 FIR 的人照着复现。2. Kaiser 窗参数与带阻滤波器的理论映射2.1 为什么带阻要用“低通相减”来构造FIR 带阻没有现成的闭式公式常见做法是用频域拼接一个截止在 wc1 的低通加上一个截止在 π 的全通其实就是 δ 函数再减去一个截止在 wc2 的低通。文档里这一行是关键bd ideal_lp(wc1,N) ideal_lp(pi,N) - ideal_lp(wc2,N);逻辑上ideal_lp(pi,N)在 n0 处为 1、其余为 0相当于全通ideal_lp(wc1,N)保留 0~wc1减去ideal_lp(wc2,N)就把 wc1~wc2 这段挖掉了。这样得到的bd是理想带阻的冲激响应但它无限长且非因果必须加窗截断。2.2 Kaiser 窗的 beta 与 As 的定量关系Kaiser 窗的可调性全在 beta 上。文档给出的经验式是Beta 0.1102*(As-8.7);这个式子只在 As 50dB 时成立。As 在 21~50dB 区间要用0.5842*(As-21)^0.4 0.07886*(As-21)As 21dB 时 beta 取 0。很多人直接套 0.1102 那一行结果 As40 时过渡带偏宽就是没分段。参数含义本例取值影响N滤波器长度阶数145越大过渡带越窄计算量越大As阻带最小衰减60 dB决定 beta进而决定旁瓣betaKaiser 窗形状参数0.1102*(60-8.7)≈5.65越大旁瓣越低、主瓣越宽wc1下通带截止π/3阻带左边界wc2上通带截止2π/3阻带右边界提示N45 是奇数alpha(N-1)/222为整数ideal_lp里mn-alphaeps加 eps 是为了避免 nalpha 时除零。如果换成偶数 N群延迟会落在半整数点线性相位仍然成立但绘图时 stem 的对称中心会偏移半个点。2.3 过渡带宽度与 N 的取舍Kaiser 窗的过渡带近似满足Δw ≈ (As-8)/(2.285*N)。把 As60、N45 代进去Δw≈0.507 rad换算成归一化频率约 0.08。也就是说 wc1 和 wc2 之间至少要留出这么宽否则阻带边缘会翘起来。文档选 π/3 和 2π/3间隔 π/3≈1.047 rad远大于 0.507所以过渡带是够用的。如果想把阻带压得更窄要么加大 N要么降低 As 要求二者不能同时满足。3. 三个子程序的 MATLAB 实现与逐行拆解3.1 ideal_lp理想低通的 sinc 生成function hd ideal_lp(wc, N) % 生成截止频率为 wc、长度为 N 的理想低通冲激响应 alpha (N-1)/2; % 滤波器中心保证线性相位 n 0:N-1; m n - alpha eps; % 加 eps 防止 m0 时除零 hd sin(wc*m) ./ (pi*m); % sinc 形式pi*m 是归一化 end这里wc用弧度制sin(wc*m)/(pi*m)就是理想低通的标准形式。eps是 MATLAB 的最小浮点数加在分母上让 nalpha 那一点不报 NaN实际值趋近 wc/pi。调用时ideal_lp(pi,N)会得到 nalpha 处为 1、其余接近 0 的序列正好当全通用。3.2 freqz.m频率响应与群延迟的封装function [db,mag,pha,grd,w] freqz_m(b,a) % db 相对振幅(dB)mag 绝对振幅pha 相位grd 群延迟w 频率采样点 [H,w] freqz(b,a,1000,whole); % 取整个单位圆 1000 点 H H(1:1:501); % 只取 0~pi 半圈 w w(1:1:501); mag abs(H); db 20*log10((mageps)/max(mag)); % 归一化到 0dB 峰值 pha angle(H); grd grpdelay(b,a,w); % 群延迟衡量相位线性度 endfreqz的whole参数让采样覆盖 0~2π取前 501 点就是 0~π。db里除以max(mag)是为了把通带峰值归一到 0dB看阻带衰减时直接读负值即可。grpdelay对 FIR 来说理论上应恒等于 (N-1)/222如果画出来波动大说明系数算错了。3.3 bsfilter.m主程序串起全流程N 45; As 60; n 0:N-1; beta 0.1102*(As-8.7); w_kai kaiser(N, beta); % 注意转置成列向量 wc1 pi/3; wc2 2*pi/3; bd ideal_lp(wc1,N) ideal_lp(pi,N) - ideal_lp(wc2,N); h bd .* w_kai; % 加窗截断 [db,mag,pha,grd,w] freqz_m(h,[1]);kaiser(N,beta)返回列向量ideal_lp返回行向量所以w_kai要转置否则.*会因维度不匹配报错——这是新手最常见的坑。freqz_m(h,[1])里 a1 表示 FIR分母只有常数项。画图部分用subplot(2,2,1)到(2,2,4)排四张图stem画离散冲激响应plot画连续幅度曲线。axis([0,1,-80,10])把横轴限制在归一化频率 0~1纵轴 -80~10dB方便看 60dB 的阻带深度是否达标。4. 复现时的排错与性能验证4.1 常见报错与对应修法现象原因修法Matrix dimensions must agreew_kai与bd行列不一致给kaiser结果加转置阻带衰减只有 40dBbeta 用了 0.1102 公式但 As50换分段公式算 beta群延迟曲线波动系数非对称或 N 取偶数未调整检查ideal_lp的 alphafreqz_m未定义文件名与函数名不一致存为freqz_m.m调用同名4.2 用数据验证 60dB 是否真的达到画图只能看个大概要确认指标得直接读db数组% 找阻带区间对应的索引 idx find(w/pi 1/3 w/pi 2/3); max_stop max(db(idx)); % 阻带内最大 dB 值 fprintf(阻带最大衰减: %.2f dB\n, max_stop);如果max_stop在 -60 以下说明设计达标如果在 -55 左右说明 N 偏小或 beta 偏小把 N 加到 51 或 61 再试。注意db已经归一化通带峰值是 0dB所以阻带值直接就是衰减量。4.3 过渡带与阶数的边界前面算过 Δw≈0.507实际过渡带会略宽于理论值。如果发现 wc1 附近通带就开始掉说明过渡带侵入了通带需要把 wc1 往左挪或加大 N。FIR 带阻没有免费午餐As、过渡带宽度、N 三者固定两个第三个必然被约束。课程设计里常见的要求是 As60dB、过渡带 0.1π那 N 至少要 100 以上45 是压不住的。5. 从脚本到可复用函数参数化与批量验证把bsfilter.m改造成函数才能快速扫参数function [h,db,w] design_bs(N, As, wc1, wc2) % 参数化 FIR 带阻设计返回系数与幅度响应 if As 50 beta 0.1102*(As-8.7); elseif As 21 beta 0.5842*(As-21)^0.4 0.07886*(As-21); else beta 0; end n 0:N-1; w_kai kaiser(N, beta); bd ideal_lp(wc1,N) ideal_lp(pi,N) - ideal_lp(wc2,N); h bd .* w_kai; [db,~,~,~,w] freqz_m(h,[1]); end调用时一行就能换参数[h,db,w]design_bs(61,60,pi/3,2*pi/3);。想批量看 N 对阻带的影响套个循环把max(db(idx))打出来比一张张画图快得多。实际工程里还会把h导出成系数文件给 FPGA 或 DSP 用注意 MATLAB 的freqz是归一化频率移植时要按采样率换算成实际 Hz。这套脚本跑通一次后面改指标就是改几个数字的事。本文还有配套的精品资源点击获取
返回列表