Deconvolution of system response in Python/Matlab

Question

I had two sets of data, the output function of the system (time series with a length of 1292 entries) and the transfer function (similar to a gaussian with a length of 681 entries). I would like to calculate the input function(unknown) using deconvolution.

I tried calculating the FFT for each function and then the inverse of the division of the two ffts: Cin = F^-1 [F(Cout)]/[F(E)] however FFT of the transfer function seems to be zero everywhere and the resulting recovered signal have a lot of noise. I also tried Wiener deconvolution, but once again, if I calculate the convolution of the recovered signal with the transfer function I don't obtain the output function.

For the signals I'm using which is the most appropriate method? (FFT, SciPy convolve, Wiener)

Do I need to "zero-pad" the transfer function to match the length of the output function?

Since the deconvolution process is susceptible to noise, Do I need to filter the signals (before/after)?

there is similar code in matlab or python that I could refer to?

Here the output of my code with the wiener method (output in blue, transfer in red, recovered input in green and convolution of recovered input and transfer (purple):

Answer 1

I´m not sure if this really applies to your problem since it may be another issue, but I can tell from my experience with acoustic impulse response measurements (only there you want to estimate the response from the input signal, but the deconvolution should be the same):

In case your transfer function E is somehow bandlimited and has zero (or very low) energy in some frequencies, the inverted magnitude response 1/F(E) will contain infinitely large peaks in those frequencies. For example the inverted version of a simple bandpass magnitude response (fc_low=200Hz, fc_high=18kHz) will boost everything below 200Hz and above 18kHz:

If your Cout-Signal contains any noise (e.g. quantization-noise or from the analog-digital converter) in those frequencies, this noise will be boosted.

My recommendation: Plot and look at the magnitude responses of F(E) and F(Cout). If necessary, apply a filter to the inverse of F(E) before multiplying it with F(Cout). This filter should be designed to attenuate the problematic frequencies.

Example in Python:

from numpy.fft import rfft, irfft
from scipy import signal

Nfft = 2048

# lowpass
sos_lp = signal.butter(8, 18000, 'lowpass', fs=fs, output='sos')
lp_filter = signal.sosfilt(sos_lp, signal.unit_impulse(N_fft))
lp_filter = rfft(lp_filter)

# highpass
sos_hp = signal.butter(8, 200, 'highpass', fs=fs, output='sos')
hp_filter = signal.sosfilt(sos_hp, signal.unit_impulse(N_fft))
hp_filter = rfft(hp_filter)

bandlimit_filter = hp_filter * lp_filter

# deconvolution
FE = rfft(E, Nfft)
FCout = rfft(Cout, Nfft)
FE_inv = 1 / FE
Cin = irfft(FCout * bandlimit_filter * FE_inv, Nfft)
Cin = [:(1292-681+1)] # truncate output

Answer 2

Cin = F^-1 [F(Cout)]/[F(E)] however FFT of the transfer function seems to be zero everywhere and the resulting recovered signal have a lot of noise.

What you're building is effectively an equalizer, in this case, a zero-forcing equalizer. That comes with noise amplification: Where your F(E) is small, the division leads to large numbers – exactly where there is the least actual energy.

For the signals I'm using which is the most appropriate method?

Try a different equalizer! MMSE is the classic choice here; the usual MMSE derivations assume white input signal, and that doesn't describe your signal, but seeing it's relatively broadband, this might still work. You might want to subtract the mean at the receiver first, though.