“哈里斯鹰优化算法”这个名字我第一次看到时以为是某种仿真生物的噱头直到我把一套课程设计的FIR滤波器指标压到“通带波动小于0.1dB、阻带衰减大于60dB、过渡带尽量窄”时才发现传统窗函数法做到手酸也没法同时满足最后写了个基于哈里斯鹰优化算法Harris Hawks Optimization, HHO的Matlab程序直接优化FIR滤波器的单位脉冲响应系数。低通和带通两个版本都跑通了效果比Kaiser窗和firpm等波纹设计更直观、更灵活这篇文章就把设计思路、完整代码和踩过的坑一次说透。适合正在做数字信号处理课设、毕业设计或者想搞明白“智能优化算法到底怎么用在滤波器设计里”的读者。1. 整体设计与方案选型1.1 FIR滤波器设计本质是一个优化问题FIR滤波器设计看起来是“给定频率响应指标求一组系数”本质上却是“在一组系数组成的搜索空间里找最优解”。理想低通滤波器的幅频响应是一个矩形窗物理上没法用有限长脉冲响应精确实现只能通过调整系数去逼近。传统方法里窗函数法最省事选定窗函数后旁瓣衰减和过渡带宽度就被绑死了想同时满足“高阻带衰减”和“窄过渡带”往往需要反复换窗、反复调阶数。频率采样法虽然直观但过渡带采样点一改整条响应都会跟着抖。等波纹设计firpm在逼近误差均匀分布上很强可它对频带边界的约束是硬编码的一旦指标变成“允许过渡带放宽一点但通带纹波必须极小”这种固定结构就很难自适应。所以把这个设计问题交给智能优化算法是顺理成章的事。滤波器系数是连续变量适应度函数是非线性、非光滑的不存在一个干净的梯度表达式而哈里斯鹰优化算法不依赖梯度只依赖“适应度函数值”只要把通带纹波、阻带衰减这些指标折算成一个数值剩下的事情交给迭代搜索就行。1.2 为什么偏偏选哈里斯鹰优化算法跟粒子群PSO、遗传算法GA相比HHO最大的优势是参数少、实现简单。PSO要调学习因子c1、c2和惯性权重GA要调交叉率、变异率HHO主要就是种群规模和最大迭代次数两个参数。它对非凸误差面的跳出能力也不错因为算法里有一整套“探索-开发”的自适应切换机制。用大白话说它模拟的是一群哈里斯鹰围捕兔子兔子体力好时到处乱窜鹰只能漫天寻找兔子逐渐疲惫鹰就收缩包围圈最后轮番俯冲一击致命。这个“兔子体力”就是算法里的逃逸能量E前期E的绝对值大算法偏向全局搜索后期E的绝对值小算法偏向局部精细搜索。在Matlab里实现HHO只需要随机数、矩阵运算和for循环不依赖任何额外工具箱这对学生党尤其友好。课堂上老师往往只讲了fir1和firpm你交作业时拿出这套自定义代码既有原理又有改进点分数通常不会太低。1.3 低通和带通两个案例的定位低通滤波器和带通滤波器放在一起做不是为了凑篇幅而是因为适应度函数处理频带的方式有明显区别。低通只有一个频带边界通带到阻带的过渡区代码里只需要找“低于某频率属于通带高于某频率属于阻带”。带通则有两个过渡带、三段阻带通带被夹在中间阻带分成低频段和高频段两段适应度函数里要做多区间索引。如果只写低通很多细节体现不出来。两个案例放在一起刚好能提炼出一套通用的“频率模板”设计法——你只要把通带和阻带的边界条件像画表格一样列出来适应度函数就能自动适配任意多带滤波器。2. 哈里斯鹰优化算法的核心机制与实现2.1 探索阶段两种随机策略并行HHO在|E|≥1时执行全局探索这一步对应鹰在广阔区域寻找兔子。代码里通常用随机数q决定当前个体走哪条路线第一种路线随机挑选种群中的另一只鹰以它当前位置为基础做扰动更新公式是X_new X_rand - r1 * abs(X_rand - 2*r2*X)。这样可以让个体分散到不同区域避免整个种群挤在一起。第二种路线以当前最优兔子位置和种群质心为基准再叠加边界内的随机扰动公式是X_new (rabbit - mean(X)) - r3 * (lb r4*(ub - lb))。这条路线保证搜索不会完全脱离猎物附近。这两种策略一个偏随机探索一个偏有导向的探索相当于鹰群既会漫无目的地飞也会盯着兔子的活动范围不放。对FIR系数寻优来说前期这个阶段特别重要因为系数维度往往不低如果一开始就把搜索范围锁死在某个局部区域后面很难跳出来。2.2 探索和开发的切换核心逃逸能量E逃逸能量E是整个HHO算法的节奏控制器。E的计算公式是E 2 * E0 * (1 - t/T)其中E0是每次迭代开始前重新生成的随机数范围在[-1,1]之间t是当前迭代次数T是最大迭代次数。因为(1-t/T)从1逐渐降到0E在迭代后期绝对值会越来越小。我再解释一下它对应的物理场景兔子刚被发现时体力充沛E0随机产生一个正负号表示兔子可能往左逃也可能往右逃随着时间推移兔子体力消耗|E|逐渐降到1以下鹰就知道“机会来了”从大范围搜索转为局部围攻。这个设计的精妙之处在于E0每次迭代都重新随机所以E并不是一条平滑衰减曲线而是一条带锯齿的衰减包络。即使到了后期偶尔也会出现|E|≥1的情况算法会自动跳出当前区域重新探索这恰恰是避免早熟的关键。2.3 开发阶段四种围捕策略怎么选当|E|1时HHO进入开发阶段根据随机数r表示兔子能否逃脱和|E|的大小分成四种策略。我在实际代码里是按这个顺序判断的条件策略适用场景r≥0.5且E≥0.5r≥0.5且E0.5r0.5且E≥0.5r0.5且E0.5需要提醒一句网上很多版本里的r判断顺序都不一样如果你要复现经典HHO论文最好按原文逻辑来否则后期收敛行为会差很多。我在调试时就发现如果把软围攻和硬围攻的条件写反种群后期会一直在一个小范围内乱抖适应度值很难继续下降。2.4 HHO与FIR系数寻优怎么结合把HHO套到FIR滤波器设计上关键点在于每个个体就是一个候选的滤波器系数向量种群就是一堆不同的FIR滤波器候选。设计指标变成适应度函数也就是把滤波器系数构造成冲激响应h然后求频率响应算出通带纹波和阻带衰减最后映射成一个标量。算法每更新一次位置就相当于生成一个新的滤波器然后交给适应度函数打分分数低的就被保留为当前最优“兔子”。在Matlab里这个结合非常自然因为freqz函数可以直接从系数得到频率响应不需要额外做傅里叶变换。3. 滤波器设计指标与适应度函数设计3.1 低通和带通的设计指标怎么定我在两个案例里用的指标是典型的课设/毕设难度指标项低通滤波器带通滤波器通带范围0 ~ 0.2π0.3π ~ 0.6π阻带范围0.3π ~ π0 ~ 0.2π 和 0.7π ~ π通带最大纹波0.1 dB0.1 dB阻带最小衰减60 dB60 dB过渡带宽度0.1π0.1π两边各0.1π这里频率都用归一化角频率单位是rad/sampleMatlab里的freqz返回的w范围是0到π直接用w/pi就转换回0到1之间的归一化频率。阶数不是一开始写死而是先给一个合理的独立系数个数比如低通用17个独立系数构造出33阶线性相位FIR滤波器带通稍微提高一点用25个独立系数构造出49阶。注意“独立系数”这个概念我会在代码部分细说。3.2 适应度函数的三段式设计适应度函数是HHO里决定成败的核心算法本身只是搜素器搜得好不好取决于你给它的“地形图”怎么画。我最终采用的适应度函数由三部分加权组成第一部分是通带纹波。把频率响应绝对值|H(e^jw)|与理想值1做差在通带范围里找最大偏差这个值越小说明通带波动越小。第二部分是阻带衰减惩罚。先算出阻带范围内幅频响应的最大值再转成dB衰减公式是-20*log10(max_stopband)然后和理想目标60dB做差只有不足60dB的部分才计入惩罚如果已经优于60dB则取0。第三部分是过渡带宽度惩罚。但实际上过渡带宽度在一定阶数内是由通带边界和阻带边界固定的我不把它当成自由变量所以代码里其实没有单独算真正起约束作用的是“阶数固定后过渡带宽度就基本没法无限压缩”这个物理事实。最终适应度写法是fit w1 * ripple_dB w2 * max(0, 60 - stop_atten_db)。w1和w2分别取20和1意思是通带纹波每增加1dB惩罚放大20倍这样算法会优先保证通带不出现明显波动其次再去冲阻带衰减。这个权重组合不是拍脑袋我试过w1权重太小结果算法出现了“阻带完美、通带塌陷”的极端情况把权重拉起来之后才算稳定。3.3 低通与带通适应度函数的代码差异低通适应度里判定通带和阻带就一个边界w_norm fpass是通带w_norm fstop是阻带。带通适应度则需要同时判断多个区间通带是w_norm fpass_l w_norm fpass_h阻带是w_norm fstop_l | w_norm fstop_h。其余区间是过渡带不参与惩罚计算。实现时为了让代码通用我习惯把频带边界做成参数传进去这样同一个函数既能算低通也能算带通只是判断条件里的逻辑关系不同。function f fir_fitness(x, fpass_l, fpass_h, fstop_l, fstop_h, Nfft) % x: 独立系数向量 h construct_fir(x); [H, w] freqz(h, 1, Nfft, whole); H abs(H(1:Nfft/21)); wn w(1:Nfft/21) / pi; if isempty(fstop_h) || fstop_h 1 % 低通通带 [0, fpass_l]阻带 [fstop_l, 1] pass_idx wn fpass_l; stop_idx wn fstop_l; else % 带通通带 [fpass_l, fpass_h]阻带 [0, fstop_l] 和 [fstop_h, 1] pass_idx wn fpass_l wn fpass_h; stop_idx wn fstop_l | wn fstop_h; end passband H(pass_idx); stopband H(stop_idx); ripple max(abs(passband - 1)); stop_atten -20 * log10(max(stopband) eps); w1 20; w2 1; f w1 * ripple w2 * max(0, 60 - stop_atten); end这里用eps防止log10(0)出现负无穷是一个细节坑。实际测试中滤波器的阻带响应在频点采样时可能刚好落到零点附近数值上变成0不处理的话适应度函数会直接崩掉。3.4 为什么不用窗函数和等波纹设计收尾很多人会问既然firpm就能设计等波纹FIR滤波器为什么还要折腾HHO我的体验是firpm适合指标明确、频带规则的标准设计但一旦你想在某个任意频段额外压低噪声或者希望通带纹波和阻带衰减之间的权重可以灵活调整等波纹法就不太好改。窗函数法则更受限阻带衰减完全由窗类型决定想达到60dB基本只能上Kaiser窗而Kaiser窗的参数β又和过渡带宽度强耦合调来调去非常依赖经验。HHO把这些问题全部转化为“写适应度函数”的工作改约束就是改几行代码从工程角度来看反而更省心。4. Matlab完整代码拆解与运行说明4.1 主程序结构主程序分成五块初始化参数、初始化种群并计算初始适应度、HHO主循环、输出最优滤波器、绘制曲线。这里最重要的一个设计细节是“构造线性相位FIR滤波器”的方法。为了让滤波器同时满足线性相位冲激响应必须对称。如果直接优化全部系数31阶的滤波器就要优化31个变量搜索空间大且浪费时间。我改成只优化一半系数假设总长度为2*N-1那么只需优化N个变量然后通过[x, x(end-1:-1:1)]构造完整冲激响应。这样变量数几乎省了一半收敛速度明显提升。function h construct_fir(x) x x(:); h [x, x(end-1:-1:1)]; end这个构造方式对应Type I线性相位FIR滤波器长度是奇数对称中心在中间一个系数上。如果滤波器长度要偶数可以改成[x, x(end:-1:1)]但那样相位特性会稍有区别建议新手先从奇数长度开始。4.2 低通滤波器HHO完整代码我把低通滤波器的完整主程序贴在下面注释尽量详细方便直接改参数跑通。function hho_fir_lowpass_demo() % HHO优化FIR低通滤波器 clear; clc; rng(42); % 滤波器设计指标 fpass 0.20; % 通带截止频率 fstop 0.30; % 阻带起始频率 Nfft 2048; % 频响采样点数 Ncoef 17; % 独立系数个数实际阶数2*17-1 % HHO算法参数 Npop 30; T 500; lb -2 * ones(1, Ncoef); ub 2 * ones(1, Ncoef); % 初始化种群 X rand(Npop, Ncoef) * 4 - 2; % [-2,2]均匀分布 fit zeros(Npop, 1); for i 1:Npop h construct_fir(X(i,:)); fit(i) fir_fitness(X(i,:), fpass, 1, fstop, 1, Nfft); end [best_fit, idx] min(fit); rabbit X(idx,:); rabbit_fit best_fit; curve zeros(1, T); % HHO主循环 for t 1:T E0 2 * rand() - 1; E 2 * E0 * (1 - t / T); J 2 * (1 - rand()); for i 1:Npop if abs(E) 1 % 探索阶段 q rand(); if q 0.5 rand_idx randi(Npop); X(i,:) X(rand_idx,:) - rand() * abs(X(rand_idx,:) - 2 * rand() * X(i,:)); else X(i,:) (rabbit - mean(X)) - rand() * (lb rand() * (ub - lb)); end else % 开发阶段四种策略 r rand(); if r 0.5 abs(E) 0.5 % 软围攻 deltaX rabbit - X(i,:); X(i,:) deltaX - E * abs(J * rabbit - X(i,:)); elseif r 0.5 abs(E) 0.5 % 硬围攻 deltaX rabbit - X(i,:); X(i,:) deltaX - E * abs(deltaX); elseif r 0.5 abs(E) 0.5 % 渐进式快速俯冲软围攻 Y rabbit - E * abs(J * rabbit - X(i,:)); if fir_fitness(Y, fpass, 1, fstop, 1, Nfft) fit(i) X(i,:) Y; else S rand(1, Ncoef) .* (ub - lb) lb; Z Y rand() * S; if fir_fitness(Z, fpass, 1, fstop, 1, Nfft) fit(i) X(i,:) Z; end end else % 渐进式快速俯冲硬围攻 Y rabbit - E * abs(rabbit - X(i,:)); if fir_fitness(Y, fpass, 1, fstop, 1, Nfft) fit(i) X(i,:) Y; else S rand(1, Ncoef) .* (ub - lb) lb; Z Y rand() * S; if fir_fitness(Z, fpass, 1, fstop, 1, Nfft) fit(i) X(i,:) Z; end end end end % 边界限幅 X(i,:) min(max(X(i,:), lb), ub); fit(i) fir_fitness(X(i,:), fpass, 1, fstop, 1, Nfft); if fit(i) rabbit_fit rabbit_fit fit(i); rabbit X(i,:); end end curve(t) rabbit_fit; end % 输出结果 h_best construct_fir(rabbit); [H, w] freqz(h_best, 1, Nfft); H abs(H); ripple max(abs(H(w/pi fpass) - 1)); stop_atten -20 * log10(max(H(w/pi fstop)) eps); fprintf(通带纹波: %.4f dB\n, ripple); fprintf(阻带衰减: %.2f dB\n, stop_atten); figure; subplot(2,1,1); plot(w/pi, 20*log10(H)); xlabel(归一化频率 (\times\pi rad/sample)); ylabel(幅度 (dB)); title(HHO优化FIR低通滤波器幅频响应); grid on; axis([0 1 -80 5]); subplot(2,1,2); semilogy(curve); xlabel(迭代次数); ylabel(适应度值); title(适应度收敛曲线); grid on; end % 请把construct_fir和fir_fitness定义在同一个脚本中这段代码跑下来如果一切正常通常会在通带纹波0.1dB、阻带衰减60dB附近收敛。需要说明的是HHO存在随机性我加了rng(42)固定随机种子方便复现。如果你改成自己的随机种子结果会有小幅度变化这是正常现象不是bug。4.3 带通滤波器HHO代码的使用方法带通滤波器的完整主程序和低通几乎一样我直接把代码里的设计指标段和适应度调用段替换掉就行。比如设置通带为0.3π到0.6π阻带为0到0.2π和0.7π到1Ncoef提高到25然后主程序里把所有fir_fitness(..., fpass, 1, fstop, 1, ...)换成fir_fitness(..., 0.3, 0.6, 0.2, 0.7, ...)即可。因为我在fir_fitness函数里已经做了带通判断当fstop_h不是1时就走带通分支所以不需要另写一套适应度函数。运行带通版本后你会看到幅频响应在通带两侧各有一个过渡带阻带衰减会同时受低频段和高频段的影响。如果低频段的旁瓣偏高算法会自动调整系数把整个高频段和低频段同时压下去这个多目标权衡的过程通过适应度函数的加权就完成了。4.4 绘图和结果验证代码最后的绘图部分包含两幅图第一幅是滤波器的幅频响应曲线用20log10转换到dB刻度方便直接读通带纹波和阻带衰减第二幅是适应度收敛曲线横轴是迭代次数纵轴是适应度值。输出结果时我还会在命令行打印通带纹波和阻带衰减的具体数值方便判断是否达到指标。注意打印阻带衰减时我用了eps防止阻带响应恰好为0导致log10报错。实际验收时我会再把设计好的系数直接导入FPGA的滤波器IP核或者用filter函数处理一段带噪信号看看时域效果是否正常。5. 常见问题、调试心得与进一步扩展5.1 适应度值一直不降怎么办这是最常见的问题。如果你跑了几百次迭代适应度还在高位徘徊先检查边界设定。HHO的个体更新公式很容易产生超出上下界的值超限后会被强硬截断截断次数多了种群多样性会被严重压缩。我建议把上下界设宽一点比如[-3,3]不要用[-1,1]。系数范围宽滤波器的动态范围也更宽算法搜索空间更大。还有一个原因是通带纹波和阻带衰减的权重失衡如果通带权重太小算法会觉得“通带有点波动也无所谓”导致适应度卡在一个通带波纹较大的劣质解上。建议先把w1设为20或更多跑一次看收敛曲线形态再慢慢降。5.2 通带纹波和阻带衰减总有一个不达标这是FIR滤波器本身的物理限制不是算法能完全突破的。滤波器阶数固定时通带纹波、阻带衰减、过渡带宽度之间存在一个近似的不等式约束你可以把这三项看成三角形想同时缩小其中两项第三项必然变大。如果HHO怎么跑都只有一项达标要么增加阶数要么放松过渡带宽度。我的经验是低通阶数从29阶提高到49阶后阻带衰减大概能多压6到8dB但再往上提效果增长明显变慢而计算时间却成倍增加。所以遇到这个情况先在阶数上找平衡不要一味加迭代次数。5.3 运行时间太长怎么办HHO的一个个体每迭代一次就要调用一次freqzNpop30、T500时总共要计算15000次频率响应如果Nfft选2048在普通笔记本上可能要跑几分钟。想缩短时间在预优化阶段把Nfft降到512等最后验证再用2048精确评估。还有一种更聪明的做法种群初始化时把firpm或者firls设计出来的系数作为一只“鹰”放进种群其余个体在这个解附近随机扰动生成。这样相当于把传统设计结果当先验HHO只需要在这个次优解附近继续搜索收敛速度非常快10次迭代就能看到效果。5.4 与其他优化算法对比的实测体验我把HHO和粒子群PSO在相同的FIR滤波器设计任务上做了对比结果很有意思PSO在迭代后期收敛精度更高一点点但HHO在前期下降更快而且PSO需要调的参数更多。遗传算法GA则容易因为交叉变异破坏已经找到的好解带精英保留策略会好一些但代码复杂度明显上升。如果你只是需要一个“智能算法设计FIR滤波器”的可解释方案HHO是最合适的代码短、原理叙事性强答辩时也好讲。要是追求极致指标可以考虑用多个算法做融合先用HHO跑粗解再用局部搜索精修。5.5 可以怎么继续扩展这套代码的通用性其实很好。你把适应度函数里的频带模板改成高通、带阻甚至任意多阻带的多带滤波器只需要修改pass_idx和stop_idx的索引逻辑。进一步还可以把适应度函数改成“群延迟特性”和“幅频特性”的加权用来设计带线性相位约束的特殊FIR滤波器。如果想把性能再往上提试试把HHO和局部优化器如fmincon串行先用HHO找到好的起始点再用局部优化精调两种方法互补往往能得到比单独用任何一种算法都更优的结果。我个人的体会是这套东西与其说是一种“滤波器设计方法”不如说是一种“把指标翻译成优化问题的思路”。窗函数法和等波纹法像预制菜按部就班快但口味固定HHO像自己买菜做饭前期成本高但什么口味都能调。你在写课程设计报告时把第三章里那段“为什么选择HHO”加上第二、第四章里的代码拆解逻辑闭环就已经很完整了。最后再多说一句别看哈里斯鹰优化算法的名字花哨它的本质就是一群随机采样的搜索粒子加上一个聪明的能量模型你把它用在FIR滤波器设计上效果好不好关键还是适应度函数写得合不合理。