On the behavior of deconvolution with scipy and fft with random noise

85 Views Asked by At

In the 1D histogram, I verified whether deconvolution of the same two signals, kernel and signal, with random noise on them produces the expected distribution in python.

The code used is as follows

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import deconvolve


def gaussian(x, mean, std_dev, noise_value):
    gaus = (1.0 / (std_dev * np.sqrt(2 * np.pi))) * np.exp(-0.5 * ((x - mean) / std_dev)**2)
    random_noise = np.random.uniform(-noise_value, noise_value, size=len(x))
    
    # The range outside of 3 standard deviations is set to zero
    # Add random noise to the range inside 3 standard deviations
    gaus[(x < mean - 3 * std_dev) | (x > mean + 3 * std_dev)] = 0
    gaus[(x >= mean - 3 * std_dev) & (x <= mean + 3 * std_dev)] += random_noise[(x >= mean - 3 * std_dev) & (x <= mean + 3 * std_dev)]
    return gaus


def scipy_deconvolution(kernel, signal):
    first_idx = np.argmax(kernel != 0)
    end_idx = first_idx + np.argmax(kernel[first_idx:] == 0)
    recover, remainder = deconvolve(signal, kernel[first_idx:end_idx])
    padleft = (len(remainder) - len(recover)) // 2
    padright = padleft + (len(remainder) - len(recover)) % 2
    scipy_recovered = np.pad(recover, (padleft, padright), mode='constant')
    return scipy_recovered


def fft_deconvolution(kernel, signal):
    fft_signal = np.fft.fft(signal)
    fft_kernel = np.fft.fft(kernel)
    fft_recover = fft_signal / fft_kernel
    recover = np.fft.fftshift(np.fft.ifft(fft_recover))
    fft_recovered = np.real(recover)
    return fft_recovered


# histogram parameter
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)
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2


mean = 0
std_dev = 20
noise_value = 0.000001
kernel  = gaussian(bin_centers, mean, std_dev, noise_value)
signal  = gaussian(bin_centers, mean, std_dev, noise_value)


# Each Recoverd distribution is expected to match the signal distribution
scipy_recovered = scipy_deconvolution(kernel, signal)
fft_recovered = fft_deconvolution(kernel, signal)


# plot
fig, axes = plt.subplots(nrows=1, ncols=4, sharex=True, figsize=(12, 3))
titles = ['Kernel', 'Signal', 'Scipy Deconv', 'FFT Deconv']
for i, ax in enumerate(axes.flat):
    ax.hist(bin_edges[:-1], bin_edges, weights=[kernel, signal, scipy_recovered, fft_recovered][i])
    ax.set_title(titles[i])
plt.tight_layout()
plt.show()

Without noise, the scipy and ft results are as shown in the figure below and the signal is as expected.

enter image description here

My goal is to create a code that produces results similar to the above figure when random noise is included in the data to be deconvolved. However, increasing the value of random noise yields different results than expected (see figure below) and also behaves differently in the deconvolution of scipy and fft method. enter image description here enter image description here

Question: I would like to know what causes the extra vibration component in the presence of random noise and why the difference in behavior between the two methods. And it would be great if you could tell me how to fix the current code to get the expected results with both scipy and fft methods even in situations involving random noise.

0

There are 0 best solutions below