傅立叶反卷积在 python 上无法按预期工作

Question

我创建了 (A) 从 -100 到 +100 的 bin 数组和 (B) 两个随机数直方图，它们遵循均值 0 和标准差 10 的正态分布。它们用于使用傅立叶变换执行逆卷积 ((C ）和（D））。

使用的代码如下所示：

import numpy as np 
import matplotlib.pyplot as plt 


# (A) create bin_edges 
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) create two gauss_histograms 
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) raw fourier deconvolution 
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) wiener fourier deconvolution 
# I don't know how to determine 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 和标准差 10 的两个高斯函数进行反卷积理想情况下会产生 0 处的 delta 函数，但我的代码却没有。我尝试了维纳反卷积，但结果相似，而且我一开始不知道如何确定 snr 的值。

我还确认，如果我对observed_signal和implse_response使用通用直方图（例如observed_signal = hist1，implse_response = hist1），结果将符合预期。我认为两个直方图的统计波动是原因。我想过使用窗口函数来削减两个直方图的高频成分，但我不知道，因为我还没有实现它（代码也包含在下面）。怎样才能得到想要的结果？

# frequency domain 
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() 

Answer 1

由于随机数据集中的统计差异而产生的两个直方图之间的差异意味着它们的傅里叶变换的相位将不相同。这意味着当您尝试执行以下操作时：

F_recovered = F_observed / F_implse

对于反卷积，不同的阶段会将事情搞乱。在理想情况下，数据集之间没有任何差异，

F_recovered

将仅包含复杂的

1

，然后将其反转为 delta 函数。但是，由于相位混乱，您会收到很多噪音。

正如您上一个问题的评论中所提到的，并且在尝试进行维纳反卷积时，您需要某种形式的正则化。在维纳反卷积尝试中，您可以通过计算信号和噪声功率来估计 SNR，例如：

s = np.abs(np.fft.fft(hist1))**2
n = np.abs(np.fft.fft(hist1 - hist2))**2

snr = np.clip((s / n), -1e9, 1e9)  # clip inf values in case of divide by zero

当您运行此代码时，您会得到：

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()

给出：

正如您所希望的那样，峰值为零。