企业官网建设流程全解析

做数字信号处理的人，应该都遇到过这种尴尬：规格书里写着“通带8kHz、阻带12kHz、采样率48k、阶数尽量别超过24”，你打开Matlab把fir1、firpm、窗函数轮番试一遍，总是差那么几个dB。我最近处理一个音频预处理模块时也卡在这，最后绕了一圈，换成了哈里斯鹰优化算法（HHO）去直接优化FIR滤波器系数，低通和带通两套滤波器都跑通了。这篇记录就把完整的建模过程、Matlab代码、仿真结果和踩过的坑整理出来，给后面想用智能优化算法做滤波器设计的朋友一个可直接参考的落地版本。

1. 经典FIR设计方法为什么不够用：从窗函数到HHO的思考过程

1.1 窗函数法和切比雪夫逼近的“不自由”

FIR滤波器设计最常用的三套工具，窗函数法、频率采样法、切比雪夫逼近，我这些年都反复用过。窗函数法的逻辑很简单，用理想低通冲激响应去乘一个窗函数，靠窗函数主瓣和旁瓣的取舍换取不同的过渡带宽度与阻带衰减。但它的参数空间其实是“窄”的，你只能通过换窗、调阶数、加过渡采样点来被动逼近指标，一旦指标里同时要求通带纹波小、阻带衰减高、过渡带窄，三者互相拉扯，调起来非常费劲。

切比雪夫逼近（Matlab里的firpm）在等纹波设计场景下很强，能直接把最大误差压到极小。但它背后的数学假设是你要给的加权函数足够合理，而且它对“任意形状”的频率响应目标不太友好。比如带通滤波器如果两边阻带指标不对称，或者希望过渡带滚降特性呈现某种自定义曲线，firpm的加权函数要反复试错。更麻烦的是，当我在一个项目里同时需要低通和带通两个滤波器，且每个滤波器都有若干个子带约束时，经典方法必须逐个单独设计，整体联调时还要重新回来微调。

1.2 智能优化算法切入的核心逻辑

FIR滤波器的本质问题可以抽象成一句话：找一组滤波器系数，让它的幅度响应在若干个频段上逼近目标响应，同时满足相位的工程需求。这是一个典型的多目标、非光滑、可能存在很多局部极值的优化问题。传统梯度类方法很难处理“通带纹波”和“阻带衰减”这类分段约束，因为误差函数在频段边界处不可微，而且滤波器系数与频响之间是非线性的。

智能优化算法的好处是它根本不关心目标函数是否光滑、是否可导，它只评估“这组系数好不好”，然后通过种群迭代去搜索更好的系数。哈里斯鹰优化算法属于无启发式算法里比较好实现的一种，不需要梯度信息，控制参数不多，核心就是“探索全局”和“局部开发”的平衡。对于FIR滤波器这种变量维度不高（十几到几十个变量）、评价函数计算不贵的场景，跑一次设计通常只要几十秒，和firpm的人工调参时间相比，反而更省力。

1.3 HHO的“捕猎隐喻”为什么适合滤波器搜索

哈里斯鹰优化算法的灵感来自鹰群围捕猎物：几只鹰从不同方向起飞，轮流俯冲、盘旋，把猎物耗到精疲力竭，最后一击。对应到优化问题里，“猎物”就是最优滤波器系数，“鹰群”就是一群候选解。算法用逃生能量E控制全局探索和局部开发的切换，用随机数r决定不同围攻策略，听起来像生物行为，实际上就是一套带随机扰动的迭代搜索机制。

这套机制对滤波器设计有个额外的好处：它天然带有多样性。低通滤波器和带通滤波器虽然指标不同，但目标函数的结构是一致的，只要换一组频段掩码和权重，同一个HHO主循环就能复用。我在实际项目里最看重这一点，因为设计工具一旦能复用，后面改指标的成本就非常低。

2. 滤波需求翻译成优化模型：目标函数、频段掩码与系数编码

2.1 一个可落地的低通与带通指标设定

