Get the room impulse response by using inverse filter of sin sweep signal

Question

I want to simulate the process of obtaining the impulse response of a room.

What I did:

1. Generate a room impulse response sequence by using the python module: rir-generator
2. Generate a sine sweep signal with this formula $$ f(t) = A\sin\left(\frac{2 \pi f_{1}T}{\ln\left( \frac{f_{2}}{f_{1}} \right)}\left(e^{\frac{t}{T} \ln(f_{2}/f_{1})} - 1\right)\right) $$
3. Calculate the inverse filter by using the method under this question: Calculating the inverse filter for the (exponential) sine sweep Method
4. Simulate the response of a room by using scipy.signal.lfilter
5. Do convolution between inverse filter and the response in step 4

Here's my code: (My code may look a bit strange, as I will eventually convert it to c program, so I would like to use a more straightforward expression rather than a more intuitive one)

import numpy as np
import rir_generator
from scipy import signal as si
import matplotlib.pyplot as plt


def sineSweep(fmin, fmax, duration, sampling_rate, peak):
    k = np.exp(np.log(fmax / fmin) / duration)
    data_len = duration * sampling_rate
    sinsw = np.zeros(data_len)  # array of output 
    dt = 1.0 / sampling_rate  # time between samples
    t = 0.0  # start time
    p = 2 * np.pi * fmin / np.log(k)
    for i in range(data_len):
        sinsw[i] = peak * np.sin(p * (pow(k, t) - 1))
        t += dt
    return sinsw


def inverseFilt(sinsw, fmin, fmax):
    length = len(sinsw)
    frac = 1
    k = 1 / (np.exp(np.log(fmax / fmin) / length))

    invf = np.zeros(length)

    for i in range(length):
        invf[i] = sinsw[-i] * frac
        frac *= k

    return invf


# step 1
rir = rir.generate(
    c=340,                  # Sound velocity (m/s)
    fs=48000,           # Sample frequency (samples/s)
    r=[1, 0.75, 0.5],          # Receiver position(s) [x y z] (m)
    s=[1, 1.75, 1],          # Source position [x y z] (m)
    L=[8, 8, 3],            # Room dimensions [x y z] (m)
    reverberation_time=0.4, # Reverberation time (s)
    nsample=1024,           # Number of output samples
).T[0]

f1 = 20
f2 = 20000
fs = 48000

# step 2
sinsw = sineSweep(f1, f2, 1, fs, 1)
# step 3
invf = inverseFilt(sinsw, f1, f2)

plt.subplot(2, 1, 1)
plt.plot(np.arange(48000), sinsw)
plt.title('Sine sweep singal')
plt.subplot(2, 1, 2)
plt.plot(np.arange(48000), invf)
plt.title('Inverse filter')
plt.show()

# step 4
output = si.lfilter(rir, 1, sinsw)
# step 5
ir = si.fftconvolve(output, invf, mode='same')

plt.subplot(2, 1, 1)
plt.plot(np.arange(1024), rir)
plt.title('Room impulse response')
plt.subplot(2, 1, 2)
plt.plot(np.arange(1024), ir[24001: 24001 + 1024])
plt.title('Impulse response by inverse filter')

plt.show()

The rir signal has a length of 1024. The sine sweep signal sinsw, its reverse filter invf and the output of lfilter output has a length of 48000, so as the ir calculated by fftconvolve. I noticed that the segment of ir[24001, 24001 + 1024] is familiar with the rir signal, so I put them together.

Here's What I get:

What bothers me is that even though the two arrays have similar trends, the orders of magnitude are very different. I presume that there is a problem with my inverse filter design; however, most of the articles are vague about the inverse filter design, and I couldn't find any specific formulation.

Moreover, the final impulse response seems to be more "jagged" than the rir. Is it possible to make them less different by modifying the inverse filter?

Also, at a sampling rate of 48,000 Hz, generating a 1 second sine sweep signal would produce 48,000 points. I didn't even add enough silence time to calculate the tail of the impulse response. I think 48,000 points is just too much for a normal room, especially when my simulation input is only 1024 points. When I look at the output, I find that only the 24001st to 25025th points are what I need. Is it possible for me to save complexity by getting only the specified length of data starting from point 24001 and not the rest?

Answer 1

A few things are happening here:

1. Your original impulse response is quite "ringy", i.e. the individual components are non-causal.
2. 1024 is also very short for a reverb time of 0.4s.
3. To simplify things, I would start with a wire (i.e. $H(z) =1 $) and then a very simple filter (1st order lowpass). Once that looks good, you can move on to real impulse responses
4. The inverse is only an approximate inverse. It makes no claim about scaling at all.
5. Start by simply cascading the sine sweep and it's inverse: that's the effective transfer function or your measurement process. Ideally it's flat (i.e. $H(z)$) but it isn't. It's bandpass, with about 18 dB of gain and some high frequency ringing. The latter is probably responsible for your "jaggedness".