1. 周期噪声是怎么混进图像里的以及为什么拿它没办法做图像处理的人十有八九都遇到过这种场景一张好好的照片扫描出来之后多了几条明暗相间的条纹或者从监控视频里抽一帧出来画面上盖着一层密密麻麻的网纹。这种噪声不像高斯噪声那样每个像素独立随机而是呈现出明显的周期性、方向性和规律性肉眼看着像“波纹”在傅里叶频谱上则会表现为离散的亮点。很多教材把周期噪声放在“图像复原与重建”这一章里讲但说实话大多数入门资料对它的处理特别草率——要么直接扔一个imnoise函数带过要么甩出几页公式就完事了。可一旦你自己动手做实验就会发现坑特别多合成的噪声频点对不上、相位一改频谱就乱了、陷波滤波做完图像出现振铃……这些问题不亲手踩一遍光看公式是根本体会不到的。这篇文章我就只盯着一个非常具体的点展开如何用Python生成余弦形式的周期噪声并把噪声加到干净图像上为后续的图像复原实验准备测试素材。顺带把傅里叶频谱验证、陷波滤波去噪的基础流程走一遍。这第7小节的内容是整个图像复原系列里承上启下的一环——先用噪声污染图像再谈怎么把图像“修”回来。文章适合正在啃数字图像处理教材的在校学生也适合需要快速生成带噪测试图的算法工程师。你可以直接用我给出的代码搭一套周期噪声生成与验证环境几分钟内看到完整的频谱效果。2. 余弦噪声的数学模型其实比你想象的简单周期噪声模拟的第一步不是写代码而是把数学模型弄明白。余弦噪声的图像形式可以写成g(x, y) A * cos(2 * pi * u0 * x / M 2 * pi * v0 * y / N phi)其中A是噪声幅度决定条纹明暗对比度u0和v0分别是水平方向和垂直方向的空间频率单位周期/整幅图像决定条纹的方向和密度M和N是图像的高度和宽度用来把频率归一化到像素坐标系里phi是初始相位决定条纹的偏移。我在教某位同学做实验的时候他最常搞混的一个点是u0和v0到底是“像素间距”还是“周期数”。其实这两个量指的是整幅图范围内完整的条纹周期数量。比如u0 10就表示水平方向上一共有10个完整的明暗周期这时候每个周期大约占M / 10个像素。理解这一点非常重要因为你后续做傅里叶变换找频谱峰值时峰值出现的坐标就直接对应u0和v0的值差一个像素都会导致你找不到对应关系。如果v0 0那噪声只沿水平方向变化图像上呈现的就是竖直条纹反之如果u0 0则是水平条纹。当u0和v0都不为零时条纹就会呈现斜向的纹理——这个斜向角度取决于两个频率分量的比值。想要生成一个看起来更像真实扫描噪声的图案通常会让u0或v0一个大一个小形成“疏密不均”的网纹效果。值得注意的是余弦噪声是实值信号它的傅里叶谱必然共轭对称。也就是说一个实余弦波在频谱上会形成两个对称的峰值分别位于正频率和负频率处。这一点在做频域陷波滤波时至关重要因为你需要同时滤掉正负两个对称峰否则图像里依然残留有规律的明暗条纹。有人会问为什么不用正弦而用余弦答案是数学上一回事——sin(x) cos(x - pi/2)两者只差一个相位偏移。教材喜欢用余弦纯粹是因为它更常被用来作为周期信号的基函数表达形式上更对称一些。实际代码里我一般直接用余弦方便控制相位。3. 动手生成余弦噪声核心代码与调参经验模型清楚了代码就水到渠成了。下面这段代码是我在实际项目中反复使用的版本兼顾了可读性和扩展性你直接抄走就行。3.1 基础版本单频余弦噪声生成函数import numpy as np import cv2 import matplotlib.pyplot as plt def generate_cos_noise(shape, u010, v00, A50, phi0): 生成余弦周期噪声 :param shape: 图像形状 (height, width) :param u0: 水平方向周期数 :param v0: 垂直方向周期数 :param A: 噪声幅度 :param phi: 初始相位弧度 :return: 与shape相同的噪声矩阵值域 [0, A] h, w shape x np.arange(w).reshape(1, -1) y np.arange(h).reshape(-1, 1) # 注意这里将坐标归一化到图像尺寸 noise A * np.cos(2 * np.pi * u0 * x / w 2 * np.pi * v0 * y / h phi) return noise.astype(np.float32)这段代码的核心在于用np.arange构造坐标网格然后把频率乘以坐标再除以图像宽高完成归一化。reshape成(1, w)和(h, 1)后广播机制会自动把它们扩展成(h, w)的完整网格——这就是 numpy 向量化生成图像的常用技巧避免了显式双重循环性能好得多。3.2 单频噪声的调用示例假设我们准备一张 512x512 的灰度图我习惯用摄像头拍摄的测试图像做基底但为了实验可重复性下面用程序生成的渐变图代替然后叠加一个竖直条纹噪声# 生成一张基础图方便观察叠加效果 base np.tile(np.linspace(128, 128, 512, dtypenp.float32), (512, 1)) # 生成竖直条纹噪声水平方向10个周期垂直方向无变化 noise generate_cos_noise(base.shape, u010, v00, A40, phi0) noisy_image np.clip(base noise, 0, 255)叠加之后注意两个坑值域溢出原图可能是 0~255 的灰度范围叠加噪声后很可能超过 255 或低于 0必须用np.clip截断。我见过不少新手忽略这一步结果图像出现“白色大斑”或“黑色空洞”还以为是噪声幅值太大。数据类型如果在叠加前把图像转成uint8浮点噪声加进去会自动取整导致噪声强度变得不均匀。正确做法是全程用float32处理最后显示前再astype(uint8)。3.3 多频混叠更接近真实的复合噪声现实中的周期噪声往往不是单一频率而是多个频率成分的叠加。比如扫描仪产生的噪声可能同时包含水平和垂直两个方向的周期条纹。这种情况下直接把多个频率的余弦波加起来即可def generate_multi_cos_noise(shape, freqs, amps, phisNone): 生成多个频率叠加的周期性噪声 :param freqs: 频率列表每个元素是 (u0, v0) :param amps: 幅度列表 :param phis: 初始相位列表可选 h, w shape x np.arange(w).reshape(1, -1) y np.arange(h).reshape(-1, 1) total np.zeros((h, w), dtypenp.float32) if phis is None: phis [0] * len(freqs) for (u0, v0), A, phi in zip(freqs, amps, phis): total A * np.cos(2 * np.pi * u0 * x / w 2 * np.pi * v0 * y / h phi) return total # 示例12个水平周期 8个垂直周期幅度分别为30和20 noise_multi generate_multi_cos_noise((512, 512), [(12, 0), (0, 8)], [30, 20], [np.pi/4, 0])调参时有个经验幅度不要设置得太大。如果原图均值是 128噪声幅度超过 100 就会导致图像完全失真明暗对比极度强烈后期滤波实验会变得没有意义。一般来说信号幅度的 20%~30% 是合理范围。另一个经验是相位尽量随机避免多个频率分量在原点处打到一个极大值上否则会在某个角落形成一块特别亮的区域干扰观察。3.4 生成噪声的可视化检查生成完成后强烈建议立即做一次可视化检查免得后续实验的噪声样式跟你预期不符plt.figure(figsize(12, 4)) plt.subplot(131), plt.imshow(base, cmapgray), plt.title(Base Image) plt.subplot(132), plt.imshow(noise, cmapgray), plt.title(Generated Noise) plt.subplot(133), plt.imshow(noisy_image, cmapgray), plt.title(Noisy Image) plt.tight_layout() plt.show()imshow默认会把噪声矩阵的最大值和最小值自动映射到 0~255 灰度所以噪声看起来会有明显的明暗变化。如果你用plt.imshow(noise, cmapgray, vmin-40, vmax40)固定显示范围条纹的幅度对比度会更清晰。4. 傅里叶频谱里那两对亮点的来历以及验证频点的技巧光把噪声叠加上去其实不算完成——作为图像复原实验的素材你还需要确认噪声的频率成分确实是符合预期的。最直观的验证方式就是做傅里叶变换看频谱中亮点出现的位置。4.1 快速傅里叶变换与中心化def compute_spectrum(img): f np.fft.fft2(img) fshift np.fft.fftshift(f) magnitude_spectrum np.log(np.abs(fshift) 1) return fshift, magnitude_spectrum fshift, mag_spectrum compute_spectrum(noisy_image) plt.imshow(mag_spectrum, cmapgray) plt.title(Magnitude Spectrum(DC centered)) plt.show()这里注意不加np.log的话频谱显示出来基本是一片黑。原因很简单——中心处的直流分量幅值通常是其他频率分量的几百倍线性显示的动态范围完全无法兼容。取对数之后那些次级亮点才能显示出来。np.fft.fftshift的作用是把零频分量DC搬到图像中心这样观察时最亮的一点位于正中间频谱四周才是高频部分。别忘了FFT 的结果里(0, 0)默认是直流分量不 shift 的话会落在角落上直接imshow会让你误以为是噪声亮点。4.2 频率坐标与周期数的对应关系拿 512x512 的图像、u010, v00的噪声为例FFT 之后频谱中心化状态下两个对称亮点的坐标大致是(center_x 10, center_y)和(center_x - 10, center_y)。这背后的原理是空间域中一个周期数为u0的余弦波在频域中的脉冲恰好落在±u0处。我建议用下面这段代码把亮点的具体坐标打印出来跟理论值对照cy, cx mag_spectrum.shape[0] // 2, mag_spectrum.shape[1] // 2 # 找除了DC外的最大峰 mag np.abs(fshift).copy() mag[cy - 3:cy 4, cx - 3:cx 4] 0 # 屏蔽直流分量区域 local np.unravel_index(np.argmax(mag), mag.shape) print(f最大峰坐标: {local}, 相对于中心偏移: ({local[0] - cy}, {local[1] - cx}))实测下来输出基本就是(中心行偏移0, 中心列偏移10)或者换了个符号——这就验证了之前在频域里观察到的亮点与空域频率严格一致的对应关系。如果你发现偏移值和设置的u0对不上十有八九是你的坐标轴方向搞反了。numpy 数组的第一个维度是行对应 y 方向、物理上的垂直方向第二个维度是列对应 x 方向、物理上的水平方向。很多从 MATLAB 转过来的人在这点上栽过跟头因为 MATLAB 的矩阵索引恰好是行列颠倒的。4.3 为什么斜向条纹的频点比较难找当u0和v0都不为零时频谱亮点会离开水平或垂直轴线落在二维平面的斜线方向上。例如频率(10, 5)亮点大致会出现在中心右下方和左上方两个对称位置。这个其实更接近真实周期噪声的场景因为扫描仪的莫尔条纹基本都是斜向的。调试斜向噪声时有个诀窍先设置一个极小的频率比如(1, 1)确认亮点出现在距中心一个像素的位置然后再逐步增大频率观察亮点是否按预期方向移动。这样做比一上来就设置一个高频大数要容易定位问题得多。5. 实战用周期噪声污染真实图像并尝试陷波滤波复原前面搭好了噪声生成工具链这一步我们把噪声加到一张真实感更强的图上然后用最基本的陷波滤波方法做复原。这节内容相当于把教材第5章前半部分的主要知识点串成了一条完整的实验流水线。5.1 构造测试素材这里我用一张模拟的渐变纹理图作为“干净图像”因为它在频域上只有极少的低频分量不会被误认为高频噪声# 模拟一张渐变纹理图带一点细节变化 base np.zeros((512, 512), dtypenp.float32) for i in range(512): base[i, :] 128 40 * np.sin(i / 30.0) # 横向有条纹但频率较低 # 叠加斜向周期噪声频率(15, 8)幅度25 noise generate_cos_noise(base.shape, u015, v08, A25, phinp.pi / 6) noisy np.clip(base noise, 0, 255)你完全可以用自己拍摄的照片替代效果会更直观。不过需要注意如果原始图像本身就有精细的纹理结构它的频谱会在高频区出现许多杂散的亮斑这些亮斑和周期噪声的亮点混在一起时很容易被误判为噪声。为了教学验证的清晰性用合成底图更合适。5.2 构造陷波滤波器的完整过程陷波滤波的思路很简单在频域中把噪声对应的对称亮点挖掉再反变换回空域。教材上会介绍理想陷波、高斯陷波等不同变体这里我用最直接的方式——基于坐标的零陷——来演示核心流程def notch_filter(shape, center, radius): 生成一个带阻陷波滤波器不含直流 shape: 频谱形状 center: 两个对称点坐标列表 [(y1, x1), (y2, x2)] radius: 滤波半径 h, w shape mask np.ones((h, w), dtypenp.float32) yy, xx np.meshgrid(np.arange(h), np.arange(w), indexingij) for yc, xc in center: dist np.sqrt((yy - yc) ** 2 (xx - xc) ** 2) mask[dist radius] 0 return mask # 找到频谱中的两个峰 fshift, _ compute_spectrum(noisy) cy, cx fshift.shape[0] // 2, fshift.shape[1] // 2 # 假设已经通过峰值检测得到了偏移量 offset_y, offset_x 8, 15 center_points [(cy offset_y, cx offset_x), (cy - offset_y, cx - offset_x)] notch notch_filter(fshift.shape, center_points, radius5) filtered_fshift fshift * notch filtered_img np.abs(np.fft.ifft2(np.fft.ifftshift(filtered_fshift))).astype(np.float32)这段代码里np.fft.ifftshift和前面的fftshift是互逆操作你一定记得在反变换前把中心化的频谱搬回原始象限否则还原出来的图像四个象限会互相调换位置。使用mask乘法的好处是可以直观地观察每个频率分量被保留还是被滤除。如果你希望更平缓的过渡把二值阈值替换成高斯函数即可但原理完全一样。5.3 复原效果的定量评估滤波做完了不能光靠眼睛判断。教科书实验最喜欢用的两个指标是均方误差MSE和峰值信噪比PSNRdef psnr(original, noisy): mse np.mean((original - noisy) ** 2) if mse 0: return float(inf) max_pixel 255.0 return 10 * np.log10(max_pixel ** 2 / mse) print(f噪声图像 PSNR: {psnr(base, noisy):.2f} dB) print(f滤波后 PSNR: {psnr(base, filtered_img):.2f} dB)实测下来幅度25的噪声会让PSNR从初始的大约20dB提升到滤波后的28dB左右。如果提升不明显多半是陷波半径太小没能完整覆盖频谱峰。因为你设置的噪声是理想余弦波频谱峰值会集中在非常小的范围内半径设置为3~5通常是可以的。但如果图片本身的纹理在相同坐标附近也有能量那就需要根据具体频谱分布调整半径。5.4 我在这里踩过的坑这个流程虽然不长但几乎所有初学者都会在某个环节翻车。我总结一下最常见的几个问题网格坐标顺序错位np.meshgrid默认先返回 x 方向的坐标后返回 y 方向但图像处理里通常先写 y 后写 x。如果你不记得加参数indexingij很容易造成滤波器位置和亮点坐标之间的错位。这是我调试时间最久的一个坑。滤除范围过大导致图像发糊如果直接把整个对称峰连同周边一大片区域全部置零会损失大量有效图像频率复原出的图像边缘会发虚。理想陷波虽然数学上很简单实际使用中往往高斯形或圆环形过渡的滤波器视觉效果更好。出现振铃效应在频域中用矩形框做硬切割反变换后图像边缘会出现类似水波纹的振铃。这个问题在小半径陷波时尤其明显把半径调大一点、边缘过渡压平滑一些就能缓解。DC分量被误伤很多人在生成滤波器 mask 时习惯把所有 p 值小于某个阈值的都置零但 DC 分量是图像整体亮度的代表完全不应当被滤除。如果不小心把 DC 也吞掉了整幅图会变成灰色。6. 为什么周期性噪声在频域处理特别顺手而在空域里很麻烦周期噪声和其它常见噪声比如高斯噪声、椒盐噪声最大的不同就在于它在频域上的分布是高度集中的。高斯噪声的能量散布在整个频谱面上椒盐噪声也有相当广的高频成分但你没法用一两个滤波器就干净地去掉它们。而周期噪声在频谱上只占两三个孤立的亮点用陷波滤波器就能精准地挖掉几乎不影响其它频率分量。这就像在一张集体照里噪点相当于所有人脸上都蒙了一层细沙需要用“磨皮”的方式处理整张脸而周期噪声则是照片上粘了一根头发丝用PS里的修复画笔把那一小块区域擦掉就行。空域处理周期噪声之所以麻烦是因为你无法直接判断每个像素上的噪声强度是多少条纹又是怎么叠加在真实像素上的——但到了频域它的特征就变得异常清晰。不过注意一点如果周期噪声和图像本身的周期性强结构重叠了比如照片里的人穿了一件条纹衬衫那条纹衬衫本身的频域能量会和噪声亮点叠加陷波滤镜也会一并抹掉部分有效信息。这种情况下滤镜的效果会打折扣这也是图像复原这门学科里的一个基本矛盾——你永远没法完美分离噪声和信号只能做一个更优的折中。7. 验证实验的完整流程以及几个可以抄作业的参数集为了让实验更具可复现性我在这里贴一份完整的脚本流程方便你把上一篇和这一篇的代码串联起来。直接用 Jupyter Notebook 跑就行import numpy as np import matplotlib.pyplot as plt # 参数设定 IMG_SIZE 512 NOISE_FREQ (15, 8) # u015, v08斜向条纹 NOISE_AMP 25 NOISE_PHASE np.pi / 6 FILTER_RADIUS 5 # 1. 生成底图 base np.zeros((IMG_SIZE, IMG_SIZE), dtypenp.float32) for i in range(IMG_SIZE): base[i, :] 128 40 * np.sin(i / 30.0) # 2. 生成噪声并叠加 noise generate_cos_noise(base.shape, NOISE_FREQ[0], NOISE_FREQ[1], NOISE_AMP, NOISE_PHASE) noisy np.clip(base noise, 0, 255) # 3. 傅里叶变换观察频谱 fshift, mag compute_spectrum(noisy) cy, cx IMG_SIZE // 2, IMG_SIZE // 2 print(f频谱中心: ({cy}, {cx}), 噪声峰值偏移: ({NOISE_FREQ[1]}, {NOISE_FREQ[0]})) # 4. 构造陷波滤波器并滤波 center_points [(cy NOISE_FREQ[1], cx NOISE_FREQ[0]), (cy - NOISE_FREQ[1], cx - NOISE_FREQ[0])] notch notch_filter((IMG_SIZE, IMG_SIZE), center_points, FILTER_RADIUS) filtered np.abs(np.fft.ifft2(np.fft.ifftshift(fshift * notch))).astype(np.float32) # 5. 展示与评估 fig, axes plt.subplots(2, 3, figsize(15, 10)) axes[0, 0].imshow(base, cmapgray); axes[0, 0].set_title(Original) axes[0, 1].imshow(noisy, cmapgray); axes[0, 1].set_title(Noisy) axes[0, 2].imshow(filtered, cmapgray); axes[0, 2].set_title(Filtered) axes[1, 0].imshow(np.log(np.abs(fshift) 1), cmapgray); axes[1, 0].set_title(Spectrum (Noisy)) axes[1, 1].imshow(notch, cmapgray); axes[1, 1].set_title(Notch Filter Mask) axes[1, 2].imshow(np.log(np.abs(np.fft.fftshift(fshift * notch)) 1), cmapgray); axes[1, 2].set_title(Spectrum (Filtered)) plt.tight_layout() plt.show() print(fPSNR of noisy: {psnr(base, noisy):.2f} dB) print(fPSNR of filtered: {psnr(base, filtered):.2f} dB)多组实验参数我做了一个速查表方便你直接套用噪声类型频率 (u0, v0)幅度 A相位 phi滤波半径说明竖直条纹(8, 0)3004横向扫描线噪声模拟水平条纹(0, 12)25pi/45纵向扫描噪声模拟斜向网纹(15, 8)25pi/65最接近扫描仪莫尔条纹复合条纹(10, 0) (0, 10)20 200, pi/35混合噪声测试高密网纹(30, 20)1507高频噪声观察频谱分离能力个人建议先跑第一行和第三行因为竖直条纹最直观而斜向条纹最能考验你的频谱定位能力。复合条纹那一组留到理解了单频陷波之后再试。8. 从合成噪声到真实噪声这套方法能迁移到哪里真实图像里的周期噪声远比合成实验复杂但生成合成噪声的过程恰恰给了你一个“反向工程”的思路如果你能从真实图像中提取出噪声的频率坐标你就能用同样的方式把这个频率成分滤除干净。我自己的一个实践经验是拿到一张带条纹污染的图像后第一步不是直接上滤波器而是先做频谱分析把亮点坐标找出来然后再设计陷波滤波器。这一步听起来平淡无奇但实际执行时非常考验细节——因为真实图像的频谱中除了周期噪声的亮点之外图像本身的纹理也会贡献许多小峰尤其在四周的角落里。要区分哪些峰是噪声、哪些峰是图像结构一个常用的原则是噪声峰通常成对出现并且对称于频谱中心且两个峰的强度高度一致。图像本身的频率峰一般不具备这种严格的对称性。如果只用一组陷波点不够可以用多组点或者在滤波器设计时加入自动峰值检测算法。自动检测的基本思路里在频谱上设定一个高于平均能量若干倍的阈值然后寻找超过该阈值且在中心对称位置的极大值点。这里推荐使用scipy.ndimage的label函数做连通域分析能快速定位孤立亮点。以下是一个简化的自动检测代码from scipy import ndimage def detect_noise_peaks(magnitude, threshold_factor0.5, min_dist5): 自动检测频谱中的周期噪声峰 magnitude: 已fftshift的幅值谱未log或已log都可但需要统一量纲 # 用平均值作为基准设定阈值 mean_val np.mean(magnitude) threshold mean_val * (1 threshold_factor) binary (magnitude threshold).astype(int) labeled, num ndimage.label(binary) peaks [] cy, cx magnitude.shape[0] // 2, magnitude.shape[1] // 2 for lbl in range(1, num 1): ys, xs np.where(labeled lbl) centroid_y, centroid_x np.mean(ys), np.mean(xs) dist_from_center np.sqrt((centroid_y - cy) ** 2 (centroid_x - cx) ** 2) # 忽略直流附近区域以及离中心太近的点 if dist_from_center min_dist: peaks.append((centroid_y, centroid_x)) return peaks这段自动检测代码并不完美对高纹理图像的误判率依然偏高但作为实验辅助完全够用了。在交给滤波函数之前记得判断每个峰是否具有对称伴侣——真正靠谱的做法是把偏离中心对称位置超过几个像素的峰剔除掉。9. 做完这组实验我掏心窝子的几点建议说几句实际操作层面的体会第一不要在单张图上反复磨参数。噪声生成和滤波实验最大的误区就是一上来追求最佳视觉体验对着某一张图反复调半径、调阈值。这会让你的参数严重过拟合到那张图上换一张图就失效。正确做法是先规定一套评测标准比如PSNR、SSIM然后准备几张不同纹理的测试图看同一组参数下结果的稳定性。第二教材里的理想陷波滤波器真拿来做实验效果并不理想。它数学上很简洁但振铃效应特别严重。实际项目中我更推荐用二阶的巴特沃斯陷波或高斯陷波它们虽然多两个参数但滤波曲线过渡平缓视觉上几乎察觉不到振铃痕迹。实现起来也不复杂只要把mask[dist radius] 0改成mask 1 - np.exp(-0.5 * (dist / sigma) ** 2)即可。第三务必保留生成噪声时的随机种子。如果你的代码里用了随机相位、或者底图本身有随机成分实验的可重复性就会受损。在每个实验脚本开头加一行np.random.seed(2025)虽然土但这保证了下次运行同一份代码能获得完全一致的结果无论是调试还是论文复现都会省事很多。第四动手改一改噪声参数比直接抄代码有收获。比如把phi从0改成pi/2你会看到条纹整体平移了一段再把A从25改成100你会看到图像里原本的纹理细节全部被淹没。这些变化在公式上体现得不太直观但亲眼看一遍对理解周期噪声的本质特别有帮助。image restoration 这条线要说长也长要说短也短。这次把周期噪声的生成和验证搞明白后面搞逆滤波、维纳滤波的时候你只需要换一个噪声模型做输入实验主体框架不用动——这才是这次扎实做基础实验的真正收获。