为何用Scipy的FFT逆卷积恢复指数衰减信号会引入大量噪声？

我用Scipy的fft和ifft函数实现指数衰减信号与高斯仪器响应函数（IRF）的卷积：

convolved = ifft(fft(decay) * fft(gaussian)).real

但尝试用以下代码恢复原始指数衰减信号时：

decay_recovered = ifft(fft(convolved) / fft(gaussian)).real

得到的恢复信号噪声远高于原始信号。请问这是什么原因？能否通过处理卷积信号的FFT结果来避免该问题？

附图：红色为原始衰减信号

完整复现代码如下：

# 导入必要库
import numpy as np
import pandas as pd
from scipy.fftpack import fft, ifft
import random
import matplotlib.pyplot as plt

def create_decay(x,I,t,q):
    """根据x轴值(x)、初始荧光强度(I)、寿命分布中心(t)和异质性参数(q)创建幂律荧光衰减曲线"""
    decay = ((2-q)/t)*(1-((1-q)*(x/t)))**(1/(1-q))
    decay = decay * I
    return decay

def convolve_IRF(irf,decay):
    """对仪器响应函数(IRF)和荧光衰减曲线进行卷积"""
    conv = ifft(fft(decay) * fft(irf)).real
    return conv

def create_gaussian_curve(x):
    """根据x轴值(x)创建高斯曲线，中心值(mu)和sigma值随机生成"""
    sig = round(random.uniform(0.07,0.10),2)
    mu = round(random.uniform(1.10,1.40),2)
    gauss = np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))
    return gauss

def add_Poisson_noise(decay):
    """为曲线添加泊松噪声"""
    noise_mask = np.random.poisson(decay)
    decay = decay + noise_mask
    return decay

def create_decay_real(x,irf):
    """创建真实荧光衰减曲线：输入x轴值(x)和IRF，随机生成初始强度(I)、寿命中心(t)、异质性参数(q)，
    卷积后添加泊松噪声并归一化"""
    I = random.randint(200,1000)
    t = round(random.uniform(1.250,8.000),3)
    q = round(random.uniform(0.950,1.990),3)
    if q != 1.000:
        decay_raw = create_decay(x,I,t,q)
        decay_conv = convolve_IRF(irf,decay_raw)
        decay_noisy = add_Poisson_noise(decay_conv)
        decay_noisy = decay_noisy/np.max(decay_noisy)
    return decay_raw, decay_noisy, t, q

# 创建x轴值
x = np.arange(0.04848485,49.6,0.09696969)

# 创建IRF
irf = create_gaussian_curve(x)

# 创建衰减曲线
decay_raw, decay, t, q = create_decay_real(x,irf)

# 计算恢复原始衰减所需的FFT结果
fft_decay = fft(decay)
fft_irf = fft(irf)

# 恢复原始衰减曲线
pure_decay = ifft(fft_decay / fft_irf).real

# 归一化两条曲线
decay_raw = decay_raw / np.max(decay_raw)
pure_decay = pure_decay / np.max(pure_decay)

# 绘图
plt.plot(x,pure_decay,'k--',x,decay_raw,'r')
plt.savefig('test.png')

核心原因 #

1. IRF频域幅值的高频衰减：高斯信号的频谱呈钟形分布，高频段幅值会趋近于0。直接用卷积信号的FFT除以IRF的FFT时，高频处的极小幅值会将噪声大幅放大（相当于除以一个接近0的数），最终在时域表现为严重噪声。
2. 泊松噪声的全频段分布：你给卷积信号添加的泊松噪声是全频段分布的，高频噪声被IRF的小FFT值放大后，对恢复信号的干扰远大于原始信号。

解决方法：频域正则化/滤波 #

通过对IRF的FFT结果做处理，避免除以过小的值，从而抑制噪声放大。以下是两种实用方案：

方案1：截断小幅值的IRF FFT分量 #

设定一个阈值，当IRF的FFT幅值小于阈值时，跳过除法操作（将分母设为1），避免放大高频噪声：

# 计算FFT结果
fft_decay = fft(decay)
fft_irf = fft(irf)

# 设置阈值（取IRF FFT最大幅值的1%）
threshold = np.max(np.abs(fft_irf)) * 0.01
# 替换过小的分母，避免除以0或极小值
fft_irf_safe = np.where(np.abs(fft_irf) > threshold, fft_irf, 1)

# 恢复信号
pure_decay = ifft(fft_decay / fft_irf_safe).real

方案2：Tikhonov正则化（更稳健） #

给IRF的FFT幅值平方添加一个小正则化参数，避免分母过小，这是Wiener滤波的简化形式，能在信号保真和噪声抑制间取得平衡：

# 计算FFT结果
fft_decay = fft(decay)
fft_irf = fft(irf)

# 设置正则化参数（取IRF FFT最大幅值的0.01倍）
lambda_reg = np.max(np.abs(fft_irf)) * 0.01
# 正则化处理分母，同时用共轭项优化信号恢复
fft_recovered = fft_decay * np.conj(fft_irf) / (np.abs(fft_irf)**2 + lambda_reg**2)

# 恢复信号
pure_decay = ifft(fft_recovered).real

优化后完整代码示例（以Tikhonov正则化为例） #

# 创建x轴值
x = np.arange(0.04848485,49.6,0.09696969)

# 创建IRF
irf = create_gaussian_curve(x)

# 创建衰减曲线
decay_raw, decay, t, q = create_decay_real(x,irf)

# 计算FFT结果
fft_decay = fft(decay)
fft_irf = fft(irf)

# Tikhonov正则化处理
lambda_reg = np.max(np.abs(fft_irf)) * 0.01
fft_recovered = fft_decay * np.conj(fft_irf) / (np.abs(fft_irf)**2 + lambda_reg**2)

# 恢复原始衰减曲线
pure_decay = ifft(fft_recovered).real

# 归一化两条曲线
decay_raw = decay_raw / np.max(decay_raw)
pure_decay = pure_decay / np.max(pure_decay)

# 绘图
plt.plot(x,pure_decay,'k--',x,decay_raw,'r')
plt.savefig('test_denoised.png')

你可以通过调整正则化参数lambda_reg或阈值大小优化结果：参数越大，噪声抑制越强，但信号细节损失越多；参数越小，信号细节保留越好，但噪声可能更明显。

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