二维Morlet小波图像去噪：时频局部化降噪原理与工程实现

简介：本资源是一套面向信号处理与图像分析初学者及进阶学习者的Morlet小波实验实践包，聚焦二维Morlet小波在图像多尺度分解与信号去噪中的核心应用。内容涵盖理论原理、MATLAB代码实现、可视化结果与实测数据，适用于数字图像处理、遥感/医学影像预处理、课程设计及科研入门场景。压缩包共9个文件（2.81MB），含4个MATLAB源码（如b.m、imageScaleT.m等，实现三级小波分解与去噪流程）、2张原始/处理后JPG图像、1张PNG效果图、1个MAT文件（7.31.mat，存储实验数据）及1个FIG图形文件，完整呈现从一维信号到二维图像的Morlet小波变换全流程。已有557人学习下载，用户可直接运行代码复现实验，获取带注释的去噪脚本、多尺度系数可视化方法、阈值选取参考及Morlet图像生成逻辑，快速掌握小波去噪的关键参数调优与结果评估技巧。

1. Morlet小波为什么是二维图像去噪的“隐形主力”：它不靠卷积核大小赢，而靠时频局部化赢

你有没有试过用高斯滤波或中值滤波处理一张带纹理的医学CT切片，结果边缘糊成一片、细小血管直接消失？或者在遥感图像里，想压制条带噪声又怕把农田边界也抹平？这时候翻开源码看别人怎么做的，十有八九会撞见morlet——不是作为某个深度学习模块的装饰，而是真正在底层扛起时频分析大旗的实战组合。Morlet小波不是“更高级的滤波器”，它是把图像当成二维非平稳信号来解构：既关心某块区域“能量强不强”（幅值），也死磕“这股能量集中在哪个尺度、哪个方向”（频率+相位）。标题里反复出现的“二维_Morlet图像_信号去噪”，说的就是这件事：用复数Morlet小波在图像平面做连续小波变换（CWT），把噪声和结构分别钉在不同尺度-方向通道里，再做阈值裁剪——这不是图像处理，是信号处理思维在像素阵列上的落地。适合谁？不是只调cv2.bilateralFilter参数的初学者，而是手上有低剂量CT、红外热成像、显微电镜图、SAR遥感图等信噪比吃紧、结构细节敏感的一线算法工程师；也适合正被“传统滤波保边难、深度学习缺标注、小波包分解维度爆炸”三重卡脖子的团队。它不承诺端到端PSNR暴涨5dB，但能给你可解释、可调控、不依赖大数据集的确定性降噪路径。

2. 从一维Morlet到二维Morlet：为什么不能直接把1D公式套进图像？

2.1 一维Morlet小波的“血统”与局限：复指数+高斯窗的物理直觉

Morlet小波本质是一个复指数载波被高斯窗调制的结果。标准一维形式为：

$$ \psi(t) = \pi^{-1/4} e^{i \omega_0 t} e^{-t^2 / 2} $$

其中 $\omega_0$ 是中心角频率（通常取5~6以保证时频分辨率平衡），$\pi^{-1/4}$ 是归一化系数。关键点在于：它是个复函数，输出包含实部（cosine-like）和虚部（sine-like），合起来能同时捕获信号的幅度和相位信息。这对一维信号去噪极有用——比如心电图R波检测，相位突变比幅值变化更鲁棒。但直接把它当卷积核在图像上滑动？会出大问题。原因有三：

1. 各向同性陷阱：1D Morlet沿时间轴延展，但图像有x/y两个空间维度。若简单用 $ \psi(x) \cdot \psi(y) $ 做可分离乘积，得到的是圆对称小波，无法区分水平边缘、垂直纹理、45°裂缝——而真实图像结构高度方向敏感；
2. 尺度耦合失效：1D中缩放参数 $a$ 控制单一尺度，但在2D中，仅缩放x/y相同倍数（各向同性缩放）会丢失“长条状噪声”（如CT扫描线）的定向抑制能力；
3. 相位信息冗余：图像灰度是实值场，1D Morlet的复输出在2D中会产生四组冗余分量（实/虚 × x/y），徒增计算且无物理意义。

