How to use deconvolution code with an impulse response to achieve the original signal

Question

I am trying to write a code for my thesis to deconvolve a recording with an impulse response so that I can achieve the original audio signal. I have written a simple code so that I can implement this but it is not working. Here is the code:

import numpy as np

import matplotlib.pyplot as plt

import scipy.signal

import scipy.io.wavfile as wavfile

IR_PATH = "case11_ir.wav"

AUDIO_EXCERPT_PATH = "case1_gtr.wav"

#%% Loading the IRs

fs, case_ir_sig = wavfile.read(IR_PATH)

case_ir_sig = case_ir_sig[case_ir_sig != 0]

rec_fs, rec_sig = wavfile.read(AUDIO_EXCERPT_PATH)

assert(fs == rec_fs)

dry_sig, r = scipy.signal.deconvolve(rec_sig, case_ir_sig)

dry_sig = dry_sig / np.max(np.abs(dry_sig))

wavfile.write("case11.wav", fs, dry_sig)

---

The problem seems to be at the deconvolution line, whereas I cannot see the dry_sig in the variable explorer. I have used float 32 format for all wav files. And here are the wav files. case1_gtr is the convolved audio, case11_ir is the impulse response, and case_original is the original recording that I am trying to achieve (at least something similar to it) as the product of deconvolution.

Any suggestions?

Answer 1

scipy.signal.deconvolve performs polynomial division and I don't think it's useful for long sequence. You can try frequency-domain deconvolution. As we know (circular) convolution in time domain is equivalent to multiplication in frequency domain: $$ Y(k) = X(k) H(k) $$ where $Y(k)$, $X(k)$ and $H(k)$ are the $N$-point DFT of output, input and the system impulse response, where $N\geq L_x+L_h-1$. One can derive the DFT of input signal as $$ X(k) = \frac{Y(k)}{H(k) + \lambda(k)} $$ where $\lambda(k)$ is a frequency-dependent regularization parameter in order to avoid division by zero, as your impulse response has a bandpass magnitude response. $\lambda(k)$ should be chosen as zero in the pass band and some small values in the stop band.

Simple code would be like:

L = lenx + lenh - 1
threshold = -30
regu = 1e-5

Y = np.fft.fft(y, L)
H = np.fft.fft(h, L)
lam = [0 if 20*log10(abs(i))> threshold else regu for i in H]
X = Y / (H + lam)
x = np.fft.ifft(X, L)
x = x[:lenx]

Answer 2

Let's take a look at your transfer function

There is a steep highpass at about 1kHz and weird spike at 344 Hz. There are a lot of very narrow dips at higher frequencies. This is way too complicated a filter for straight forward deconvolution.

The best you can do here is to carefully construct an inverse filter manually based the signal to noise ratio of the original recording and the ideal gain of the inverse filter.

It's also important the set the requirements correctly: Do you actually need the original waveform (which would extremely difficult) or something just "sounds as close possible to the the original". The latter one has much higher chance of success.