Python中傅里叶反卷积未达预期效果的解决求助

背景与当前实现 #

我完成了以下操作：

• (A) 创建范围为-100到+100的分箱数组
• (B) 生成两组服从均值0、标准差10的正态分布随机数，以此得到两个直方图
• (C) 基于傅里叶变换实现原始反卷积
• (D) 尝试维纳反卷积

当前使用的代码如下：

import numpy as np 
import matplotlib.pyplot as plt 

# (A) 创建分箱边界
bin0 = 100 
bin_width = 1 
eps = 1e-10 
min_value = - (bin0 * bin_width + round(bin_width/2, 2)) 
max_value = + (bin0 * bin_width + round(bin_width/2, 2)) + eps 
bin_edges = np.arange(min_value, max_value, bin_width)

# (B) 生成两个高斯直方图
np.random.seed(42) 
random_data1 = np.random.normal(loc=0, scale=10, size=80000) 
random_data2 = np.random.normal(loc=0, scale=10, size=80000) 
hist1, _ = np.histogram(random_data1, bins=bin_edges) 
hist2, _ = np.histogram(random_data2, bins=bin_edges) 

plt.hist(bin_edges[:-1], bin_edges, weights=hist1, alpha=0.5, label='hist1') 
plt.hist(bin_edges[:-1], bin_edges, weights=hist2, alpha=0.5, label='hist2') 
plt.legend() 
plt.show()

# (C) 原始傅里叶反卷积
observed_signal = hist1 
implse_response = hist2 
F_observed = np.fft.fft(observed_signal) 
F_implse = np.fft.fft(implse_response) 
F_recovered = F_observed / F_implse 
recovered_signal = np.fft.ifft(F_recovered) 
shifted_recovered_signal = np.fft.fftshift(recovered_signal) 

plt.hist(bin_edges[:-1], bin_edges, weights=np.real(shifted_recovered_signal), alpha=0.5, label='recovered_signal') 
plt.legend() 
plt.show()

# (D) 维纳反卷积
# 不知道如何确定SNR
snr = 100 
F_wiener = np.conj(F_implse) / (np.abs(F_implse)**2 + 1/snr) 
F_recovered = F_observed * F_wiener 
recoverd_signal = np.fft.ifft(F_recovered) 
shifted_recovered_signal = np.fft.fftshift(recoverd_signal) 

plt.hist(bin_edges[:-1], bin_edges, weights=np.real(shifted_recovered_signal), alpha=0.5, label='wiener_recovered_signal') 
plt.legend() 
plt.show() 

问题描述 #

理论上，对两个参数完全相同的高斯分布进行反卷积，结果应该是位于0点的狄拉克δ函数，但当前代码无法得到该预期结果。尝试维纳反卷积后结果类似，且不知道如何合理设置SNR（信噪比）的取值。

已验证：若使用同一直方图同时作为观测信号和脉冲响应（如observed_signal=hist1、implse_response=hist1），结果符合预期。因此推测问题源于两个独立直方图的统计波动，曾考虑用窗函数截断直方图的高频分量，但尚未实现相关代码（以下是初步的频域分析代码）：

# 频域分析
F_hist1 = np.fft.fft(hist1) 
F_abs1 = np.abs(np.fft.fftshift(F_hist1)) 
F_hist2 = np.fft.fft(hist2) 
F_abs2 = np.abs(np.fft.fftshift(F_hist2)) 

plt.plot(F_abs1, label='freq_hist1') 
plt.plot(F_abs2, label='freq_hist2') 
plt.legend() 
plt.show() 

解决方案 #

1. 问题根源 #

两个独立生成的直方图存在随机统计误差，在频域中表现为高频噪声。原始反卷积直接做除法会放大这些噪声，导致结果偏离δ函数。

2. 窗函数滤波实现 #

通过在频域施加低通窗函数，截断高频噪声分量，抑制统计波动的影响。常用的窗函数有汉宁窗、汉明窗等，这里以汉宁窗为例。

3. 维纳反卷积的SNR确定 #

SNR可通过信号功率与噪声功率的比值计算：

• 信号功率：取直方图低频分量的平均功率
• 噪声功率：取直方图高频分量的平均功率