提示：别被“二维小波”字面迷惑——真正有效的2D Morlet不是1D的简单外积，而是构造方向选择性的复数基函数。这是所有后续操作的起点。

2.2 二维Morlet小波的工程化定义：方向+尺度+偏移三要素

工业界和论文中广泛采用的2D Morlet定义（如Torrence & Compo, 1998）是：

$$ \psi_{a,\theta}(x,y) = \frac{1}{a^2} \pi^{-1/2} e^{i \omega_0 \left( \frac{x \cos\theta + y \sin\theta}{a} \right)} e^{-\left[ \left( \frac{x \cos\theta + y \sin\theta}{a} \right)^2 + \left( \frac{-x \sin\theta + y \cos\theta}{a} \right)^2 \right] / 2} $$

这个式子看着吓人，拆解后就是三个可控旋钮：

• 尺度参数 $a$：控制小波在主方向（$\theta$）上的伸展长度。$a$ 越大，感受野越宽，对应低频（粗结构）；$a$ 越小，聚焦越细，对应高频（噪声/边缘）。实践中 $a$ 取 2^k 形式（k=0,1,2,...）形成对数尺度序列；
• 方向参数 $\theta$：决定小波的“朝向”。$\theta=0^\circ$ 捕捉水平结构，$\theta=90^\circ$ 捕捉垂直结构，$\theta=45^\circ$ 捕捉斜向纹理。典型设置为 $\theta \in {0^\circ, 45^\circ, 90^\circ, 135^\circ}$，共4个方向；
• 旋转坐标系：式中 $x \cos\theta + y \sin\theta$ 是沿 $\theta$ 方向的投影（主轴），$-x \sin\theta + y \cos\theta$ 是垂直方向（副轴）。高斯窗在主轴方向按 $a$ 缩放，在副轴方向也按 $a$ 缩放——这是各向同性缩放；若要各向异性（如拉长副轴以增强线状特征），需额外引入副轴缩放因子 $b$（$b \neq a$），但会显著增加参数调优成本，Morlet图像去噪中95%场景用各向同性已足够。

2.3 在Python中手搓二维Morlet小波核：避开SciPy的坑

很多工程师第一反应是查scipy.signal.morlet2，但注意：morlet2返回的是1D Morlet在指定尺度下的采样，不是2D核！它设计初衷是给1D信号做CWT，强行reshape成2D会得到错误的方向响应。正确做法是自己生成2D网格并代入公式：

import numpy as np import matplotlib.pyplot as plt def morlet2d(shape, scale, theta, omega0=5.0): """ 生成二维Morlet小波核 :param shape: (height, width) 输出核尺寸，建议为奇数（如33x33） :param scale: 尺度参数 a > 0 :param theta: 方向角（弧度） :param omega0: 中心频率，默认5.0（保证时频局部化） :return: 复数二维数组 (H, W) """ h, w = shape # 创建中心对齐的坐标网格（-h//2 到 h//2-1） y = np.arange(-h//2, h//2).reshape(-1, 1) # (h, 1) x = np.arange(-w//2, w//2).reshape(1, -1) # (1, w) # 旋转坐标系：u = x*cosθ + y*sinθ, v = -x*sinθ + y*cosθ u = x * np.cos(theta) + y * np.sin(theta) v = -x * np.sin(theta) + y * np.cos(theta) # Morlet公式：π^(-1/2) * exp(i*ω0*u/a) * exp(-(u²+v²)/(2a²)) # 注意：这里省略了1/a²归一化（因后续做卷积时会由conv2d自动处理） psi = (np.pi**(-0.5) * np.exp(1j * omega0 * u / scale) * np.exp(-(u**2 + v**2) / (2 * scale**2))) return psi # 示例：生成一个33x33、尺度a=4、方向0°的Morlet核 kernel_0deg = morlet2d((33, 33), scale=4, theta=0) print(f"Kernel shape: {kernel_0deg.shape}, dtype: {kernel_0deg.dtype}") # 输出：Kernel shape: (33, 33), dtype: complex128

这段代码的关键逻辑说明：