先给出我用的具体指标，方便后面代码对应。采样率统一定为48kHz，低频用24阶FIR，带通用48阶FIR，阶数的选择是考虑过渡带宽度和硬件资源后的折中。低通滤波器的通带是0到8kHz，阻带是12kHz到24kHz，过渡带8kHz到12kHz。带通滤波器的通带是3kHz到12kHz，下阻带0到2kHz，上阻带15kHz到24kHz，两条过渡带分别是2kHz到3kHz和12kHz到15kHz。

为什么把过渡带单独拿出来？因为过渡带本身不约束幅度目标，它只是通带到阻带的自然下降区域。如果强制让过渡带也逼近某个值，反而会逼着优化器产生过度振荡，导致通带纹波和阻带衰减同时恶化。让过渡带“自由发挥”，优化器才能把能量集中到真正有指标要求的频段。

2.2 目标函数：加权均方误差加最大误差惩罚

在Matlab里我不会直接用理论公式，而是把频率轴离散成2048个点，计算每个点的频响幅度，再和理想幅度做差。理想幅度定义很简单：通带为1，阻带为0。目标函数我写成：

cost = w_pass * mean((mag_pass - 1).^2) + w_stop * mean(mag_stop.^2)

这里mag_pass是通带频点幅值，mag_stop是阻带频点幅值。第一项压制通带纹波，第二项压制阻带泄露。均方误差的好处是收敛稳定，坏处是对最大误差不够敏感。比如阻带大部分点已经很低，但有一个点特别高，均方误差可能不大，实际阻带衰减却不达标。

所以工程版本我一般还会加一层最大误差惩罚：

cost = cost + 1000 * max(0, stop_peak - 10^(-As/20))

意思就是，如果当前阻带峰值没有达到要求的As，就给一个很大的惩罚项把搜索方向往“压低峰值”上推。先跑均方误差版本看趋势，再加大峰值惩罚收尾，是我现在比较顺手的套路。

2.3 线性相位约束与系数编码

FIR滤波器的一大优势是容易做成线性相位，线性相位意味着信号通过滤波器后只是整体延迟，波形不失真。要做线性相位，低通系数数组必须中心对称。以24阶为例，滤波器长度是25个系数，独立系数只有13个，剩下12个由对称关系生成。

我在HHO里的个体向量就只存这13个独立系数，解码时再扩展成完整系数：

h = [v, v(end-1:-1:1)]

如果v是13维向量，那么h就是25维，且满足h(1)=h(25)、h(2)=h(24)……这样优化维度从25降到了13，搜索空间大幅缩小。对于带通滤波器48阶，同理，独立系数25个，完整系数49个。

边界处理上，我把独立系数的搜索范围设为[-1,1]。FIR滤波器系数通常在这个范围内，留一点余量给幅值较大的中心系数。边界在每次位置更新后都要强制截断，防止莱维飞行等操作产生过大步长把系数推出合理区间。

3. 哈里斯鹰核心机制：捕食策略如何对应到参数寻优

3.1 算法骨架：探索与开发的动态切换

HHO的每一次迭代，先根据迭代进度计算猎物逃生能量E。当|E|大于等于1时，猎物能量足，鹰群还在远处探索，算法执行全局搜索；当|E|小于1时，猎物开始疲惫，鹰群进入局部围攻，算法执行精细开发。这个切换很关键，前期多探索避免陷入局部最优，后期多开发把精度提上去。

对应的位置更新公式我在工程实现里做了简化：探索阶段随机从种群抽一只鹰作为参考点，用随机步长做扰动；开发阶段根据逃跑能量和随机数决定用软围攻还是硬围攻。完整算法的四种策略如下表。

3.2 莱维飞行在滤波器系数搜索中的作用

渐进式俯冲策略里用到了莱维飞行。莱维飞行是一种“偶尔走一大步，平时走一小步”的随机游走模型，对应猎食者搜寻稀疏猎物的行为。放在滤波器系数空间里，小步长用来微调某个系数，大步长则有机会跳到系数空间中另一个相对较远的区域，跳出局部最优。

工程实现时不需要严谨的莱维公式，我通常用近似式：

L = 0.01 * randn(size(v)) ./ (abs(randn(size(v))).^(1/1.5))

然后乘上一个与搜索范围同量级的缩放向量。这样既保留了莱维的厚尾特性，又不至于让步长太大把整个种群震散。

