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