1. 项目整体思路拆解一次搞定FFT去噪和数值积分先把这个标题拆开看里面其实是两个完全独立的经典小项目一个是图像去噪利用傅里叶变换把噪声在频域里“抠掉”另一个是求圆周率利用定积分把半圆面积算出来。它们唯一的共同点是都依赖Scipy这个库一个用fft模块一个用integrate模块。这种组合式Demo特别适合正在学Python科学计算的人因为两个任务都足够短小但背后牵扯的原理一个比一个深。先说图像去噪。传统空间域处理比如中值滤波、均值滤波都是直接滑动窗口去改像素值。傅里叶变换的思路完全不一样它先把整张图从“像素空间”变到“频率空间”噪声在频率空间里通常表现为密密麻麻的高频分量而图像主体往往是低频分量。我们在频域里把高频部分压掉再逆变换回空间域噪声就没了。这个思路一旦打开你会发现它不只是去噪还能用来做边缘提取、图像压缩甚至水印去除。再说定积分求圆周率。圆周率π是圆的周长和直径之比看起来和积分八竿子打不着。但如果你把单位圆的上半部分看成函数y sqrt(1 - x²)从-1到1对这个函数做定积分算出来的正好是半个圆的面积也就是π/2。那么π就等于这个积分结果乘以2。这个方案最妙的地方在于它把几何问题转化成了数值计算问题让我们可以顺理成章地引入scipy.integrate.quad这个自适应积分器。两个任务放在一起还有个额外好处FFT本身是一种精确变换误差可控而数值积分天生就是近似计算误差可以量化。这一精一近似形成对比能让初学者更直观地理解“数值计算到底在干什么”。我建议你按顺序先做去噪再做积分因为去噪的代码更短视觉反馈也更直观做完会特别有成就感。1.1 核心需求解析标题背后到底想考你什么这个标题看着像两个小练习拼在一起其实考察点非常明确。第一部分考的是你是否理解傅里叶变换的物理意义而不是只会调用fft2然后瞎调参数。第二部分考的是你对定积分基本概念和自适应性数值算法的掌握程度。如果你只是把代码跑通那你最多只完成了30%剩下70%在于你能不能解释清楚为什么噪声在高频为什么quad能算这种带端点的开根号函数为什么返回值里还带着一个误差估计从实操角度看还有一个隐含考察点环境搭建。Scipy不是Python自带库需要安装numpy、scipy、matplotlib三个基础包。很多人卡在去噪显示图像这一步其实只是没装matplotlib。我建议直接用pip install numpy scipy matplotlib一把梭然后用Jupyter Notebook来做因为每一步都可以即时看到频谱图和滤波结果比写一个完整脚本调试要舒服太多。1.2 为什么是Scipy选型背后的真实原因有人可能会问FFT用numpy.fft不就行了吗为什么非要Scipy答案要从两方面看。numpy.fft确实提供了基础的fft2和ifft2够用但接口比较糙而scipy.fft设计得更现代它对FFT的封装支持worker线程、计划缓存并且和老代码的兼容性更好尤其在处理大规模数组时速度优势明显。scipy.fft还有一个加分项它包含rfft2、irfft2这种专门处理实数输入的快速路径图像恰恰是实数的用实数FFT能省下一半内存和可观的计算时间。积分这边就更不用说了。scipy.integrate.quad包装的是经典的数值积分库内部采用自适应Gauss-Kronrod方法。你不用关心它内部怎么选点它自己会在函数变化剧烈的地方加密采样在平缓的地方减少采样最终返回一个高精度的近似值和误差估计。你要是自己写个矩形法或者辛普森法想要达到quad的精度计算量可能要差出好几个数量级。我给一个非常直观的建议先跟着本文用scipy.fft和scipy.integrate.quad跑通流程之后再去试试numpy.fft和手写积分对比一下性能和代码复杂度你会对“为什么科学计算界普遍使用成熟库”有切身体会。2. 傅里叶变换去噪从数学到代码的完整落地2.1 图像傅里叶变换的核心直觉想搞懂图像FFT去噪先建立直觉别急着看公式。把图像想象成一张起伏不平的地形图像素亮度就是海拔。傅里叶变换做的事情就是把这整片地形拆解成无数不同频率、不同方向的正弦波叠加。低频正弦波是那种变化缓慢的大尺度起伏对应图像里渐变区域和主体轮廓高频正弦波是那种变化剧烈的细小毛刺对应图像边缘、纹理还有噪声。噪声为什么基本都在高频因为噪声的特点是随机、突兀、局部突变比如椒盐噪声里那些黑白点和周围像素亮度差异极大。从波的角度看这种突变只能靠高频率的正弦波才能拟合出来。所以去噪的基本策略也就清楚了在频域里把高频成分衰减掉保留低频主体再逆变换回去。这里有一个很多人忽略的点图像边缘本身也是高频。低通滤波去噪的同时必然会让边缘变模糊。这是个物理上无法避免的权衡你能做的只是控制衰减强度在去噪和保留细节之间找一个平衡点。实操中我建议你每次改变滤波半径后都输出一张滤波后图像观察边缘变化而不是只看一个PSNR数值。2.2 生成一张带噪声的测试图先造一张干净的测试图再加噪声这样才能对比去噪前后的效果。我一般用圆盘图原因是它有明显的边缘和平坦区域去噪效果肉眼一看便知。下面这段代码生成一张256×256的圆盘图然后加上椒盐噪声。import numpy as np import matplotlib.pyplot as plt from scipy import fft # 生成一个圆盘图像 x np.linspace(-2, 2, 256) X, Y np.meshgrid(x, x) img 255 * ((X ** 2 Y ** 2) 1).astype(float) # 添加椒盐噪声 rng np.random.default_rng(42) noise rng.random(img.shape) img_noisy img.copy() img_noisy[noise 0.01] 0 # 随机黑点 img_noisy[noise 0.99] 255 # 随机白点 plt.figure(figsize(10, 5)) plt.subplot(1, 2, 1) plt.imshow(img, cmapgray) plt.title(Original) plt.axis(off) plt.subplot(1, 2, 2) plt.imshow(img_noisy, cmapgray) plt.title(Noisy) plt.axis(off) plt.show()我没有用高斯噪声而选椒盐噪声是因为椒盐噪声在频域里的特征更典型稀疏但强烈的冲激频谱上呈现为散布在整个高频区域的亮点。高斯噪声则是均匀地铺在所有频率上处理起来更麻烦。对教学来说椒盐噪声的滤波效果更明显能让你一眼看到魔法发生。2.3 频域滤波核心代码每个参数都有讲究傅里叶变换去噪的标准流程五步二维FFT、频谱中心化、构造掩膜、频域滤波、逆变换。每一步都有讲究我拆开说。# 1. 二维FFT spectrum fft.fft2(img_noisy) # 2. 频谱中心化把低频移到图像中心 spectrum_shifted fft.fftshift(spectrum) # 3. 构造圆形低通掩膜 rows, cols img.shape crow, ccol rows // 2, cols // 2 radius 40 yy, xx np.ogrid[:rows, :cols] mask ((xx - ccol) ** 2 (yy - crow) ** 2) radius ** 2 # 4. 频域滤波 filtered_spectrum spectrum_shifted * mask # 5. 逆变换取实部 filtered_img fft.ifft2(fft.ifftshift(filtered_spectrum)).real plt.figure(figsize(12, 4)) plt.subplot(1, 3, 1) plt.imshow(np.log(1 np.abs(spectrum_shifted)), cmapgray) plt.title(Spectrum) plt.axis(off) plt.subplot(1, 3, 2) plt.imshow(mask, cmapgray) plt.title(Mask (radius40)) plt.axis(off) plt.subplot(1, 3, 3) plt.imshow(filtered_img, cmapgray) plt.title(Filtered) plt.axis(off) plt.show()为什么要fftshift因为fft2算出来的频谱低频在四个角落高频在中心附近这不符合我们“中心放低频”的视觉习惯。fftshift把频谱四角的低频挪到中心方便画圆掩膜。做完掩膜乘法后必须用ifftshift把频谱还原回原来的排列方式再做逆变换这一步顺序反了会得到一张莫名其妙的花屏。为什么要圆形掩膜而不是方形因为方形掩膜会在滤波后产生明显的振铃效应图像边缘出现一圈圈水波纹。圆形掩膜对应的频域边界是连续旋转对称的振铃会小得多。这是图像处理里的老经验新手经常在这里翻车。radius40这个值不是拍脑袋定的它对应的是保留频率范围。图像尺寸256中心化后最大频率半径约为128。圆盘图的主体是低频半径40已经能覆盖绝大多数主体频率了。你可以试试调成20看效果圆盘边缘会变糊但噪声几乎全没了再调成80边缘清晰了但噪声点会残留不少。这就是上一节说的权衡。2.4 去噪效果评估别只顾着看漂不漂亮最关键的一步是评估效果。我常用的方法是同时输出三张图原图、去噪图、噪声残差图。残差就是原图减去去噪图理想情况下残差里只能看到噪声和一些边缘细节如果残差里出现明显的圆形轮廓说明低通滤波把边缘信息也抹掉了。定量评估可以用信噪比SNR公式是信号能量除以噪声能量。不过对于教学项目肉眼判断反而更直接。我最喜欢的一个做法是放大局部区域对比比如用plt.xlim/ylim截取圆盘边界附近放大能很清楚看到噪声点是否被抹掉、边缘是否变虚。还有一个小技巧低通滤波可以再做一次高通实验来加深理解。把掩膜取反只保留高频再逆变换你会发现得到的图几乎是纯边缘和噪声的叠加。这就能直观印证“细节和噪声都在高频”这个结论。我强烈建议你跑一下这个对比实验比你单独做低通滤波理解要深得多。3. 定积分求圆周率quad函数背后的数值艺术3.1 数学原理半圆面积就是π/2圆周率与定积分的关系本质上是圆的面积公式。单位圆的方程是x² y² 1取上半部分得到y sqrt(1 - x²)。对这个函数在[-1, 1]上做定积分得到的是上半圆的面积也就是π/2。所以π 2 × ∫_{-1}^{1} sqrt(1 - x²) dx。如果你觉得从-1到1麻烦也可以只算第一象限的部分x从0到1积分得到四分之一圆面积π/4最后乘4。两种写法数值上没差别但四分之一圆的写法积分区间更短对数值算法来说通常更友好。这个题目看起来有点“绕”但它的教学价值极高它让你意识到很多看起来和微积分无关的常数都可以通过数值积分算出来。你还可以扩展出很多变体比如用同样的方法算椭圆面积、球体体积只要写出被积函数quad就能处理。3.2 quad函数实操一行代码算出半圆面积直接上代码这是整个项目里最短的一段但信息量一点不少。from scipy.integrate import quad import numpy as np def half_circle(x): return np.sqrt(1 - x ** 2) area_half, error_estimate quad(half_circle, -1, 1) pi_approx 2 * area_half print(f半圆面积 {area_half:.16f}) print(f误差估计 {error_estimate:.2e}) import math print(fπ 近似值 {pi_approx:.16f}) print(fπ 真实值 {math.pi:.16f}) print(f绝对误差 {abs(pi_approx - math.pi):.2e})我跑这段代码得到的结果是半圆面积输出1.5707963267948966π近似值3.141592653589793误差在10的负16次方量级。这个精度已经接近双精度浮点数的极限了相当惊人。重点看quad的第二个返回值error_estimate这是算法自己估算的绝对误差。数值积分最妙的地方在于它不但给你结果还附带一个可信度指标。实际使用中我会习惯性把error_estimate打出来看一眼如果这个误差估计很大说明被积函数可能有问题比如奇异点、震荡剧烈或者区间选择不当。3.3 为什么quad这么准自适应算法揭秘也许你会想不就是求积分吗我自己把区间切成10000份用矩形法累加也能算出差不多的π吧。没错但效率差太多了。我做过一个简单测试用矩形法计算同样的积分要达到误差10的负8次方需要大约50万个月。而quad内部采用自适应Gauss-Kronrod积分它在函数变化缓的地方少采样变化剧的地方多采样。举例来说sqrt(1 - x²)在x接近±1时函数曲线几乎垂直导数趋向无穷。普通均匀采样在端点附近会很吃亏因为这一段函数变化太快均匀采样的密度根本不够。quad的聪明之处就在这里它会先在整体上采样一组点然后检查哪些区域的误差贡献大再对误差大的区域递归加密采样直到整个区间都满足精度要求。这种“把力气花在刀刃上”的策略正是自适应算法的精髓。理解了这一点你就能解释为什么同样的积分quad耗时不到一毫秒精度却碾压手写矩形法几个数量级。3.4 扩展用quad做更复杂的圆面积计算既然标题是“定积分求圆周率”我建议你做完基本版后顺手做两个扩展练习。第一个扩展是算任意半径R的圆面积。只需要在half_circle里加入半径参数用quad(half_circle, -R, R, args(R,))传入。这里args的用法很实用后续做参数扫描时会经常用到。def circle_upper(x, R): return np.sqrt(R ** 2 - x ** 2) for R in [1, 2, 0.5]: area, err quad(circle_upper, -R, R, args(R,)) print(f半径{R}的圆面积: {2 * area:.10f}, 误差估计: {err:.1e})第二个扩展是修改被积函数算半个椭圆的面积。椭圆方程x²/a² y²/b² 1上半部分y b * sqrt(1 - x²/a²)积分区间[-a, a]。结果应该是πab除以2。这个扩展能让你真正理解“定积分求面积”的本质边界函数求积分得到的是曲线下的面积如果把边界函数改一下就能算各种稀奇古怪的图形面积。4. 实战中的坑与排查技巧4.1 FFT谱图里那些奇怪的亮点和条纹第一次把频谱用imshow显示出来的时候很多人会吓一跳为什么中心一个超亮的大斑点四周还有横竖十字形的亮线中心大斑点其实是直流分量也就是全图亮度的平均值。如果是一个全黑背景的圆盘图直流分量会很高这是正常的。横竖亮线通常来自图像边缘的不连续性图像左右边界、上下边界在傅里叶眼里都是突变所以会在对应的频率轴上产生一条亮线。处理这类问题一个经典做法是给图像做边缘平滑或者加窗。不过在我们这个去噪项目里并不需要特别处理。只要你清楚这些现象的来源就不会被它们带偏思路。4.2 逆变换结果出现复数怎么办调用fft.ifft2之后返回结果其实是复数数组。这是因为浮点运算误差理论上实部才是真实图像虚部会有极小的残留。很多人忘了取实部直接plt.imshow复数数组结果要么报错要么显示出一片乱码。正确做法是调用.real取实部代码里我已经写了千万别省。还有一个常见坑是把ifftshift放在fftshift前面。正确的顺序是fft2-fftshift- 滤波 -ifftshift-ifft2。如果你把ifftshift和fftshift搞混唯一的表现就是输出图像错乱这时候先检查这一步。4.3 quad函数报错和异常情况处理quad最常见的报错是ValueError: Inf or NaN这通常是被积函数在积分区间内出现无穷大。比如你把被积函数写成1 / sqrt(1 - x**2)在x ±1处分母为零函数值无穷大但恰恰这个积分是收敛的值等于π。quad其实能处理一部分端点奇异的情况但需要你显式设置points参数把奇异点指出来。# 端点奇异但积分收敛的例子 integral, err quad(lambda x: 1 / np.sqrt(1 - x**2), -1, 1, points[-1, 1]) print(integral) # 接近 pi设置了points之后quad会在这些点附近加密采样从而正确处理奇异行为。这个知识点在很多工程计算里都会用到值得记一下。另一个我遇到过的问题是积分区间过大导致被积函数在中间出现剧烈震荡。这时候quad可能报“Maximum number of subdivisions exceeded”的警告解决办法是把大区间拆成多个小区间分别积分再把结果相加。实际工程里遇到这类问题永远优先考虑分区间而不是盲目调limit参数。4.4 常见问题速查表症状可能原因解决办法频谱中心亮点巨大图像直流分量过高正常现象可做对数变换再显示滤波后图像有波纹使用了方形掩膜换成圆形掩膜降低振铃逆变换输出乱码ifftshift位置错误或没取实部检查顺序加.real去噪后图像过糊滤波半径太小增大radius重新尝试quad报Inf或NaN被积函数含奇异点用points参数标示奇异点quad耗时长且警告积分区间过大或震荡剧烈手动拆分积分区间分别计算这张表是我把实际踩过的坑整理出来的。做图像处理时遇到视觉异常第一反应先看频谱和掩膜形状对不对做数值积分时遇到报错第一反应先把被积函数图像画出来。把函数可视化很多问题一眼就能定位。我个人在这次实践中最大的收获是真正建立了“频域思维”以前遇到图像噪声第一反应是换滤波器窗口大小现在会先想到去频谱里看一眼噪声长什么样。这种思维方式的转变远比记住几个函数调用有价值。至于定积分求π看起来像个玩具但当你意识到quad能处理带奇异点的函数时你就会发现工程里很多复杂定积分都可以放心交给它前提是你要理解算法适应的边界在哪里。希望这两个小项目能帮你打开科学计算的大门。 
长河君 · 2026/10/16 8:26:45 · 阅读 1.2k Scipy实战：傅里叶变换去噪与定积分求圆周率
1. 项目整体思路拆解一次搞定FFT去噪和数值积分先把这个标题拆开看里面其实是两个完全独立的经典小项目一个是图像去噪利用傅里叶变换把噪声在频域里“抠掉”另一个是求圆周率利用定积分把半圆面积算出来。它们唯一的共同点是都依赖Scipy这个库一个用fft模块一个用integrate模块。这种组合式Demo特别适合正在学Python科学计算的人因为两个任务都足够短小但背后牵扯的原理一个比一个深。先说图像去噪。传统空间域处理比如中值滤波、均值滤波都是直接滑动窗口去改像素值。傅里叶变换的思路完全不一样它先把整张图从“像素空间”变到“频率空间”噪声在频率空间里通常表现为密密麻麻的高频分量而图像主体往往是低频分量。我们在频域里把高频部分压掉再逆变换回空间域噪声就没了。这个思路一旦打开你会发现它不只是去噪还能用来做边缘提取、图像压缩甚至水印去除。再说定积分求圆周率。圆周率π是圆的周长和直径之比看起来和积分八竿子打不着。但如果你把单位圆的上半部分看成函数y sqrt(1 - x²)从-1到1对这个函数做定积分算出来的正好是半个圆的面积也就是π/2。那么π就等于这个积分结果乘以2。这个方案最妙的地方在于它把几何问题转化成了数值计算问题让我们可以顺理成章地引入scipy.integrate.quad这个自适应积分器。两个任务放在一起还有个额外好处FFT本身是一种精确变换误差可控而数值积分天生就是近似计算误差可以量化。这一精一近似形成对比能让初学者更直观地理解“数值计算到底在干什么”。我建议你按顺序先做去噪再做积分因为去噪的代码更短视觉反馈也更直观做完会特别有成就感。1.1 核心需求解析标题背后到底想考你什么这个标题看着像两个小练习拼在一起其实考察点非常明确。第一部分考的是你是否理解傅里叶变换的物理意义而不是只会调用fft2然后瞎调参数。第二部分考的是你对定积分基本概念和自适应性数值算法的掌握程度。如果你只是把代码跑通那你最多只完成了30%剩下70%在于你能不能解释清楚为什么噪声在高频为什么quad能算这种带端点的开根号函数为什么返回值里还带着一个误差估计从实操角度看还有一个隐含考察点环境搭建。Scipy不是Python自带库需要安装numpy、scipy、matplotlib三个基础包。很多人卡在去噪显示图像这一步其实只是没装matplotlib。我建议直接用pip install numpy scipy matplotlib一把梭然后用Jupyter Notebook来做因为每一步都可以即时看到频谱图和滤波结果比写一个完整脚本调试要舒服太多。1.2 为什么是Scipy选型背后的真实原因有人可能会问FFT用numpy.fft不就行了吗为什么非要Scipy答案要从两方面看。numpy.fft确实提供了基础的fft2和ifft2够用但接口比较糙而scipy.fft设计得更现代它对FFT的封装支持worker线程、计划缓存并且和老代码的兼容性更好尤其在处理大规模数组时速度优势明显。scipy.fft还有一个加分项它包含rfft2、irfft2这种专门处理实数输入的快速路径图像恰恰是实数的用实数FFT能省下一半内存和可观的计算时间。积分这边就更不用说了。scipy.integrate.quad包装的是经典的数值积分库内部采用自适应Gauss-Kronrod方法。你不用关心它内部怎么选点它自己会在函数变化剧烈的地方加密采样在平缓的地方减少采样最终返回一个高精度的近似值和误差估计。你要是自己写个矩形法或者辛普森法想要达到quad的精度计算量可能要差出好几个数量级。我给一个非常直观的建议先跟着本文用scipy.fft和scipy.integrate.quad跑通流程之后再去试试numpy.fft和手写积分对比一下性能和代码复杂度你会对“为什么科学计算界普遍使用成熟库”有切身体会。2. 傅里叶变换去噪从数学到代码的完整落地2.1 图像傅里叶变换的核心直觉想搞懂图像FFT去噪先建立直觉别急着看公式。把图像想象成一张起伏不平的地形图像素亮度就是海拔。傅里叶变换做的事情就是把这整片地形拆解成无数不同频率、不同方向的正弦波叠加。低频正弦波是那种变化缓慢的大尺度起伏对应图像里渐变区域和主体轮廓高频正弦波是那种变化剧烈的细小毛刺对应图像边缘、纹理还有噪声。噪声为什么基本都在高频因为噪声的特点是随机、突兀、局部突变比如椒盐噪声里那些黑白点和周围像素亮度差异极大。从波的角度看这种突变只能靠高频率的正弦波才能拟合出来。所以去噪的基本策略也就清楚了在频域里把高频成分衰减掉保留低频主体再逆变换回去。这里有一个很多人忽略的点图像边缘本身也是高频。低通滤波去噪的同时必然会让边缘变模糊。这是个物理上无法避免的权衡你能做的只是控制衰减强度在去噪和保留细节之间找一个平衡点。实操中我建议你每次改变滤波半径后都输出一张滤波后图像观察边缘变化而不是只看一个PSNR数值。2.2 生成一张带噪声的测试图先造一张干净的测试图再加噪声这样才能对比去噪前后的效果。我一般用圆盘图原因是它有明显的边缘和平坦区域去噪效果肉眼一看便知。下面这段代码生成一张256×256的圆盘图然后加上椒盐噪声。import numpy as np import matplotlib.pyplot as plt from scipy import fft # 生成一个圆盘图像 x np.linspace(-2, 2, 256) X, Y np.meshgrid(x, x) img 255 * ((X ** 2 Y ** 2) 1).astype(float) # 添加椒盐噪声 rng np.random.default_rng(42) noise rng.random(img.shape) img_noisy img.copy() img_noisy[noise 0.01] 0 # 随机黑点 img_noisy[noise 0.99] 255 # 随机白点 plt.figure(figsize(10, 5)) plt.subplot(1, 2, 1) plt.imshow(img, cmapgray) plt.title(Original) plt.axis(off) plt.subplot(1, 2, 2) plt.imshow(img_noisy, cmapgray) plt.title(Noisy) plt.axis(off) plt.show()我没有用高斯噪声而选椒盐噪声是因为椒盐噪声在频域里的特征更典型稀疏但强烈的冲激频谱上呈现为散布在整个高频区域的亮点。高斯噪声则是均匀地铺在所有频率上处理起来更麻烦。对教学来说椒盐噪声的滤波效果更明显能让你一眼看到魔法发生。2.3 频域滤波核心代码每个参数都有讲究傅里叶变换去噪的标准流程五步二维FFT、频谱中心化、构造掩膜、频域滤波、逆变换。每一步都有讲究我拆开说。# 1. 二维FFT spectrum fft.fft2(img_noisy) # 2. 频谱中心化把低频移到图像中心 spectrum_shifted fft.fftshift(spectrum) # 3. 构造圆形低通掩膜 rows, cols img.shape crow, ccol rows // 2, cols // 2 radius 40 yy, xx np.ogrid[:rows, :cols] mask ((xx - ccol) ** 2 (yy - crow) ** 2) radius ** 2 # 4. 频域滤波 filtered_spectrum spectrum_shifted * mask # 5. 逆变换取实部 filtered_img fft.ifft2(fft.ifftshift(filtered_spectrum)).real plt.figure(figsize(12, 4)) plt.subplot(1, 3, 1) plt.imshow(np.log(1 np.abs(spectrum_shifted)), cmapgray) plt.title(Spectrum) plt.axis(off) plt.subplot(1, 3, 2) plt.imshow(mask, cmapgray) plt.title(Mask (radius40)) plt.axis(off) plt.subplot(1, 3, 3) plt.imshow(filtered_img, cmapgray) plt.title(Filtered) plt.axis(off) plt.show()为什么要fftshift因为fft2算出来的频谱低频在四个角落高频在中心附近这不符合我们“中心放低频”的视觉习惯。fftshift把频谱四角的低频挪到中心方便画圆掩膜。做完掩膜乘法后必须用ifftshift把频谱还原回原来的排列方式再做逆变换这一步顺序反了会得到一张莫名其妙的花屏。为什么要圆形掩膜而不是方形因为方形掩膜会在滤波后产生明显的振铃效应图像边缘出现一圈圈水波纹。圆形掩膜对应的频域边界是连续旋转对称的振铃会小得多。这是图像处理里的老经验新手经常在这里翻车。radius40这个值不是拍脑袋定的它对应的是保留频率范围。图像尺寸256中心化后最大频率半径约为128。圆盘图的主体是低频半径40已经能覆盖绝大多数主体频率了。你可以试试调成20看效果圆盘边缘会变糊但噪声几乎全没了再调成80边缘清晰了但噪声点会残留不少。这就是上一节说的权衡。2.4 去噪效果评估别只顾着看漂不漂亮最关键的一步是评估效果。我常用的方法是同时输出三张图原图、去噪图、噪声残差图。残差就是原图减去去噪图理想情况下残差里只能看到噪声和一些边缘细节如果残差里出现明显的圆形轮廓说明低通滤波把边缘信息也抹掉了。定量评估可以用信噪比SNR公式是信号能量除以噪声能量。不过对于教学项目肉眼判断反而更直接。我最喜欢的一个做法是放大局部区域对比比如用plt.xlim/ylim截取圆盘边界附近放大能很清楚看到噪声点是否被抹掉、边缘是否变虚。还有一个小技巧低通滤波可以再做一次高通实验来加深理解。把掩膜取反只保留高频再逆变换你会发现得到的图几乎是纯边缘和噪声的叠加。这就能直观印证“细节和噪声都在高频”这个结论。我强烈建议你跑一下这个对比实验比你单独做低通滤波理解要深得多。3. 定积分求圆周率quad函数背后的数值艺术3.1 数学原理半圆面积就是π/2圆周率与定积分的关系本质上是圆的面积公式。单位圆的方程是x² y² 1取上半部分得到y sqrt(1 - x²)。对这个函数在[-1, 1]上做定积分得到的是上半圆的面积也就是π/2。所以π 2 × ∫_{-1}^{1} sqrt(1 - x²) dx。如果你觉得从-1到1麻烦也可以只算第一象限的部分x从0到1积分得到四分之一圆面积π/4最后乘4。两种写法数值上没差别但四分之一圆的写法积分区间更短对数值算法来说通常更友好。这个题目看起来有点“绕”但它的教学价值极高它让你意识到很多看起来和微积分无关的常数都可以通过数值积分算出来。你还可以扩展出很多变体比如用同样的方法算椭圆面积、球体体积只要写出被积函数quad就能处理。3.2 quad函数实操一行代码算出半圆面积直接上代码这是整个项目里最短的一段但信息量一点不少。from scipy.integrate import quad import numpy as np def half_circle(x): return np.sqrt(1 - x ** 2) area_half, error_estimate quad(half_circle, -1, 1) pi_approx 2 * area_half print(f半圆面积 {area_half:.16f}) print(f误差估计 {error_estimate:.2e}) import math print(fπ 近似值 {pi_approx:.16f}) print(fπ 真实值 {math.pi:.16f}) print(f绝对误差 {abs(pi_approx - math.pi):.2e})我跑这段代码得到的结果是半圆面积输出1.5707963267948966π近似值3.141592653589793误差在10的负16次方量级。这个精度已经接近双精度浮点数的极限了相当惊人。重点看quad的第二个返回值error_estimate这是算法自己估算的绝对误差。数值积分最妙的地方在于它不但给你结果还附带一个可信度指标。实际使用中我会习惯性把error_estimate打出来看一眼如果这个误差估计很大说明被积函数可能有问题比如奇异点、震荡剧烈或者区间选择不当。3.3 为什么quad这么准自适应算法揭秘也许你会想不就是求积分吗我自己把区间切成10000份用矩形法累加也能算出差不多的π吧。没错但效率差太多了。我做过一个简单测试用矩形法计算同样的积分要达到误差10的负8次方需要大约50万个月。而quad内部采用自适应Gauss-Kronrod积分它在函数变化缓的地方少采样变化剧的地方多采样。举例来说sqrt(1 - x²)在x接近±1时函数曲线几乎垂直导数趋向无穷。普通均匀采样在端点附近会很吃亏因为这一段函数变化太快均匀采样的密度根本不够。quad的聪明之处就在这里它会先在整体上采样一组点然后检查哪些区域的误差贡献大再对误差大的区域递归加密采样直到整个区间都满足精度要求。这种“把力气花在刀刃上”的策略正是自适应算法的精髓。理解了这一点你就能解释为什么同样的积分quad耗时不到一毫秒精度却碾压手写矩形法几个数量级。3.4 扩展用quad做更复杂的圆面积计算既然标题是“定积分求圆周率”我建议你做完基本版后顺手做两个扩展练习。第一个扩展是算任意半径R的圆面积。只需要在half_circle里加入半径参数用quad(half_circle, -R, R, args(R,))传入。这里args的用法很实用后续做参数扫描时会经常用到。def circle_upper(x, R): return np.sqrt(R ** 2 - x ** 2) for R in [1, 2, 0.5]: area, err quad(circle_upper, -R, R, args(R,)) print(f半径{R}的圆面积: {2 * area:.10f}, 误差估计: {err:.1e})第二个扩展是修改被积函数算半个椭圆的面积。椭圆方程x²/a² y²/b² 1上半部分y b * sqrt(1 - x²/a²)积分区间[-a, a]。结果应该是πab除以2。这个扩展能让你真正理解“定积分求面积”的本质边界函数求积分得到的是曲线下的面积如果把边界函数改一下就能算各种稀奇古怪的图形面积。4. 实战中的坑与排查技巧4.1 FFT谱图里那些奇怪的亮点和条纹第一次把频谱用imshow显示出来的时候很多人会吓一跳为什么中心一个超亮的大斑点四周还有横竖十字形的亮线中心大斑点其实是直流分量也就是全图亮度的平均值。如果是一个全黑背景的圆盘图直流分量会很高这是正常的。横竖亮线通常来自图像边缘的不连续性图像左右边界、上下边界在傅里叶眼里都是突变所以会在对应的频率轴上产生一条亮线。处理这类问题一个经典做法是给图像做边缘平滑或者加窗。不过在我们这个去噪项目里并不需要特别处理。只要你清楚这些现象的来源就不会被它们带偏思路。4.2 逆变换结果出现复数怎么办调用fft.ifft2之后返回结果其实是复数数组。这是因为浮点运算误差理论上实部才是真实图像虚部会有极小的残留。很多人忘了取实部直接plt.imshow复数数组结果要么报错要么显示出一片乱码。正确做法是调用.real取实部代码里我已经写了千万别省。还有一个常见坑是把ifftshift放在fftshift前面。正确的顺序是fft2-fftshift- 滤波 -ifftshift-ifft2。如果你把ifftshift和fftshift搞混唯一的表现就是输出图像错乱这时候先检查这一步。4.3 quad函数报错和异常情况处理quad最常见的报错是ValueError: Inf or NaN这通常是被积函数在积分区间内出现无穷大。比如你把被积函数写成1 / sqrt(1 - x**2)在x ±1处分母为零函数值无穷大但恰恰这个积分是收敛的值等于π。quad其实能处理一部分端点奇异的情况但需要你显式设置points参数把奇异点指出来。# 端点奇异但积分收敛的例子 integral, err quad(lambda x: 1 / np.sqrt(1 - x**2), -1, 1, points[-1, 1]) print(integral) # 接近 pi设置了points之后quad会在这些点附近加密采样从而正确处理奇异行为。这个知识点在很多工程计算里都会用到值得记一下。另一个我遇到过的问题是积分区间过大导致被积函数在中间出现剧烈震荡。这时候quad可能报“Maximum number of subdivisions exceeded”的警告解决办法是把大区间拆成多个小区间分别积分再把结果相加。实际工程里遇到这类问题永远优先考虑分区间而不是盲目调limit参数。4.4 常见问题速查表症状可能原因解决办法频谱中心亮点巨大图像直流分量过高正常现象可做对数变换再显示滤波后图像有波纹使用了方形掩膜换成圆形掩膜降低振铃逆变换输出乱码ifftshift位置错误或没取实部检查顺序加.real去噪后图像过糊滤波半径太小增大radius重新尝试quad报Inf或NaN被积函数含奇异点用points参数标示奇异点quad耗时长且警告积分区间过大或震荡剧烈手动拆分积分区间分别计算这张表是我把实际踩过的坑整理出来的。做图像处理时遇到视觉异常第一反应先看频谱和掩膜形状对不对做数值积分时遇到报错第一反应先把被积函数图像画出来。把函数可视化很多问题一眼就能定位。我个人在这次实践中最大的收获是真正建立了“频域思维”以前遇到图像噪声第一反应是换滤波器窗口大小现在会先想到去频谱里看一眼噪声长什么样。这种思维方式的转变远比记住几个函数调用有价值。至于定积分求π看起来像个玩具但当你意识到quad能处理带奇异点的函数时你就会发现工程里很多复杂定积分都可以放心交给它前提是你要理解算法适应的边界在哪里。希望这两个小项目能帮你打开科学计算的大门。