提到 NumPy 实现卷积很多人的第一反应是“这不就是一层循环吗”。可等你真去写会发现边界模式、核翻转、输出尺寸和频域混叠每一处都能把你折腾半天。最近我把手头一个图像预处理工具收敛到纯 NumPy 层面只用数组操作完成了卷积的两种实现一种直接滑窗相乘再相加另一种用 FFT 把空间域卷积转成频域乘法。这篇把两种方法的思路、完整代码和踩过的坑放在一起讲清楚。如果你只是想在项目里临时算一个卷积当然可以直接调现成的图像库函数。但如果想理解卷积到底在做什么、想自定义核、想把算法移植到只有 NumPy 的轻量环境里或者想为后面的深度学习反向传播打底那这两种手写实现就很有价值。文章按“基础概念 → 方法一滑动窗口 → 方法二 FFT → 对比和踩坑”的顺序展开适合读过一点 NumPy、但还没自己动手实现过卷积的开发者。1. 先搞清楚你写的到底是卷积还是互相关1.1 一维到二维的数学约定离散卷积的标准写法是y[n] sum_k x[k] h[n-k]注意这里的h[n-k]意味着卷积核必须翻转一次再做加权求和。也就是说数学意义上严格卷积计算的是“输入与翻转后的核”的滑动相乘。但我们在图像处理和深度学习里常用的卷积层其实一般不翻转核严格说应该叫“互相关”。比如一个 3×3 的拉普拉斯核你写出来是什么样算的时候就是什么样不会先上下翻转再左右翻转。对于对称核翻转不翻转没区别但对于非对称核比如 Sobel 算子到底翻不翻结果方向可能完全不同。我下面两个实现默认都用“深度学习里的卷积”方式核不翻转直接取窗口做点积。这样写出来的结果和 PyTorch 的卷积行为是一致的大家也更容易对照。如果你非要标准数学卷积也很简单每个实现里把核先做一次kernel[::-1, ::-1]再传入即可改造成本几乎为零。1.2 输出尺寸和三种边界模式手写卷积最容易出问题的地方不是乘加计算而是输出到底多大、边界怎么补。常见有三种模式valid不做任何填充核只在图像内部滑动输出尺寸是H-kh1、W-kw1。图像会变小。same输出尺寸保持和输入一样大。做法是在输入四周补一圈 0补多少取决于核大小。通常要求核是奇数比如 3×3 就上下左右各补 1 像素。full输出尺寸是Hkh-1、Wkw-1。这种模式下图像会被放大核每次移动一步直到只露出一像素在对面。不少资料会把“same”和“full”混着说但实际编码时差异很大。建议在动手前先把目标输出尺寸写在注释里否则后面调试时很容易越调越乱。1.3 为什么要自己实现而不是直接调包一句话为了掌握底层也为了减少环境依赖。很多嵌入式或者轻量分析环境里并没有完整的图像处理库但装一个 NumPy 要容易得多。另外图像处理里的高斯模糊、边缘检测、锐化本质上都是小核卷积。如果只用一个 3×3 核直接手写循环可能比反复调库函数还更快更可控。反过来说如果你以后想实现卷积层的反向传播或想理解卷积怎么转成矩阵乘法这篇里提到的滑动窗口和 FFT 思路就是很好的热身。基础不打牢后面写加速代码会非常痛苦。2. 第一种方法滑动窗口法直接算2.1 最直观的三层循环版把卷积“翻译”成人话第一种方法最简单把卷积公式直接翻译成循环。对二维单通道图像用两层循环遍历每个输出像素再对每个窗口做一次加权求和。import numpy as np def conv_direct(img, kernel, modevalid): H, W img.shape kh, kw kernel.shape if mode same: p_h, p_w kh // 2, kw // 2 img_pad np.pad(img, ((p_h, p_h), (p_w, p_w))) return conv_direct(img_pad, kernel, modevalid) if mode full: img_pad np.pad(img, ((kh - 1, kh - 1), (kw - 1, kw - 1))) return conv_direct(img_pad, kernel, modevalid) out_h H - kh 1 out_w W - kw 1 out np.zeros((out_h, out_w), dtypenp.float64) for i in range(out_h): for j in range(out_w): out[i, j] np.sum(img[i:i kh, j:j kw] * kernel) return out这个函数我故意用递归处理了same和full。第一次进来时先对图像做np.pad再把 padding 后的图以valid模式重新调用。这样避免了在每个分支里重复写一堆输出索引逻辑也非常直观。核心在于这一行out[i, j] np.sum(img[i:i kh, j:j kw] * kernel)img[i:ikh, j:jkw]取出以(i, j)为左上角的窗口然后和核逐元素相乘最后累加。这正是“卷积是加权求和”的翻译。这里有几个注意点。第一内层如果再用一层 for 去遍历核元素性能会很差直接用切片加np.sum已经是最简写法。第二输出数组建议显式dtypenp.float64否则输入是 uint8 时结果会被截断。第三modesame的填充方式我用的是kh // 2这默认核是奇数如果是偶数核需要额外讨论前后 padding 怎么分配建议日常使用时尽量选奇数核。2.2 用 sliding_window_view 干掉显式循环上面的三层循环虽然好懂但索引写法不够优雅。NumPy 从 1.20 开始提供了sliding_window_view可以一次把所有窗口都“看”出来。注意它返回的是视图不是复制副本所以内存占用比你想的小得多。from numpy.lib.stride_tricks import sliding_window_view def conv_sliding(img, kernel): kh, kw kernel.shape windows sliding_window_view(img, (kh, kw)) # windows 的形状 (H - kh 1, W - kw 1, kh, kw) out np.tensordot(windows, kernel, axes([2, 3], [0, 1])) return out这段代码理解起来分两步。sliding_window_view把输入切成了无数个和核一样大的窗口windows的后两个维度正好就是窗口内容。然后np.tensordot把后两个维度[2, 3]和核的两个维度[0, 1]做点积相当于“窗口内一一对应相乘再求和”。和手工三层循环对比算术量没有变化还是每个窗口做kh * kw次乘法。但代码很短而且底层利用了步长视图索引开销小很多。在实际项目中我更常用这个版本。2.3 边界模式封装把 valid/same/full 统一起来为了公平对比我封装一个统一入口。这里我仍然推荐奇数核same模式下才不需要写一堆难懂的偏移。def conv_sliding_mode(img, kernel, modevalid): kh, kw kernel.shape if mode valid: return conv_sliding(img, kernel) if mode same: assert kh % 2 1 and kw % 2 1, same模式请使用奇数核 top, left kh // 2, kw // 2 img_pad np.pad(img, ((top, top), (left, left))) return conv_sliding(img_pad, kernel) if mode full: img_pad np.pad(img, ((kh - 1, kh - 1), (kw - 1, kw - 1))) return conv_sliding(img_pad, kernel) raise ValueError(f未知模式: {mode})np.pad里的((top, top), (left, left))表示在第一个维度高的前后各补top第二个维度宽的前后各补left。边界补零是最常用的做法有时也会用常量填充或镜像填充卷积实现里一般先用补零。这里有一个容易踩的小坑不要盲目把same的 padding 设置成(kh-1)//2然后期待kh是奇数就行。如果核尺寸是偶数kh//2和(kh-1)//2结果不同输出会偏移一个像素。我会在最后一节再展开。2.4 一个小实验验证两种滑动窗口写法完全一致写实现时不要直接拿大图去跑先用随机二维数组做一致性测试。比如rng np.random.default_rng(0) img rng.random((32, 32)) kernel np.array([[0, -1, 0], [-1, 4, -1], [0, -1, 0]], dtypenp.float64) out1 conv_direct(img, kernel, modesame) out2 conv_sliding_mode(img, kernel, modesame) print(np.max(np.abs(out1 - out2))) # 应该是 0.0我这里用的是拉普拉斯核也就是一个锐化算子。它对中心点和上下左右邻居做差分能突出图像中的细边缘。测试结果里两个版本逐像素完全一致说明索引和边界处理没有写错。3. 第二种方法FFT 频域乘法3.1 卷积定理一句话的大白话版如果只用一句话描述卷积定理空间域里的卷积等于频率域里的逐元素乘法。更严格地说两个信号的傅里叶变换相乘后再做逆傅里叶变换得到的结果就是这两个信号的卷积。这个定理的价值在于复杂度。直接滑窗法的时间复杂度是O(H * W * kh * kw)图像越大、核越大耗时越明显。FFT 的复杂度大约在O(N log N)当核足够大时频域乘法会快非常多。用生活里的例子理解你有一长串数据想用一个 100 长度的核去滑动求平均。滑窗法每个位置都要做 100 次乘法而 FFT 相当于先把整串数据变成频率成分把核也做同样的变换两者在频域“各频点相乘”最后再重新合回空间域。虽然前后各做了一次变换但整体计算量反而更小。3.2 FFT 卷积的完整实现下面的函数是核心def next_power_of_two(n): return 1 (n - 1).bit_length() def conv_fft(img, kernel, modevalid): H, W img.shape kh, kw kernel.shape # 因为我们想和滑动窗口法的互相关结果一致先翻转核 kernel_flipped kernel[::-1, ::-1] # 线性卷积的 full 输出尺寸 out_h_full H kh - 1 out_w_full W kw - 1 # 为了避免循环卷积产生混叠至少要补到 full 尺寸 size_h next_power_of_two(out_h_full) size_w next_power_of_two(out_w_full) img_fft np.fft.fft2(img, s(size_h, size_w)) kernel_fft np.fft.fft2(kernel_flipped, s(size_h, size_w)) conv_full np.real(np.fft.ifft2(img_fft * kernel_fft)) conv_full conv_full[:out_h_full, :out_w_full] if mode full: return conv_full if mode valid: # 去掉 full 结果中核覆盖不到的边缘区域 return conv_full[kh - 1 : out_h_full - (kh - 1), kw - 1 : out_w_full - (kw - 1)] if mode same: assert kh % 2 1 and kw % 2 1, same模式请使用奇数核 top (kh - 1) // 2 left (kw - 1) // 2 return conv_full[top:top H, left:left W] raise ValueError(f未知模式: {mode})注意第 2 行到第 4 行我把核翻转了。为什么因为滑动窗口法用的是“互相关”核不翻转而 FFT 的卷积定理天然计算的是标准数学卷积标准卷积本身要求核翻转。所以我这里先手动翻转一次两者利益就对齐了。如果你想用 FFT 实现严格数学卷积不翻转核那一步就可以去掉直接用kernel做fft2。3.3 几个容易让人懵的细节第一个细节是填充尺寸。为什么是H kh - 1因为卷积的 full 输出大小就是这个值。如果填充不够FFT 做完的是圆循环卷积而不是线性卷积结果会出现“周期性环绕”边缘那一圈会被旁边贴过来。这种问题很难一眼发现表现就是输出边缘出现不自然的条纹。第二个细节是我把尺寸又拉到了 2 的幂次。next_power_of_two这个辅助函数把完整尺寸扩大到比它大的最小 2 的幂。比如out_h_full513它就返回1024。NumPy 的 FFT 对任意尺寸都能算但 2 的幂次长度在底层实现里常常更快内存也不是问题。如果你做数组比较在意结果和滑窗法完全一致不需要非得是 2 的幂只要大于等于out_h_full即可。第三个细节是np.fft.fft2(img, s(size_h, size_w))里的s参数。它会在做 FFT 前自动把输入 pad 到指定尺寸省掉了手动np.pad。但要注意这种 pad 默认是在数组尾部补零而频域计算要求的是“周期延伸”两者只要尺寸足够大就等价放心使用。第四个细节是np.real。理论上图像和核都是实数组逆 FFT 的结果虚部应该接近零但由于浮点舍入误差会残留非常小的虚部。如果不取实部后续做图像显示或者和滑动窗口结果比较时类型和数值都会非常难受。3.4 FFT 什么时候会更快很多人以为 FFT 是万能加速实际用小核的时候它反而更慢。比如 3×3 高斯核滑动窗口法每输出一个像素只做 9 次乘法FFT 却要先对整张图做一次正变换、一次核变换、一次逆变换中间还会产生复数数组整体开销更大。但核一大效果就不一样了。比如同样在 1000×1000 图像上做滤波核从 3×3 变成 31×31滑动窗口法的乘法次数会从 900 万跳到 9.6 亿FFT 的复杂度主要取决于图像尺寸核变大基本不增加计算量。我实际测试下来核直径超过 15×15 以后FFT 版本通常更划算。到 31×31 时两者可能差一个数量级。如果你的场景是批量处理同一个核的多张图像FFT 优势还会更明显因为核的 FFT 只需要计算一次后面每张图只做一次正变换和一次逆变换。4. 两种方法的性能对比和选型经验4.1 复杂度对比表方法时间复杂度额外内存适合场景滑动窗口直接法O(H * W * kh * kw)极低视图不复制数据小核 3×3、5×5单张小图基于 kerneltensordot 的滑动窗口同上windows 是视图tensordot 会产生临时数组代码简洁推荐作为默认实现FFT 频域法O(HW log(HW))需要复数数组内存明显更高大核、大批量、超大图像表格里说的临时数组指的是sliding_window_view只有视图很省内存但np.tensordot在累加时会生成一个输出数组同时乘加过程中会有内部临时结果。这个开销远小于复制整堆窗口通常不用担心。4.2 我自己跑出来的体感数没有固定机器说数字意义不大但我可以给一个容易复现的脚本思路import time def bench(fn, img, kernel, modesame, repeat5): best float(inf) for _ in range(repeat): t0 time.perf_counter() fn(img, kernel, mode) best min(best, time.perf_counter() - t0) return best rng np.random.default_rng(42) img rng.random((512, 512)) small_kernel np.ones((3, 3), dtypenp.float64) / 9 large_kernel np.ones((31, 31), dtypenp.float64) / (31 * 31) print(bench(conv_sliding_mode, img, small_kernel)) print(bench(conv_fft, img, small_kernel)) print(bench(conv_sliding_mode, img, large_kernel)) print(bench(conv_fft, img, large_kernel))我这边跑下来的规律很稳定3×3核时滑动窗口法更快而且明显31×31核时 FFT 几乎要快一个数量级。在7×7到15×15之间两者可能接近具体看机器和图像尺寸。要注意FFT 在边界模式上也比滑窗版本更吃内存。假设输入和核都很大复数数组占用的空间大概是滑窗法的好几倍。如果你在内存紧张的设备上处理一张几千万像素图像优先选滑动窗口法。4.3 选型建议不复杂我自己的选型逻辑是核是 7×7 以下无脑用conv_sliding_mode核超过 15×15 或者要处理大批量图像切到conv_fft。中间范围不好说就在实际数据上跑一次基准测试。另外如果你主要做深度学习里的卷积层实现滑动窗口法更贴近后面的反向传播逻辑FFT 版的梯度推导会麻烦很多因为损失对频域乘法求导不太直观。还有一点如果输入是整数类型比如 uint8 图像建议在函数入口统一转成float64。整数类型在np.tensordot以及 FFT 中都可能带来溢出或截断尤其是图像像素值接近边界时边缘信息会直接丢。5. 实操中常见的坑和排查技巧5.1 两种方法结果对不上先查核到底翻没翻这是我在一些开源实现里见过最多的问题。滑动窗口法直接用核FFT 法标准卷积定理要求核翻转两边如果不统一结果就是错位的。尤其当核是非对称的比如 Sobel 水平边缘检测结果差一正负。排查方法很简单拿一个非对称核比如[[1, 2], [3, 4]]在随机小数组上分别用两个方法跑full模式再用np.allclose对比。对不上的时候手动把其中一个核翻转一下。如果是数学卷积直接保持一致即可。5.2 FFT 输出边缘出现鬼影多半是填充不够另一种典型症状输出中心区域和滑动窗口法很接近但边缘凹陷或者出现异常条纹。这基本是 FFT 尺寸没有填到H kh - 1。只要尺寸不足频域乘法会把本来不存在的“周期性邻居”混进来产生边缘混叠。解决办法就是那句老话先算out_h_full H kh - 1和out_w_full W kw - 1再拿去作为最小 FFT 尺寸。不要在输入上只补kh-1就完事。这里留一个检查技巧直接把 FFT 结果在full模式下和滑动窗口法对比看 max error 是否在1e-10级别。如果边缘差了几个像素值就是混叠。5.3 虚部、浮点数误差和数值类型每次做np.fft.ifft2都会得到 complex128 数组。直接那它去对比滑动窗口结果可能因为虚部而到处报错。正确做法是先np.real再比较。浮点数误差也非常正常。FFT 里面有多轮蝶形运算误差大概在1e-12到1e-14量级不要用比较用np.allclose并设置合理的atol比如1e-8。对于图像这种视觉数据完全看不出差别。另外如果输入是uint8在做乘法前一定要转float64。uint8相加会溢出变成 0 或 255FFT 里也会得到不可预测的复数结果错误会很隐蔽。5.4 用脉冲图像快速验证实现是否正确我调试卷积实现时最常用的一招是拿一张“只有中心点有值其他位置都是 0”的图去和一个全 1 核做卷积。作用相当于把核本身“复印”到输出上。def test_conv_valid(): img np.zeros((9, 9)) img[4, 4] 1.0 kernel np.ones((3, 3)) out_slide conv_sliding_mode(img, kernel, modefull) out_fft conv_fft(img, kernel, modefull) assert np.allclose(out_slide, out_fft, atol1e-8) print(full 模式一致性通过)因为脉冲图只有一个非零像素卷积结果本质就是核本身的位置偏移。这样不仅能验证 FFT 和滑动窗口的一致性还能直观看到 padding 和 offset 有没有写错。如果你把核换成类似[[1, 2], [3, 4]]的非对称核还能判断翻转方向是否正确。我自己的工作流里默认小核用滑动窗口法一旦核超过 15×15 或者输入超过 512×512我就切到 FFT。两种实现放在同一个工具模块里用同一个mode参数后面调试非常省心。你现在也可以拿自己的数组或图像去跑一下用np.allclose对比两个输出如果出现差异从核翻转和边界填充两处入手90% 的问题都能在这里找到。