[FIXED] Fourier deconvolution doesn't work as expected on python

Issue

I created (A)array of bins from -100 to +100 and (B)two histograms from random numbers following a normal distribution with mean 0 and standard deviation 10. They were used to perform inverse convolution using the Fourier transform ((C) and (D)).

The code used is shown below:

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

I would think that deconvolution two Gaussians with mean 0 and standard deviation 10 would ideally yield a delta function at 0, but my code does not. I tried wiener deconvolution, but with similar results, and I didn't know how to determine the value of snr in the first place.

I have also confirmed that if I use a common histogram for observed_signal and implse_response (e.g. observed_signal = hist1, implse_response = hist1), the results are as expected. I believe that the statistical fluctuations in the two histograms are the cause. I thought about using a window function to cut the high-frequency component of the two histograms, but I don't know because I haven't implemented it (The code are also included below). How can I get the desired results?

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

Solution

The differences between your two histograms that come about due to the statistical differences in you random data sets will mean that the phases of their Fourier transforms will not be the same. This means that when you try to do:

F_recovered = F_observed / F_implse

for the deconvolution, the different phases will scramble things up. In ideal circumstances, without any differences between your datasets, F_recovered would just contain complex 1s, which would then invert into a delta function. But, with the scrambled phases you'll get a lot of noise.

As mentioned in the comment to your previous question, and with your attempt at Wiener deconvolution, you need some form of regularisation. In your Wiener deconvolution attempt you could estimate the SNR by working out the signal and noise power, e.g.,:

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

which when you run through this code you get:

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

gives:

which peaks at zero as you would hope.

Answered By - Matt Pitkin