• y和x使用arange(-h//2, h//2)确保核中心在(0,0)，这对保持卷积的空间对齐至关重要；
• u/v的旋转计算必须严格按公式，任何符号错误（如v的负号漏掉）会导致方向响应完全错乱；
• omega0=5.0是经验值：小于4则高斯窗太宽，时域定位差；大于7则复指数振荡过密，频域泄漏严重；
• 返回complex128类型，因为后续CWT需要保留相位信息用于重构。

注意：此核是复数，不能直接用cv2.filter2D（它只支持实数核）。必须用scipy.signal.convolve2d或 PyTorch 的F.conv2d（输入转为复数张量）。

3. 二维Morlet连续小波变换（CWT）实战：如何把一张图变成多尺度-多方向特征图？

3.1 CWT流程图：不是一次卷积，而是“尺度×方向”的全排列扫描

对一张灰度图像 $I(x,y)$ 做2D Morlet CWT，本质是：对每个预设尺度 $a_k$ 和每个预设方向 $\theta_m$，用对应的2D Morlet核 $\psi_{a_k,\theta_m}(x,y)$ 与图像做卷积，得到该尺度-方向下的复数响应 $W_{a_k,\theta_m}(x,y)$。整个过程可理解为构建一个4D张量：(尺度数, 方向数, 高度, 宽度)。例如，取4个尺度（a=2,4,8,16）和4个方向（0°,45°,90°,135°），最终得到16张复数特征图。每张图的模长|W|表示该位置在该尺度-方向下的能量强度，相位angle(W)表示结构走向。去噪的核心就藏在这里：噪声在所有尺度-方向上呈现均匀、无结构的“毛刺”能量，而真实结构只在特定尺度-方向上形成连贯的高能量脊线。

3.2 用Scipy实现高效CWT：避免for循环的向量化技巧

直接写四层嵌套for循环（尺度×方向×图像高×图像宽）会慢到无法忍受。正确姿势是：预生成所有核，堆叠成4D张量，再用scipy.signal.convolve2d批量卷积。但注意：convolve2d不支持批量核，所以得用scipy.ndimage.convolve配合np.stack：

from scipy import ndimage import numpy as np def cwt_2d_morlet(image, scales, thetas, omega0=5.0, kernel_size=33): """ 对图像执行2D Morlet连续小波变换 :param image: 2D numpy array (H, W)，灰度图 :param scales: 尺度列表，如 [2,4,8,16] :param thetas: 方向列表（弧度），如 [0, np.pi/4, np.pi/2, 3*np.pi/4] :return: 4D complex array (len(scales), len(thetas), H, W) """ h, w = image.shape # 预生成所有核并堆叠：(S, T, K, K) kernels = [] for a in scales: for theta in thetas: kernel = morlet2d((kernel_size, kernel_size), scale=a, theta=theta, omega0=omega0) kernels.append(kernel) kernels = np.stack(kernels) # (S*T, K, K) # 将图像扩展为 (1, H, W) 以便广播 image_3d = image[np.newaxis, ...] # (1, H, W) # 批量卷积：对每个核，与图像做2D卷积 # 注意：ndimage.convolve默认用'constant'填充，边界效应需后续处理 cwt_result = np.zeros((len(scales), len(thetas), h, w), dtype=np.complex128) idx = 0 for i, a in enumerate(scales): for j, theta in enumerate(thetas): # 卷积输出与输入同尺寸（mode='same'） conv_real = ndimage.convolve(image, np.real(kernels[idx]), mode='constant', cval=0.0) conv_imag = ndimage.convolve(image, np.imag(kernels[idx]), mode='constant', cval=0.0) cwt_result[i, j] = conv_real + 1j * conv_imag idx += 1 return cwt_result # 示例调用 img = np.random.rand(256, 256) # 模拟含噪图像 scales = [2, 4, 8, 16] thetas = [0, np.pi/4, np.pi/2, 3*np.pi/4] cwt_out = cwt_2d_morlet(img, scales, thetas) print(f"CWT output shape: {cwt_out.shape}") # (4, 4, 256, 256)

这段代码的性能关键点：

• kernel_size=33是经验值：太大（如65）导致核内大部分值趋近于0，纯属算力浪费；太小（如15）则无法覆盖Morlet的有效支撑域（约±3σ），造成截断误差；
• mode='constant', cval=0.0是最稳妥的边界填充，避免reflect或wrap引入虚假周期性；
• 分开计算实部/虚部卷积，是因为ndimage.convolve不支持复数核——这是Scipy的硬限制，绕不开；
• 输出cwt_out[i,j]是复数矩阵，后续所有操作（阈值、重构）都基于其模长和相位。

3.3 可视化CWT结果：看懂“能量脊线”才是去噪的开始

光有数据不够，得会读图。以下代码将CWT结果中某尺度-方向的模长图可视化，并叠加原始图像对比：

def plot_cwt_slice(cwt_result, scale_idx, theta_idx, original_img, title_suffix=""): """绘制单个尺度-方向的CWT模长图""" magnitude = np.abs(cwt_result[scale_idx, theta_idx]) fig, axes = plt.subplots(1, 2, figsize=(12, 5)) # 左图：原始图像 axes[0].imshow(original_img, cmap='gray') axes[0].set_title(f'Original Image {title_suffix}') axes[0].axis('off') # 右图：CWT模长（归一化到0-1） mag_norm = (magnitude - magnitude.min()) / (magnitude.max() - magnitude.min() + 1e-8) im = axes[1].imshow(mag_norm, cmap='jet') axes[1].set_title(f'CWT Magnitude (scale={scales[scale_idx]}, θ={int(np.degrees(thetas[theta_idx]))}°)') axes[1].axis('off') plt.colorbar(im, ax=axes[1], fraction=0.046, pad=0.04) plt.tight_layout() plt.show() # 绘制尺度2、方向0°的响应 plot_cwt_slice(cwt_out, scale_idx=0, theta_idx=0, original_img=img, title_suffix="(noisy)")

观察重点：

• 在干净区域（如均匀背景），模长图呈现低幅值、无规律的“雪花噪点”；
• 在边缘/纹理处，模长图出现连续、高亮的线条（脊线），其走向与边缘方向一致；
• 在噪声密集区（如椒盐噪声点），模长图出现孤立、尖锐的亮点，但无延伸性。
  这就是去噪的判据：保留脊线，抑制孤立点。下一章的阈值策略，全基于这个视觉直觉。

4. 小波系数阈值策略：为什么全局阈值是玄学，而尺度-方向自适应才是正解

4.1 全局阈值的三大翻车现场：它为何在Morlet CWT中彻底失效

很多教程直接套用Donoho的VisuShrink公式：threshold = σ * sqrt(2*log(N))（N为像素总数，σ为噪声标准差）。但在2D Morlet CWT中，这招大概率翻车：

• 现象1：边缘断裂。全局阈值一刀切，把弱边缘（如CT中早期微钙化灶）的脊线能量误判为噪声削掉；
• 现象2：伪影残留。噪声在某些尺度-方向上能量意外地高（如传感器固定模式噪声），全局阈值不够狠，残留条带；
• 现象3：纹理失真。自然纹理（如木材年轮、织物经纬）在多个尺度上都有响应，全局阈值无法区分“结构”和“噪声”的能量分布形态。

根本原因：Morlet CWT的系数统计特性随尺度和方向剧烈变化。小尺度（a=2）下，系数近似高斯白噪声；大尺度（a=16）下，系数呈现长程相关性（结构主导）。用同一阈值处理，等于让小学生和博士生考同一张数学卷。

4.2 尺度-方向自适应阈值：用局部方差估计噪声强度

工业级做法是：对每个尺度 $a_k$ 和每个方向 $\theta_m$，独立估计该通道的噪声标准差 $\sigma_{k,m}$，再计算对应阈值。核心思想是——噪声在CWT域中仍近似白噪声，其方差可用系数的局部统计量估计。常用方法：

• 中位绝对偏差（MAD）法：对|W_{k,m}|的所有像素，计算MAD = median(| |W| - median(|W|) |)，则 $\sigma \approx MAD / 0.6745$；
• 鲁棒中位法：取|W_{k,m}|的低百分位（如第1%）像素值作为噪声基线，再向上浮动2~3倍；
• 我们推荐的混合策略（兼顾鲁棒性与效率）：

def estimate_sigma_per_channel(magnitude_map, method='mad'): """ 为单个CWT通道的模长图估计噪声标准差 :param magnitude_map: 2D array, |W_{k,m}(x,y)| :param method: 'mad' or 'percentile' :return: scalar sigma """ if method == 'mad': # MAD法：对所有像素计算MAD med = np.median(magnitude_map) mad = np.median(np.abs(magnitude_map - med)) sigma = mad / 0.6745 else: # percentile法：取1%分位数，再乘系数 p1 = np.percentile(magnitude_map, 1) sigma = p1 * 2.5 # 经验系数，可根据图像类型微调 return max(sigma, 1e-6) # 防止sigma为0 def adaptive_thresholding(cwt_result, method='mad', threshold_factor=1.2): """ 对CWT结果进行尺度-方向自适应阈值 :param cwt_result: 4D complex array (S, T, H, W) :param method: 阈值估计方法 :param threshold_factor: 阈值放大系数（>1.0） :return: 阈值后的4D complex array """ S, T, H, W = cwt_result.shape cwt_thresh = np.zeros_like(cwt_result) for i in range(S): for j in range(T): mag = np.abs(cwt_result[i, j]) sigma = estimate_sigma_per_channel(mag, method=method) thresh = threshold_factor * sigma # 软阈值（更平滑）：W_thresh = sign(W) * max(|W| - thresh, 0) phase = np.angle(cwt_result[i, j]) mag_thresh = np.maximum(mag - thresh, 0) cwt_thresh[i, j] = mag_thresh * (np.cos(phase) + 1j * np.sin(phase)) return cwt_thresh # 应用自适应阈值 cwt_thresh = adaptive_thresholding(cwt_out, method='mad', threshold_factor=1.2)

参数说明：

• threshold_factor=1.2是起点：太小（1.0）去噪不足，太大（1.5）易伤结构。实际项目中，我们总在1.1~1.3间微调；
• 用软阈值而非硬阈值：软阈值让系数平滑过渡到0，避免硬截断引入吉布斯振铃；
• estimate_sigma_per_channel中method='mad'更鲁棒，'percentile'在强结构图像中更快（因只算分位数）。

4.3 避坑：Morlet CWT去噪的5个致命误区与血泪经验

误区1：直接对复数系数做阈值，忽略相位一致性

• 现象：去噪后图像出现诡异的“彩虹色条纹”或大面积模糊。
• 原因：对复数W = A*exp(iφ)直接W[W < thresh] = 0，破坏了A和φ的耦合关系。当A被置零但φ未同步清零，逆变换时相位混乱。
• 解决：永远先算mag = |W|，阈值作用于mag，再用原φ重建W_thresh = mag_thresh * exp(iφ)。代码中phase = np.angle(...)正是为此。

误区2：CWT后不做系数重构，以为模长图就是去噪结果

• 现象：输出的“去噪图”全是彩色斑点，完全不像原图。
• 原因：CWT系数是中间表示，不是图像。必须通过小波逆变换（ICWT）把阈值后的系数映射回像素域。Morlet的ICWT有解析解，但工程中更常用重构核法（下一章详解）。
• 解决：把cwt_thresh当作新特征图，必须走完整重构流程，不可跳步。

误区3：尺度数量太少（<3）或太多（>8），导致频带覆盖不全

• 现象：小尺度噪声没压住，或大尺度结构（如器官轮廓）被过度平滑。
• 原因：尺度序列应覆盖图像的主要频率成分。太少则频带缺口；太多则计算爆炸且小尺度噪声与大尺度结构混叠。
• 解决：用对数尺度scales = [2**i for i in range(min_power, max_power+1)]。对512x512图，min_power=1, max_power=5（即2,4,8,16,32）是黄金组合。

误区4：方向数固定为4，忽视图像内容特异性

• 现象：处理文字扫描件时，水平/垂直方向效果好，但45°方向全是噪声；处理织物图时，45°/135°方向反而最关键。
• 原因：方向数应与图像主结构方向匹配。通用图选4方向（0/45/90/135），但若已知主方向（如CT扫描线为水平），可精简为2方向（0/90）提速。
• 解决：先用cv2.Canny或梯度直方图粗估主方向，再定制thetas。

误区5：忽略CWT的冗余性，用conv2d后直接拼接，导致内存溢出

• 现象：cwt_out占用GB级内存，程序崩溃。
• 原因：CWT是冗余变换（系数数 > 像素数）。4尺度×4方向×256×256 = 1MB，但若用64尺度×8方向，直接飙到16MB。
• 解决：
  1. 用np.float32存储模长（重构时再转复数）；
  2. 对每个尺度-方向单独处理（不用堆叠4D张量）；
  3. 用dask.array或分块计算（对超大图）。

5. 从CWT系数到去噪图像：Morlet逆变换（ICWT）的两种落地路径

5.1 理论逆变换的困境：为什么Morlet没有完美解析ICWT

Morlet小波不是正交基，也不是双正交基，因此不存在严格的、能量守恒的解析逆变换公式。文献中常写的ICWT积分式：

$$ I(x,y) = \frac{1}{C_\psi} \int_0^\infty \int_{-\infty}^\infty \int_{-\infty}^\infty W_{a,\theta}(x',y') , \psi_{a,\theta}\left(\frac{x-x'}{a}, \frac{y-y'}{a}\right) , dx' dy' \frac{da}{a^3} d\theta $$

其中 $C_\psi$ 是容许性常数。但这个三重积分在离散图像上无法精确实现：

• 连续尺度 $a$ 必须离散化，引入近似误差；
• 方向 $\theta$ 离散化后，旋转核的插值带来失真；
• 数值积分精度受网格密度制约，计算量爆炸。

所以工程实践必须妥协：用重构核（Reconstruction Kernel）替代理论ICWT。核心思想是——既然正向CWT是卷积，那逆变换就该是某种“反卷积”，而Morlet的重构核，就是其自身共轭翻转（conjugate and flip）的归一化版本。

5.2 重构核法（Reconstruction Kernel Method）：稳定、快速、可复现

这是工业界首选方案。步骤清晰：

1. 对每个尺度 $a_k$ 和方向 $\theta_m$，生成其对应的重构核 $g_{a_k,\theta_m}(x,y) = \frac{1}{a_k^2} \psi_{a_k,\theta_m}^*(-x,-y)$；
2. 将阈值后的系数 $W_{a_k,\theta_m}^{thresh}(x,y)$ 与 $g_{a_k,\theta_m}$ 做卷积；
3. 对所有尺度-方向的结果求和，再除以总能量归一化因子。

关键洞察：由于Morlet是复数，其重构核必须是共轭翻转（不是简单翻转）。代码实现：

def reconstruction_kernel_2d(shape, scale, theta, omega0=5.0): """ 生成2D Morlet重构核：g(x,y) = (1/a²) * ψ*(-x,-y) """ h, w = shape y = np.arange(-h//2, h//2).reshape(-1, 1) x = np.arange(-w//2, w//2).reshape(1, -1) # 翻转坐标：-x, -y u_flip = (-x) * np.cos(theta) + (-y) * np.sin(theta) v_flip = -(-x) * np.sin(theta) + (-y) * np.cos(theta) # 共轭：exp(i*...) -> exp(-i*...) psi_conj = (np.pi**(-0.5) * np.exp(-1j * omega0 * u_flip / scale) * np.exp(-(u_flip**2 + v_flip**2) / (2 * scale**2))) # 乘以1/a²归一化 g = (1.0 / (scale**2)) * psi_conj return g def icwt_reconstruct(cwt_thresh, scales, thetas, kernel_size=33, image_shape=None): """ 用重构核法从阈值CWT系数重建图像 :param cwt_thresh: 4D complex array (S, T, H, W) :param image_shape: 原图尺寸 (H, W)，用于初始化输出 :return: 2D real array (H, W) """ S, T, H, W = cwt_thresh.shape if image_shape is None: image_shape = (H, W) # 初始化重建图像 recon_img = np.zeros(image_shape, dtype=np.complex128) # 对每个尺度-方向，生成重构核并卷积 for i, a in enumerate(scales): for j, theta in enumerate(thetas): # 生成重构核 g = reconstruction_kernel_2d((kernel_size, kernel_size), scale=a, theta=theta) # 对该通道系数做卷积（注意：cwt_thresh[i,j]是复数） conv_real = ndimage.convolve(np.real(cwt_thresh[i, j]), np.real(g), mode='constant', cval=0.0) conv_imag = ndimage.convolve(np.imag(cwt_thresh[i, j]), np.imag(g), mode='constant', cval=0.0) conv_complex = conv_real + 1j * conv_imag recon_img += conv_complex # 归一化：除以总能量（经验系数） # 理论上应除以 C_psi，但实践中用均值归一化更鲁棒 recon_img = np.real(recon_img) # 取实部（虚部应接近0） recon_img = (recon_img - recon_img.min()) / (recon_img.max() - recon_img.min() + 1e-8) return recon_img.astype(np.float32) # 执行重构 denoised_img = icwt_reconstruct(cwt_thresh, scales, thetas, kernel_size=33, image_shape=img.shape) print(f"Denoised image shape: {denoised_img.shape}, dtype: {denoised_img.dtype}")

这段代码的生存指南：

• reconstruction_kernel_2d中u_flip/v_flip的推导必须严格，任何坐标符号错误都会导致重构图像整体偏移或模糊；
• ndimage.convolve再次被使用，因为它支持实数核与复数输入的卷积（实部/虚部分开算）；
• 最终np.real(recon_img)是必须的——理论上虚部应为0，但数值误差会残留微小虚部；
• 归一化用(x-min)/(max-min)而非除以C_psi，因为C_psi依赖于连续积分，离散化后无精确值，经验归一化更稳定。

5.3 验证去噪效果：不止看PSNR，更要盯住“结构保真度”

PSNR/SSIM是必要但不充分指标。我们坚持三个验证动作：

1. 残差图可视化：residual = |original - denoised|，理想情况是残差集中在噪声位置，结构区域残差≈0；
2. 频谱对比：对原图、去噪图、残差图分别做2D FFT，看高频噪声是否被压制，而中频结构频谱是否保留；
3. 关键结构ROI放大检查：如CT中的血管分叉点、SAR中的道路交叉口，手动放大100%，确认边缘是否锐利、无振铃、无伪影。

def validate_denoising(original, denoised, title="Denoising Validation"): """三合一验证：残差图、频谱、ROI放大""" residual = np.abs(original - denoised) # 计算FFT（中心化） fft_orig = np.fft.fftshift(np.fft.fft2(original)) fft_deno = np.fft.fftshift(np.fft.fft2(denoised)) fft_res = np.fft.fftshift(np.fft.fft2(residual)) # ROI：取中心64x64区域放大 h, w = original.shape roi_orig = original[h//2-32:h//2+32, w//2-32:w//2+32] roi_deno = denoised[h//2-32:h//2+32, w//2-32:w//2+32] fig, axes = plt.subplots(2, 3, figsize=(15, 10)) # 行1：原图、去噪图、残差图 axes[0,0].imshow(original, cmap='gray'); axes[0,0].set_title('Original') axes[0,1].imshow(denoised, cmap='gray'); axes[0,1].set_title('Denoised') im3 = axes[0,2].imshow(residual, cmap='hot'); axes[0,2].set_title('Residual'); plt.colorbar(im3, ax=axes[0,2]) # 行2：频谱（取log10(abs+1)增强可视性） axes[1,0].imshow(np.log10(np.abs(fft_orig)+1), cmap='viridis'); axes[1,0].set_title('FFT Original') axes[1,1].imshow(np.log10(np.abs(fft_deno)+1), cmap='viridis'); axes <p> <a href="https://download.csdn.net/download/weixin_42696271/25533828" style="color:#ec7500;font-size:14px;"> 本文还有配套的精品资源，点击获取 </a> <img alt="menu-r.4af5f7ec.gif" src="https://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif" style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;"> </p>