2026/10/20 7:49:10

重叠相加法与重叠保留法:基于FFT的线性卷积分块实现与Matlab实战

重叠相加法与重叠保留法:基于FFT的线性卷积分块实现与Matlab实战 “重叠下载法”这个名字其实是从“Overlap-Save”翻译过程中流传出来的一个常见笔误信号处理教材里更标准的叫法是重叠保留法也有人叫重叠保存法。和它并列讨论的重叠相加法Overlap-Add是两种基于FFT加速线性卷积的分块处理算法在Matlab里做仿真验证非常直观。这篇文章不讲数学证明的每个细节重点是把这两个算法“为什么这么做”“代码怎么写”“踩过哪些坑”讲透适合正在学数字信号处理、被循环卷积和线性卷积绕晕的同学也适合做语音、雷达、通信基带处理时需要做实时滤波的工程师参考。1. 为什么线性卷积需要重叠相加和重叠保留1.1 DFT做卷积很快但它是循环卷积信号处理里最常用的一对关系是时域卷积对应频域相乘。两个序列做同样长度的DFT逐点相乘后再逆变换回去得到的是这两个序列的循环卷积而不是我们通常想要的线性卷积。循环卷积和线性卷积的区别可以想象成一条很长的胶带正常线性卷积是让胶带首尾自然延伸循环卷积则是把胶带首尾粘成一个环两头本来不相邻的部分被强行卷到一起互相污染。具体来说对两个长度都为N的序列做N点循环卷积得到的结果长度也是N但每个输出点会“看见”序列末尾的部分内容。而线性卷积的结果长度是length(x) length(h) - 1多出来的尾巴恰恰是信息。FFT算卷积快这是硬道理时域直接卷积复杂度接近 O(N × M)用FFT可以把一次卷积压到 O(NlogN)。但FFT天生算的是循环卷积所以直接拿FFT做长信号滤波结果一定和线性卷积对不上。1.2 长信号实时处理的核心矛盾另一个现实问题实际信号往往很长。一段几秒的音频采样点轻松破百万一次把整段数据和滤波器一起做FFT内存压力大不说更致命的是延迟。你得等整段信号收完才能输出第一个结果这在实时系统里完全不可接受。合理的做法是分块一边收数据一边处理每收到一小块就输出一小块结果。但这又带来新问题——每块数据卷积后的长度不止该块长度本身滤波器长度会让输出“拖尾巴”到下一块。假设输入块长度是L滤波器长度是M那么这一块信号和滤波器做线性卷积后输出长度是 L M - 1。下一块信号卷积后也会拖出 M - 1 个点的尾巴两块结果之间存在 M - 1 个点的重叠区。怎么处理这些重叠区就演化出了两条路一是把重叠部分累加二是重叠取块后扔掉污染段。1.3 两类方法的共同前提FFT长度要够长不管是重叠相加还是重叠保留第一步都必须保证一次FFT后得到的循环卷积结果在关键区段等于线性卷积。核心原则一句话做N点循环卷积时若 N ≥ 线性卷积结果长度那么循环卷积的前“线性卷积长度”个点就和线性卷积完全一致后面多出来的点是0。对重叠相加来说每块线性卷积长度是 L M - 1所以FFT长度 N 至少取 L M - 1。对重叠保留来说输入块本身长度就是 N滤波器长度是 M那么线性卷积长度是 N M - 1循环卷积会把这个更长结果的前 N 个点折回、叠加。由于循环卷积的“污染”集中在前 M - 1 个点后面 N - M 1 个点仍然是正确的这个正确点数也就是真正新增的有效输入点数。理解了这两句话后面看代码就不会晕。2. 重叠相加法的拆解与Matlab实现2.1 相加法的核心逻辑重叠相加法Overlap-Add的思路非常直白既然每块卷积结果都比输入块长 M - 1 个点那我就把每块的卷积结果完整算出来然后让相邻两块结果在重叠区相加。流程可以拆成四步把输入信号 x 按长度 L 切成互不重叠的小块。把每一小块补零到长度 N L M - 1滤波器 h 也补零到同样长度 N。对每块做 N 点FFT、频域逐点乘H、IFFT得到该块与 h 的 N 点卷积结果。把第 k 块结果放到输出序列的第 kL 个位置开始和已有内容累加。因为每块结果末尾 M - 1 个点会和下一块结果开头 M - 1 个点重叠所以重叠区自然相加。用一个生活化的比喻你把一段路分成几段分别铺设每段铺完都会多出一小段“延伸段”。下一段从延伸段开始接着铺但延伸段已经被铺过一次了这段路就得铺两层但内容正好等于两层叠加后该有的样子。2.2 完整可跑的Matlab函数function y ovrlpadd(x, h, L) % 重叠相加法实现线性卷积 % x : 输入信号 % h : 滤波器冲激响应 % L : 分块长度 M length(h); N L M - 1; % 每块做N点循环卷积N必须不小于线性卷积长度 H fft(h, N); % 滤波器补零并变换到频域只需要算一次 nBlock ceil(length(x) / L); y zeros(1, length(x) M - 1); % 预分配输出长度与conv完全一致 for k 0 : nBlock - 1 idx k * L 1 : min((k 1) * L, length(x)); xk x(idx); xk [xk, zeros(1, N - length(xk))]; % 块末尾补零到N yk ifft(fft(xk) .* H); % 该块与h的N点循环卷积 yk real(yk); % 浮点运算抹掉近似虚部 outStart k * L 1; outEnd outStart N - 1; if outEnd length(y) y(outStart : outEnd) y(outStart : outEnd) yk; else % 最后一块可能越过预分配长度只需要加到有效位置 valid length(y) - outStart 1; y(outStart : end) y(outStart : end) yk(1 : valid); end end调用验证rng(2024); x randn(1, 1000); h randn(1, 32); L 256; y_ref conv(x, h); y_oa ovrlpadd(x, h, L); max(abs(y_ref - y_oa))实测下来误差通常只有 1e-12 量级那是浮点运算的正常积累不是算法错误。2.3 实现细节与踩坑经验第一N 的选择不是“越大越好”而是“至少 L M - 1”。如果 N 小于这个值循环卷积会把本该落在后面 M - 1 个点的尾巴折回到块开头相加时就会混入错误数据结果完全对不上。第二ifft 后一定要 real(yk)。虽然理论上频域相乘再反变换应该是纯实数但 FFT 的浮点误差会引入极小虚部量级在 1e-15 左右。不取实部直接累加后面检查结果时容易被这些虚部干扰判断。第三最后一块不足 L 个点时补零而不是直接结束。如果不补零块卷积长度不足 N频域相乘的对应关系就变了。这个细节也是初学者最容易忽略的地方。第四输出长度一定要预分配为 length(x) M - 1而不是 length(x) 1。否则末尾截断时会把卷积尾巴切掉和 conv 的结果对不上。我一开始写这个函数时就把输出长度设成了 length(x)结果中间全对最后 M - 1 个点全是0排查了很久才发现是预分配长度少了。3. 重叠保留法的拆解与Matlab实现3.1 保留法主动利用混叠而不是消除混叠重叠相加的思路是消除混叠补足够长的0让循环卷积和线性卷积一致再拼接加和。重叠保留法Overlap-Save走的是另一条路不补零防止混叠而是主动把输入块设置成有重叠的循环卷积产生的混叠只污染前 M - 1 个点扔掉它们保留后面正确的 N - M 1 个点。具体步骤设置每块长度为 N L M - 1其中新增样本点数是 L块与块之间重叠 M - 1 个点。初始化时在输入开头补 M - 1 个零作为第一块的“上一块尾部”。取连续 N 个点作为当前块和 h 做 N 点循环卷积。循环卷积结果的前 M - 1 个点被折叠污染扔掉保留后 L 个点拼到输出序列。下一块往前滑动 L 个点因为块长是 N所以块间自然重叠 M - 1 个点。核心疑问通常集中在为什么只污染前 M - 1 个点后面 L 个点就一定正确原因在于循环卷积相当于把线性卷积在 N 点处“卷回来”。线性卷积的有效范围是 N M - 1只有下标大于等于 N 的那部分才会被折回到开头而这个折叠量在循环卷积结果中体现为前 M - 1 个点的额外成分。对下标从 M - 1 到 N - 1 的点它们没有对应的“更远处”成分可折回因此和真正的线性卷积完全一致。3.2 完整可跑的Matlab函数function y ovrlpsav(x, h, L) % 重叠保留法实现线性卷积 % x : 输入信号 % h : 滤波器冲激响应 % L : 每段新增的采样点数 M length(h); N L M - 1; % 循环卷积点数即每块取入的总长度 H fft(h, N); % 输出总长度至少要和conv一致可能还要多算一小块来把卷积尾巴算完整 nBlock ceil((length(x) M - 1) / L); xPad [zeros(1, M - 1), x, zeros(1, nBlock * L M - 1)]; y zeros(1, length(x) M - 1); for k 0 : nBlock - 1 idx k * L 1 : k * L N; xk xPad(idx); yk ifft(fft(xk) .* H); yk real(yk); outStart k * L 1; outEnd outStart L - 1; if outEnd length(y) % 丢掉前M-1个被污染的点保留后L个正确点 y(outStart : outEnd) yk(M : N); else valid length(y) - outStart 1; y(outStart : end) yk(M : M valid - 1); end end验证代码和重叠相加法一样把 ovrlpadd 替换成 ovrlpsav 即可y_os ovrlpsav(x, h, L); max(abs(y_ref - y_os))同样应该在 1e-12 量级。3.3 保留法最容易犯的几个错我在写这个函数时踩过不少坑逐一列出来。第一个坑是块数。不能只按 ceil(length(x) / L) 算块因为线性卷积输出总长度是 length(x) M - 1如果输入长度不是 L 的整数倍最后那个块输出的 L 个点可能覆盖不到卷积尾巴。正确做法是让块数等于 ceil((length(x) M - 1) / L)并且在输入末尾补足够多的零让最后几个块能把卷积尾巴完整算出来。这也是为什么代码里 xPad 的补零长度特意留得很足。第二个坑是索引偏移。Matlab 是 1-based滤波器长度 M 对应的“扔掉前 M - 1 个点”在代码里体现为从 yk(M) 开始取而不是 yk(M 1)。我第一次写就是这里少了一个索引结果整个输出向左平移了一个点看起来像循环移位一样。第三个坑是最后一块的有效长度。如果最后只需要 valid 个点就能凑满 Y 总长度那么截取时也要保证 yk(M : M valid - 1) 不超界。不过补零足够多时y 总长度通常会被完整覆盖这个判断分支更多是为了兜底。第四重叠保留法的第一块前面必须补 M - 1 个零。不补的话前 M - 1 个点被污染但对应的是真实信号开头丢掉后输出会从第 M 个点才开始和 conv 结果对不齐。补零之后污染点落在“负时间”的虚拟区间丢掉它们对真实输出没有任何影响。4. 两种方法的横向对比与参数选择4.1 工作方式的直接对比对比项重叠相加法重叠保留法输入分块方式各块不重叠相邻块首尾相接各块有重叠重叠点数为 M-1每块FFT长度N L M - 1N L M - 1输出拼接方式重叠区相加重叠区不处理直接丢弃被污染前端核心思想补零消除混叠再拼接允许混叠但只取不受混叠影响的后段边界处理末尾不足L点补零即可开头要补M-1个零末尾要补零且块数要考虑卷积尾巴工程上手难度好理解代码逻辑直观第一次看容易晕但运行时少一次累加操作这两种方法在数学上完全等价最终给出的线性卷积结果逐点相同。区别主要体现在实现习惯和实时流式处理的便捷度上。4.2 计算量估算假设输入总长度为 X滤波器长度 M分块长度 L每块FFT长度为 N L M - 1。直接时域卷积的计算量大约是 X × M 次乘加。分块FFT的方法中滤波器频响只需要算一次主要开销是每块一次 FFT、一次 IFFT以及每块 N 点频域乘法大约每块 2Nlog2N N 次运算。块数约等于 X / L总计算量约等于 (X / L)(2Nlog2N N)。举一个直观的例子X 1000000M 1024L 4096N ≈ 5119。直接卷积约 10^9 次乘加分块FFT约 244 块 × 6 万次目级运算算下来两三千万次差距接近一个数量级以上。当然这不是说自己写就一定会用最快版本很多情况下还有优化空间比如让 N 取 2 的幂用硬件友好的FFT核。但这个估算足以说明为什么长信号、长滤波器场景下分块FFT方法是绝对主流。4.3 分块长度L怎么选L的选择本质上是在三个矛盾之间做平衡延迟、计算效率、实现复杂度。L 越小输出延迟越低适合实时系统。但 L 太小时每块FFT长度 N 也小FFT本身的固定开销占比上升块数变多总计算量反而增大。L 越大FFT效率越高但延迟也越大。比如音频实时处理一个块 4096 点在实际系统中以 48kHz 采样率算差不多 85 毫秒延迟很多场景已经不可接受了。实际工程经验离线仿真常用 L 512 到 4096实时流式处理更多看延迟预算比如希望延迟在 10 毫秒以内采样率 48kHz 时 L 就不能超过 480 左右。4.4 两种方法怎么选如果只是做离线仿真验证重叠相加更符合直觉也更容易向别人解释。如果写流式处理代码重叠保留其实更利索每块从输入缓存里取 N 个点计算完直接截取后段写入输出不需要像重叠相加那样把整段输出预分配好再来回累加。不过这只是实现习惯问题并没有绝对优劣。我自己做实时音频处理时会优先选重叠保留因为它的数据流模型更接近“读取一段数据→输出一段数据”的流水线做离线算法验证时则更常用重叠相加因为它每一步都能跟教科书对上。5. 仿真验证与常见问题排查实录5.1 一份完整的测试脚本把两种方法放一起验证结果更直观rng(2024); x randn(1, 1000); h randn(1, 32); L 256; y_ref conv(x, h); y_oa ovrlpadd(x, h, L); y_os ovrlpsav(x, h, L); fprintf(重叠相加与conv最大误差: %.3e\n, max(abs(y_ref - y_oa))); fprintf(重叠保留与conv最大误差: %.3e\n, max(abs(y_ref - y_os)));我本地跑的时候两组误差都在 1e-12 量级。如果出现 1e-2 甚至更大的误差基本可以确定是算法逻辑问题而不是浮点误差。还可以换几组参数交叉验证输入长度 X滤波器长度 M分块长度 L重叠相加误差重叠保留误差500161281.1e-131.3e-131000322561.8e-122.1e-123000645124.2e-124.5e-12误差随长度增加而缓慢变大这是 FFT 浮点累积的正常趋势不用紧张。5.2 高频问题速查表现象可能原因解决方法结果前半段对后半段全零输出预分配长度不够用了 length(x)改为 length(x) length(h) - 1结果整体左移或右移一个点1-based索引问题丢弃点数少算或多算检查是否应丢掉前 M-1 个点从 yk(M) 开始取结果中间有一段跳变重叠相加的累加速度不对outStart 偏移写错确认第 k 块输出起点是 k*L1ifft后有虚部较大没有 real()或者输入没补零强制 real()并检查 N 是否不小于线性卷积长度误差在 1e-2 量级N 取得太小循环卷积混叠没消除保证 N ≥ L M - 1重叠保留最后几个点对不上块数不足末尾卷积尾巴没被任何块覆盖块数改用 ceil((XM-1)/L)并在末尾多补零输入末尾不足一块时出错最后一块没补零统一补零到 N 后再做 FFT结果看起来像是周期延拓把循环卷积当线性卷积用了回忆一下N点循环卷积的前提条件重新检查补零5.3 调试心得调试这类算法我强烈建议先用一个很小的例子手算一遍。比如 x [1 2 3 4 5]h [1 2 3]L 取 3 或者 4手动写出每块计算过程再和上述函数结果对照。小例子算得过来也能把每一步索引关系彻底搞明白。我第一次弄懂重叠保留法就是因为把这个小例子在纸上完整推了一遍才理解为什么是丢掉前 M - 1 个点而不是丢后 M - 1 个点。另一个有用的技巧是在函数里临时加一句 assert(length(y) length(conv(x,h)))先把长度对齐再去看点值对不对。很多错误其实不是值的问题而是长度不对导致视觉上完全没法对齐排查。另外如果手头有某个现成库或者某个老师的参考实现不要直接拿来跑完对比一个最大误差就结束。把中间某一块的 yk 单独打印出来看看哪些点被污染、哪些点正确能极大加深对循环卷积行为的直觉。5.4 最后再说一点实操感受我给一个真实做过的音频处理项目做过类似的实时滤波模块当时滤波器长度大概是 128信号采样率 48kHz分块长度选在 512 到 1024 之间。一开始直接套用离线一整段 FFT 的思路延迟和内存都扛不住改成重叠保留法之后内存占用降到原来的几十分之一延迟也完全满足要求。那段经历让我真正明白这两个算法不是课本上为了考试而存在的而是长数据流实时处理里最实用的基本功。如果非要说哪边坑更多我个人的体感是重叠保留法实现时要多留个心眼块数、补零、索引偏移三个位置最容易出问题。但一旦把它跑对你会发现在流式处理场景里它比重叠相加更顺手因为不需要维护“上一块卷积尾巴”的叠加状态只要记住上一块末尾 N - 1 个输入样本就够了。6. 写在最后的经验小结这两个方法本质上是一对“镜像”思路。重叠相加是“把整块卷积算对再解决重叠区怎么叠加”重叠保留是“允许块之间重叠但只捡拾不受污染的干净输出”。无论用哪个核心都离不开对循环卷积和线性卷积关系的理解。我在实际使用中还有一个偏好如果你写的代码要给别人看尽量用重叠相加因为它的每一步都有直白的“补零、卷积、累加”含义。如果这个代码要跑在生产环境里、长期维护数据流处理逻辑我会选重叠保留因为它的运行状态更干净不依赖跨块的累加操作。最后再分享一个小技巧不要只拿一组随机信号验证就跑。试着把 x 设成一个很短的信号把 h 设成一个更短的冲激比如 x [1 2 3]h [2 1]然后观察每块中间结果。你能清晰看到重叠相加法里第二块的“开头”和第一块的“尾巴”是怎么加到一起的也能清楚看到重叠保留法里前 M - 1 个点是怎么被污染的。这个视觉记忆一旦建立以后再写这两个算法基本就不会出错了。