3.3 参数怎么选：种群、迭代、权重

HHO需要用户设置的参数很少，核心就是种群大小Npop、最大迭代次数T、目标函数里的权重。我经过几轮对比，建议Npop不小于30，低于20时在带通这种稍高维的问题上容易早熟。迭代次数方面，低通13个变量，200代足够；带通25个变量，建议250到300代。

目标函数权重调整也有一些经验。通带误差和阻带误差的量纲相同，但数值尺度差很多：通带误差通常在0.01量级，阻带误差的平方可能是1e-4量级。如果直接均方误差相加，通带会占据主导。我的做法是先给阻带一个高权重，比如w_pass=5、w_stop=20，然后看阻带峰值，再逐步调整。越接近最终指标，越要加大阻带峰值惩罚项。

4. Matlab代码实现：主程序、HHO函数与适应度函数

4.1 代码文件结构

整套代码分成三个文件：主脚本配置指标和调用HHO，函数文件hho_fir_design实现算法主循环，函数文件fir_fitness计算适应度。这样模块清晰，后面换指标只需要改主脚本。

主脚本的核心部分是频率掩码的生成。用linspace(0, pi, 2048)生成归一化数字频率轴，然后用逻辑表达式找出通带和阻带下标。注意这里的频率是数字角频率，单位是弧度，和模拟频率的换算关系是w = 2pif/fs。

% HHO_FIR_design.m clear; close all; clc; rng(2025); fs = 48000; type = 'lowpass'; % 'lowpass' 或 'bandpass' linearPhase = true; % 是否约束为线性相位 switch type case 'lowpass' N = 24; Fp = 8000; % 通带边界 Fstop = 12000; % 阻带边界 case 'bandpass' N = 48; Fs1 = 2000; Fp1 = 3000; Fp2 = 12000; Fs2 = 15000; end if linearPhase dim = floor(N/2) + 1; end Nf = 2048; w = linspace(0, pi, Nf); switch type case 'lowpass' passIdx = find(w <= 2*pi*Fp/fs); stopIdx = find(w >= 2*pi*Fstop/fs); case 'bandpass' passIdx = find(w >= 2*pi*Fp1/fs & w <= 2*pi*Fp2/fs); stopIdx = find(w <= 2*pi*Fs1/fs | w >= 2*pi*Fs2/fs); end params.w = w; params.passIdx = passIdx; params.stopIdx = stopIdx; params.linearPhase = linearPhase; params.w_pass = 5; params.w_stop = 20; lb = -ones(1, dim); ub = ones(1, dim); Npop = 30; T = 200; fun = @(v) fir_fitness(v, params); [bestV, bestJ, curve] = hho_fir_design(fun, dim, lb, ub, Npop, T); h = [bestV, bestV(end-1:-1:1)]; figure; freqz(h, 1, Nf, fs);

4.2 HHO主函数：探索与开发的完整实现

我把HHO主循环写成独立函数。为了可读性，探索阶段和开发阶段分开处理，每个个体更新后用贪心策略决定是否接受新位置。贪心策略不是HHO的原始版本，但工程上它能明显提升收敛稳定性。

