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

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

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.

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?

• As already stated in another answer, your IR has a very short length. Generally speaking, allowing for a long IR vector does no harm in most (if not all) situations (for LTI systems). The worst thing that could happen is to get a whole bunch of zeros (or very small values) at a big part of your vector which you could easily discard in the post-processing. Mar 24 at 17:12
• The result in this post of a properly compensated chirp with FFT for channel estimation may be useful to you so wanted share: dsp.stackexchange.com/a/66545/21048 Mar 24 at 22:37

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
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". 