Fourier deconvolution with numpy

Question

I am attempting to remove my probes function from a signal using Fourier deconvolution, but I can not get a correct output with test signals.

t = np.zeros(30)
t = np.append(t, np.arange(0, 20, 0.1))

sigma = 2
mu = 5.
g = 1/np.sqrt(2*np.pi*sigma**2) * np.exp(-(np.arange(mu-3*sigma,mu+3*sigma,0.1)-mu)**2/(2*sigma**2))

def pad_signals(s1, s2):
    size = t.size +g.size - 1
    size = int(2 ** np.ceil(np.log2(size)))
    s1 = np.pad(s1, ((size-s1.size)//2, int(np.ceil((size-s1.size)/2))), 'constant', constant_values=(0, 0))
    s2 = np.pad(s2, ((size-s2.size)//2, int(np.ceil((size-s2.size)/2))), 'constant', constant_values=(0, 0))
    return s1, s2

def decon_fourier_ratio(signal, removed_signal):
    signal, removed_signal = pad_signals(signal, removed_signal)
    recovered = np.fft.fftshift(np.fft.ifft(np.fft.fft(signal)/np.fft.fft(removed_signal)))
    return np.real(recovered)

gt = (np.convolve(t, g, mode='full') / g.sum())[:230]
tr = decon_fourier_ratio(gt, g)

fig, ax = plt.subplots(nrows=2, ncols=2, sharex=True)
ax[0,0].plot(np.arange(0,np.fft.irfft(np.fft.rfft(t)).size), np.fft.irfft(np.fft.rfft(t)), label='thickness')
ax[0,1].plot(np.arange(0,np.fft.irfft(np.fft.rfft(g)).size), np.fft.irfft(np.fft.rfft(g)), label='probe shape')
ax[1,0].plot(np.arange(0,gt.size),gt, label='recorded signal')
ax[1,1].plot(np.arange(0,tr.size),tr, label='deconvolved signal')
plt.show()

enter image description here

The above script creates a demo sample (t), and a probe with Gaussian shape (g). Then, it convolves them to a signal gt, which is what a sample would look like when probed. I pad the signal to the nearest 2^N with pad_signals(), for efficiency and to fix any non-periodicity. Then I try to remove the gaussian probe with decon_fourier_ratio(). As is clear from the images, I do not recover the initial thickness gradient. Any ideas why the deconvolution is not working?

Note: I have also tried SciPy's deconvolve. But, this function only works for gaussians of certain widths.

Any help is greatly appreciated,

Eric

Jacques Kvam · Answer 1 · 2018-04-15T00:33:50.320

1

Any reason you are not doing the full convolution? If I change the construction of gt to:

g /= g.sum()  # so the deconvolved signal has the same amplitude
gt = np.convolve(t, g, mode='full')

Then I get the following plots:

I can't quite tell you why your seeing this behavior, other than the partial convolution is probably altering the frequency content. Alternatively, you can pad your input signal with zeros if you want to get the same behavior and use same.

edited Apr 15 '18 at 00:33

answered Apr 15 '18 at 00:15

Jacques Kvam

2,856
1
26
31

I just tried this does fix the issue of recovering an incorrect signal. My guess is because the convolved signal is continuous now. The only issue is that the thickness is "infinitely increasing", ie my probe will only ever see increasing thickness, so the convolved signal should not decrease. I have posted a corrected image above and corrected the code to crop at 230. – eric hoglund Apr 15 '18 at 19:04

Fourier deconvolution with numpy

1 Answers1