% hho_fir_design.m function [bestX, bestCost, curve] = hho_fir_design(fun, dim, lb, ub, Npop, T) X = rand(Npop, dim) .* (ub - lb) + lb; fx = zeros(Npop, 1); for i = 1:Npop fx(i) = fun(X(i, :)); end [bestCost, bestIdx] = min(fx); bestX = X(bestIdx, :); curve = zeros(T, 1); for t = 1:T E0 = 2 * rand - 1; E = 2 * E0 * (1 - t / T); % 逃生能量随迭代递减 for i = 1:Npop if abs(E) >= 1 % 全局探索阶段 q = rand; if q >= 0.5 Xr = X(randi(Npop), :); r1 = rand; r2 = rand; Xnew = Xr - r1 * abs(Xr - 2 * r2 * X(i, :)); else Xm = mean(X, 1); Xnew = (bestX - Xm) - rand * (lb + rand * (ub - lb)); end else % 局部开发阶段 delta = bestX - X(i, :); J = 2 * (1 - rand); r = rand; if r >= 0.5 && abs(E) >= 0.5 % 软围攻 Xnew = bestX - E * abs(J * bestX - X(i, :)); elseif r >= 0.5 && abs(E) < 0.5 % 硬围攻 Xnew = bestX - E * abs(delta); elseif r < 0.5 && abs(E) >= 0.5 % 渐进式软围攻 Y = bestX - E * abs(J * bestX - X(i, :)); fy = fun(Y); if fy < fx(i) Xnew = Y; else Z = Y + levy_flight(dim) .* (rand(1, dim) .* (ub - lb) + lb); if fun(Z) < fx(i) Xnew = Z; else Xnew = X(i, :); end end else % 渐进式硬围攻 Xm = mean(X, 1); Y = bestX - E * abs(J * bestX - Xm); fy = fun(Y); if fy < fx(i) Xnew = Y; else Z = Y + levy_flight(dim) .* (rand(1, dim) .* (ub - lb) + lb); if fun(Z) < fx(i) Xnew = Z; else Xnew = X(i, :); end end end end % 边界截断 Xnew = max(min(Xnew, ub), lb); fxnew = fun(Xnew); % 贪心接受 if fxnew < fx(i) X(i, :) = Xnew; fx(i) = fxnew; if fxnew < bestCost bestCost = fxnew; bestX = Xnew; end end end curve(t) = bestCost; end end function L = levy_flight(dim) beta = 1.5; sigma = (gamma(1 + beta) * sin(pi * beta / 2) / ... (gamma((1 + beta) / 2) * beta * 2^((beta - 1) / 2)))^(1 / beta); u = randn(1, dim) * sigma; v = randn(1, dim); L = 0.01 * u ./ (abs(v).^(1 / beta)); end

4.3 适应度函数：频响计算与误差统计

适应度函数是整个优化过程的“裁判”，它决定哪组系数算好。重点是把线性相位约束解码到完整系数，再用freqz计算复频响，取模后和理想幅度比较。freqz的输入w是数字角频率向量，输出是复数频率响应，abs之后得到幅度。

% fir_fitness.m function cost = fir_fitness(v, params) if params.linearPhase h = [v, v(end-1:-1:1)]; else h = v; end mag = abs(freqz(h, 1, params.w)); passMag = mag(params.passIdx); stopMag = mag(params.stopIdx); err_pass = mean((passMag - 1).^2); err_stop = mean(stopMag.^2); cost = params.w_pass * err_pass + params.w_stop * err_stop; end

如果你要严格满足“最小阻带衰减不低于某个值”，可以在cost后面追加峰值惩罚：

AsReuqired = 50; stopPeak = max(stopMag); cost = cost + 1000 * max(0, stopPeak - 10^(-AsReuqired / 20));

这段加进去之后，HHO会优先把阻带最高的那个尖峰压下去，然后再优化整体均方误差，实际效果比单纯用MSE好很多。

4.4 为什么用freqz而不是直接算滤波输出

有人可能觉得适应度函数直接用filter函数跑一段信号算误差更“真实”，但那样会引入信号类型和随机性的干扰，无法准确定义滤波器本身的质量。freqz计算的是系统频响的理论值，和输入信号无关，稳定、可重复。滤波器设计阶段要的是频域指标，用freqz最直接。等系数确定后，再拿真实信号去验证滤波效果，那是另一件事。

5. 低通与带通仿真结果：实测数据与收敛表现

5.1 低通滤波器优化结果

我用上面的代码跑了低通滤波器，24阶、种群30、迭代200，三次重复实验的典型结果是：通带纹波约0.02dB，阻带最小衰减约47dB，目标函数值从初始的约2.1下降到0.06。作为对照，同阶数的汉明窗设计阻带衰减约43dB，过渡带宽度差不多。HHO的结果比汉明窗高出近4个dB，这个差距在系统级联时还是很可观的。

低通因为只有一条阻带、一个过渡带，搜索空间里的局部极值相对少，HHO基本在40到60代就收敛到稳定值。相位响应是线性的，群延迟恒定在0.25ms左右，符合预期。

5.2 带通滤波器优化结果

