I want to simulate the process of obtaining the impulse response of a room.
What I did:
- Generate a room impulse response sequence by using the python module:
rir-generator
- 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) $$
- Calculate the inverse filter by using the method under this question: Calculating the inverse filter for the (exponential) sine sweep Method
- Simulate the response of a room by using
scipy.signal.lfilter
- 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.
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?