1
$\begingroup$

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?

$\endgroup$

2 Answers 2

3
$\begingroup$

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]
$\endgroup$
2
  • $\begingroup$ do you have any suggestions on how to turn these equations into code? i'm not very experienced with coding and i thought i could use the signal.deconvolve function for this. $\endgroup$ May 15, 2022 at 17:27
  • $\begingroup$ @MertAlperten code added, see edit $\endgroup$
    – ZR Han
    May 16, 2022 at 2:03
0
$\begingroup$

Let's take a look at your transfer function

enter image description here

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.

$\endgroup$
2
  • $\begingroup$ hmm actually this one is just the test impulse response. the impulse responses that i'm going to test for my thesis work probably will be simpler without weird peaks and steep high/low filters. the thing that i'm wondering is that why i can't have any audio results with this code? and what do i need to change for it to work albeit not fully successfully. $\endgroup$ May 15, 2022 at 17:38
  • $\begingroup$ Use a test impulse response that has a reasonably well behaved inverse. $\endgroup$
    – Hilmar
    May 16, 2022 at 12:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.