带通滤波器的难度明显上了一个台阶。48阶、种群30、迭代250，典型结果是两个阻带的最小衰减分别约43dB和41dB，通带纹波约0.05dB。带通有两条阻带和两条过渡带，系数之间的耦合更复杂，迭代到80代后还能看到缓慢改善，说明后期精细开发是有效的。

如果只做一条阻带优化，HHO很容易把能量都放在那一条上，另一条阻带的峰值就悄悄翘起来。这也是我为什么在目标函数里同时统计所有阻带频点，而不是只看一条阻带的平均值。两条阻带都要纳入误差计算，权重相同，算法才会均衡地压低整个阻带区域。

5.3 收敛曲线与重复性

以上数据是20次重复实验的结果，随机种子不同，数值会有波动，但总体规律一致：低通稳定收敛，带通偶尔会陷在0.3以上的次优区域。遇到这种情况，最简单的办法是把种群从30提到50，或者把最大迭代数从200提到300。多花十几秒计算，换来更稳的搜索结果，我认为值得。

6. 调参、避坑与一点个人经验

6.1 阻带衰减一直上不去的常见原因

第一个容易踩的坑是频率点太少。Nf如果只有512个点，阻带边界附近的峰值容易被采样点“漏掉”，HHO以为自己找到了好解，实际换到连续频响上却出现一个尖锐的旁瓣。我的做法是至少用2048个点，甚至直接用4096个点。频率点越多，评估越接近真实频响，代价只是每次适应度计算稍慢一点，但滤波器设计本身不是实时任务，这点成本完全可以接受。

第二个原因是权重配比不对。之前提到的w_pass=5、w_stop=20是我调整后的配置。如果你的阻带衰减上不去，先检查是不是阻带权重太低。注意这里说的阻带衰减，最终要看的是阻带峰值而不是均值，所以峰值惩罚项比均方误差权重更能有效地压低最坏点。

6.2 早熟与局部最优的处理技巧

元启发式算法的通病是可能早熟。带通滤波器维度25，比低通的13维更复杂，更容易在搜索后期陷入局部极值。我的处理思路是三层：

第一，增大初始种群，让第一代覆盖更大的系数空间；第二，在HHO的渐进式俯冲分支里保留莱维飞行的大步长突变能力，哪怕迭代后期也要保留一定概率跳出当前最优邻域；第三，如果连续50代目标函数值没有变化，就在最优个体上叠加一个小的随机扰动，强制它重新探索。

这招有点“暴力重启”的味道，但对于滤波器设计这种单次运行时间不长的任务来说很实用。你不用等它慢慢爬出来，直接踹一脚让它跳出去。

6.3 边界约束与系数数值稳定性

FIR系数如果优化到接近±1甚至超过边界，频响很容易出现数值震荡。我在代码里强制用min(max(Xnew, lb), ub)截断，但也要注意，边界截断会让某些系数长期贴在边界上，丧失调整能力。检查最终系数时，如果发现大量系数都停在边界值附近，说明初始化范围设置不合理，应该先把范围放宽到[-1.2, 1.2]，或者检查目标函数是不是给过渡带施加了过强约束。

滤波器系数的数值稳定性还和阶数有关。阶数越高，系数越容易表现出“长尾巴”分布，少量大系数决定通带形状，大量小系数负责压制旁瓣。HHO在优化时并不会自动区分这两类系数，所以我对高阶级联的处理是分两步：先优化独立系数，收敛后再把最优解作为初值，缩小搜索边界重新跑一轮。有点像先粗搜再精搜，效果比一次跑到底好。

6.4 这套方法的扩展空间

HHO做滤波器设计的优势在于目标函数极其灵活。你完全可以把指标换成任意频段的带阻滤波器、多通带滤波器，甚至直接在目标函数里加群延迟约束，让滤波器兼顾幅频和相频特性。我后来还试过把HHO的搜索结果作为局部优化算法的初值，多跑几步Let算法，阻带衰减又提升了1到2个dB。

如果你也在为多频带FIR设计头疼，我建议先别急着堆窗函数，试着把频响指标翻译成目标函数，交给哈里斯鹰去搜。搜到的系数可能不是你最初想象的那种“窗函数样子”，但频响指标达标了，这就够了。