完整优化代码 #

import numpy as np 
import matplotlib.pyplot as plt 

# (A) 创建分箱边界
bin0 = 100 
bin_width = 1 
eps = 1e-10 
min_value = - (bin0 * bin_width + round(bin_width/2, 2)) 
max_value = + (bin0 * bin_width + round(bin_width/2, 2)) + eps 
bin_edges = np.arange(min_value, max_value, bin_width)
num_bins = len(bin_edges) - 1

# (B) 生成两个高斯直方图
np.random.seed(42) 
random_data1 = np.random.normal(loc=0, scale=10, size=80000) 
random_data2 = np.random.normal(loc=0, scale=10, size=80000) 
hist1, _ = np.histogram(random_data1, bins=bin_edges) 
hist2, _ = np.histogram(random_data2, bins=bin_edges) 

plt.hist(bin_edges[:-1], bin_edges, weights=hist1, alpha=0.5, label='hist1') 
plt.hist(bin_edges[:-1], bin_edges, weights=hist2, alpha=0.5, label='hist2') 
plt.legend() 
plt.title("原始高斯直方图")
plt.show()

# 频域分析与窗函数生成
F_hist1 = np.fft.fft(hist1) 
F_hist2 = np.fft.fft(hist2)
# 生成汉宁窗，调整为低通窗
window = np.hanning(num_bins)
window = np.fft.fftshift(window)
# 保留前90%的有效频率，截断高频噪声
cutoff_ratio = 0.9
cutoff_idx = int(num_bins * cutoff_ratio / 2)
window[:cutoff_idx] = 0
window[-cutoff_idx:] = 0
window = 1 - window  # 反转得到低通窗

# (C) 带窗的原始傅里叶反卷积
F_observed = np.fft.fft(hist1) 
F_implse = np.fft.fft(hist2) 
# 施加窗函数抑制高频噪声，添加小值避免除零
F_observed_windowed = F_observed * window
F_implse_windowed = F_implse * window
F_recovered = F_observed_windowed / (F_implse_windowed + 1e-8)
recovered_signal = np.fft.ifft(F_recovered) 
shifted_recovered_signal = np.fft.fftshift(recovered_signal) 

plt.hist(bin_edges[:-1], bin_edges, weights=np.real(shifted_recovered_signal), alpha=0.5, label='带窗反卷积结果') 
plt.legend()
plt.title("带窗原始反卷积结果")
plt.show()

# (D) 优化的维纳反卷积
# 计算实际SNR
F_abs = np.abs(np.fft.fftshift(F_implse))
# 分离低频信号分量与高频噪声分量
signal_bins = np.concatenate([F_abs[:int(num_bins*0.1)], F_abs[-int(num_bins*0.1):]])
noise_bins = F_abs[int(num_bins*0.4):int(num_bins*0.6)]
signal_power = np.mean(signal_bins**2)
noise_power = np.mean(noise_bins**2)
snr = signal_power / noise_power if noise_power > 0 else 100

# 构建维纳滤波核并施加窗函数
F_wiener = np.conj(F_implse) / (np.abs(F_implse)**2 + noise_power/signal_power)
F_wiener = F_wiener * window
F_recovered = F_observed * F_wiener 
recoverd_signal = np.fft.ifft(F_recovered) 
shifted_recovered_signal = np.fft.fftshift(recoverd_signal) 

plt.hist(bin_edges[:-1], bin_edges, weights=np.real(shifted_recovered_signal), alpha=0.5, label='维纳反卷积结果') 
plt.legend()
plt.title("优化维纳反卷积结果")
plt.show()

# 可视化频域窗函数
plt.plot(np.fft.fftshift(window), label='低通窗函数')
plt.legend()
plt.title("频域低通窗函数")
plt.show()

结果说明 #

1. 带窗的原始反卷积通过截断高频噪声，有效抑制了统计波动的影响，结果会更接近0点的δ函数。
2. 维纳反卷积通过计算实际SNR，结合窗函数滤波，能进一步优化反卷积效果，减少噪声放大。
3. 可通过调整cutoff_ratio参数，控制保留的频率分量比例，找到最优的噪声抑制与信号保真平衡点。

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