含随机噪声时Scipy与FFT反卷积的异常问题及解决方案咨询

一、振动分量的成因 #

• 反卷积的本质是频域除法：信号频域 = (观测信号频域) / (核频域)。当核的频域分量接近0时，噪声的微小扰动会被剧烈放大（除以趋近于0的数），这些被放大的高频噪声分量转换回时域后，就表现为额外的振动伪影。
• 随机噪声是全频域分布的，信噪比越低，核频域零点/极小值处的噪声放大效应越显著，最终的振动也就越明显。

二、Scipy与手动FFT反卷积表现不同的原因 #

• Scipy内置的反卷积函数（比如scipy.signal.deconvolve）默认做了隐式的优化：要么对频域中极小值分量做了阈值过滤，要么在时域做了边界补全/截断，一定程度上抑制了噪声放大。
• 手动直接执行ifft(fft(观测信号)/fft(核))的操作，没有任何正则化处理，会完全保留所有频域的噪声放大结果，振动自然更突出。另外，两种方法的边界延拓方式（比如是否补零到2的幂次、补零位置）、相位处理逻辑不同，也会导致结果差异。

三、含噪声场景下的代码修改方案 #

核心思路是加入频域正则化抑制噪声放大，同时统一边界处理逻辑，让两种方法的结果更稳定且符合预期。

方案1：Wiener滤波（经典带噪反卷积） #

Wiener滤波通过引入信噪比参数，平衡信号还原和噪声抑制：

import numpy as np
from scipy.fft import rfft, irfft

def wiener_deconvolve(obs_signal, kernel, noise_power):
    # 转频域计算
    obs_fft = rfft(obs_signal)
    kernel_fft = rfft(kernel)
    # 计算Wiener滤波核
    kernel_conj = np.conj(kernel_fft)
    signal_power = np.mean(np.abs(obs_fft)**2)
    wiener_kernel = kernel_conj / (np.abs(kernel_fft)**2 + noise_power / signal_power)
    # 反卷积转回时域
    deconv_fft = obs_fft * wiener_kernel
    return irfft(deconv_fft)

使用时，noise_power可以通过观测信号的无信号区域方差估算，噪声越大，该值设置越大。

方案2：Tikhonov正则化（截断极小值） #

直接给核的频域分量加一个小常数，避免除以接近0的数：

import numpy as np
from scipy.fft import fft, ifft

def tikhonov_deconvolve(obs_signal, kernel, reg_param=1e-3):
    # 补零到相同长度，避免循环卷积
    target_len = len(obs_signal) + len(kernel) - 1
    obs_padded = np.pad(obs_signal, (0, target_len - len(obs_signal)), 'constant')
    kernel_padded = np.pad(kernel, (0, target_len - len(kernel)), 'constant')
    
    obs_fft = fft(obs_padded)
    kernel_fft = fft(kernel_padded)
    # 加正则项后做除法
    deconv_fft = obs_fft / (kernel_fft + reg_param)
    # 取实部（输入为实数信号）
    return np.real(ifft(deconv_fft))

reg_param需根据噪声强度调整，噪声越大，参数值要适当增大。

方案3：Scipy工具加后处理 #

用Scipy内置函数结合平滑处理，抑制振动：

from scipy.signal import deconvolve, wiener
import numpy as np

# 基础反卷积+滑动平均平滑
deconv_result, _ = deconvolve(obs_signal, kernel)
deconv_smoothed = np.convolve(deconv_result, np.ones(5)/5, mode='same')

# 直接用Scipy的Wiener滤波
noise_std = np.std(obs_signal[:100])  # 假设前100个点是纯噪声
deconv_wiener = wiener(obs_signal, mysize=5, noise=noise_std)

内容的提问来源于stack exchange，提问作者omooon