Estimating the Impulse Response of the Room Using Sweep Signal Microphone Recorded Signal (Input & Output of a Convolution)

Question

I played this signal A (a 20Hz to 20000Hz sinusoidal sweep in 10 seconds) with a studio monitor speaker in a big church, and I recorded the result B with good microphones.

The result is very reverb-ish, that's exactly what I wanted to catch.

Now a software (such as Deconvolver but non open-source) can build an impulse response from A + B, that can be later used in a convolution reverb.

It works well. But I would like to learn how to do this myself via DSP / programming.

How can I use a sweep (signal A) + recorded output (signal B) to get an impulse reponse?

---

Edit: In other words, if a is the original sweep, b the recorded output, and h the impulse response, how to get h from

$$a * h = b$$

Is this formulation correct? is the solution $h=a^{-1} * b$, where $a^{-1}$ is the inverse of $a$ for the convolution? How to compute a convolution-inverse of a discrete signal?

Answer 1

this is the two-channel FFT method of spectrum analyzer:

$$ y[n] = h[n] \ \circledast \ x[n] $$

just make sure that the length of the FFT $N$ is at least as large as the length of sound $x[n]$ plus the expected length of the impulse response $h[n]$. the length of sound $y[n]$ is also as long as the FFT. you can round $N$ up to the nearest power of two. just zero-pad everything to that length and then

$$H[k] = \frac{Y[k]}{X[k]} $$

is functionally true.

you might sometimes have to worry about division by zero, but if your driving signal $x[n]$ is sufficiently broad-banded (which a linear sweep or a maximum-length sequence is), then you don't have to worry too much about division by zero.

if you know $H[k]$, then you know $h[n]$.

Answer 2

So here is the math for how to use a pair of linearly-swept sinusoids to extract the impulse response from a room or some ostensibly linear time-invariant (LTI) system having impulse response $h(t)$ and transfer function (or "frequency response"):

$$\begin{align} H(f) &\triangleq \mathscr{F} \big\{h(t) \big\} \\ & = \int\limits_{-\infty}^{\infty} h(t) \, e^{-j 2 \pi f t} \, dt \ \end{align}$$

The "reference" or "driving" or input signal:

$$ \begin{align} x(t) &= e^{j \pi \beta t^2 } \\ &= \cos\left(\pi \beta t^2\right) \ + \ j \sin\left(\pi \beta t^2\right) \\ &= x_\text{re}(t) \ + \ j \, x_\text{im}(t) \\ \end{align}$$

So it's really two driving signals, $x_\text{re}(t)$ and $x_\text{im}(t)$.

One at a time, pass $x_\text{re}(t)$ through the system $H$ having response $y_\text{re}(t)$ followed by $x_\text{im}(t)$ (which is a real signal when the "$j$" is not attached) resulting in response $y_\text{im}(t)$. Both times the response is recorded synchronously with the input signal, so you know when $t=0$ and how to align $y_\text{re}(t)$ and $y_\text{im}(t)$ in the mind of the computer.

$\beta$ is the sweep rate of the sweep in Hz per second (if $t$ is in seconds). $\mathfrak{f}(t) \triangleq \beta t$ is the instantaneous frequency of the sweep at time $t$. You start with the instantaneous frequency at a very large negative value (like -Nyquist), sweep through $\mathfrak{f}(t)=0$ and continue to a very large positive value. Do that for both the cosine sweep and the sine sweep.

Then, in the mind of the computer, you assemble this complex response, $y(t)$ to the complex input, $x(t)$, from the two real responses $y_\text{re}(t)$ and $y_\text{im}(t)$:

$$ y(t) \triangleq y_\text{re}(t) \ + \ j \,y_\text{im}(t) $$

You can do this because the system $H(f)$ is linear and time-invariant.

So then the complex output $y(t)$, even though complex, is simply the response of the LTI system to the complex input $x(t)$.

$$\begin{align} y(t) &= y_\text{re}(t) \ + \ j \,y_\text{im}(t) \\ \\ &= h(t) \circledast x_\text{re}(t) \ + \ j \, h(t) \circledast x_\text{im}(t) \\ \\ &= \int\limits_{-\infty}^{\infty} h(u) \, x_\text{re}(t-u) \, du \ + \ j \int\limits_{-\infty}^{\infty} h(u) \, x_\text{im}(t-u) \, du \\ &= \int\limits_{-\infty}^{\infty} h(u) \, \big( x_\text{re}(t-u) \ + \ j \, x_\text{im}(t-u) \big) \, du \\ &= \int\limits_{-\infty}^{\infty} h(u) \, x(t-u) \, du \quad = \ h(t) \circledast x(t) \\ &= \int\limits_{-\infty}^{\infty} h(u) \, e^{j \pi \beta (t-u)^2 } \, du \\ &= \int\limits_{-\infty}^{\infty} h(u) \, e^{j \pi \beta (t^2-2tu+u^2) } \, du \\ &= \int\limits_{-\infty}^{\infty} h(u) \, e^{j \pi \beta t^2}e^{-j \pi \beta 2 t u}e^{j \pi \beta u^2} \, du \\ &= e^{j \pi \beta t^2} \int\limits_{-\infty}^{\infty} \big( h(u)\,e^{j \pi \beta u^2} \big) \, e^{-j 2 \pi \, \beta t \, u } \, du \\ &= e^{j \pi \beta t^2} \int\limits_{-\infty}^{\infty} \big( h(u)\,e^{j \pi \beta u^2} \big) \, e^{-j 2 \pi \, \mathfrak{f}(t) \, u } \, du \\ \end{align}$$

So now you take your output $y(t)$ (with both real and imaginary parts) and you "adjust" it by multiplying by the complex conjugate of the sweep function

$$\begin{align} \tilde{y}(t) &\triangleq y(t) \cdot e^{-j \pi \beta t^2} \\ &= \int\limits_{-\infty}^{\infty} \big( h(u)\,e^{j \pi \beta u^2} \big) \, e^{-j 2 \pi \, \beta t \, u } \, du \\ \end{align}$$

When you evaluate that "adjusted" time function, $\tilde{y}(t)$ at a time that is a frequency $f$ divided by the sweep rate $\beta$, you have in the time function a representation of the Fourier Transform of something:

$$\begin{align} \tilde{y}(t)\Bigg|_{t=\tfrac{f}{\beta}} &= \int\limits_{-\infty}^{\infty} \big( h(u)\,e^{j \pi \beta u^2} \big) \, e^{-j 2 \pi f u} \, du \\ &= \int\limits_{-\infty}^{\infty} \tilde{h}(u) \, e^{-j 2 \pi f u} \, du \\ &= \tilde{H}(f) \\ \end{align}$$

where $ \tilde{h}(t) \triangleq h(t) \, e^{j \pi \beta t^2} $ and

$$\begin{align} \tilde{H}(f) &\triangleq \mathscr{F} \big\{\tilde{h}(t) \big\} \\ &= \mathscr{F} \big\{h(t)\,e^{j \pi \beta t^2}\big\} \ \end{align}$$

This all means that the "adjusted" response in the time domain follows the frequency response of the "adjusted" LTI system:

$$ \tilde{y}(t) = \tilde{H}(\beta t) = \tilde{H}\big(\mathfrak{f}(t)\big) $$

and that your unadjusted response is another sweep, with the same sweep rate $\beta$ and same instantaneous frequency $\mathfrak{f}(t)=\beta t$, but with magnitude and phase modified by this "adjusted" frequency response:

$$ y(t) = \tilde{y}(t) \, e^{j \pi \beta t^2} = \tilde{H}(\beta t) e^{j \pi \beta t^2} $$

So what you're getting with linear swept frequency measurements are not the frequency response directly of your impulse response

$$ H(f) = \mathscr{F} \big\{h(t) \big\} $$

but the frequency response of your impulse response that is multiplied by the sweep function itself:

$$ \tilde{H}(f) = \mathscr{F} \big\{\tilde{h}(t) \big\} = \mathscr{F} \big\{h(t)\,e^{j \pi \beta t^2}\big\} $$

So there are a couple of things you can do about this.

1. Select your sweep rate $\beta$ to be so slow that this is approximately true: $$ \tilde{h}(t) \approx h(t) $$ for all values of $t$ such that $|h(t)|$ is not close to zero. That means that $$ e^{j \pi \beta t^2} \approx 1 $$ for all $|t|$ small enough that $|h(t)| \gg 0$ and that $$\begin{align} \tilde{y}(t) &= y(t) \, e^{-j \pi \beta t^2} \\ &= \tilde{H}(\beta t) \\ &\approx H(\beta t) \\ \end{align}$$ so your time-domain response represents very closely your frequency response. That's sorta the immediate motivation behind using swept-frequency sinusoids for identifying LTI systems.

or

2. Inverse Fourier Transform $\tilde{H}(f)$ (using the DFT and sufficient padding) to get $\tilde{h}(t)$ and like you did to $y(t)$, multiply $\tilde{h}(t)$ by the complex conjugate of the sweep to get $h(t)$

$$ h(t) = \tilde{h}(t) \, e^{-j \pi \beta t^2} $$

Answer 3

If the sweep is constructed well, then convolving then its frequency spectrum shows a nearly constant magnitude in the relevant band. This means that if you multiply the Fourier transform of the sweep with its complex conjugate, you get a zero phase signal with nearly constant magnitude in that same band. Interpreted as system responses, the sweep response is therefore undone by the time-reversed sweep response resulting in the identity response for the band of interest.

So in order to undo the sweep response part of your recorded room response, you simply convolve the recorded signal by the time reversed sweep signal to get the impulse response of the room.

I'll add some math if I find the time and there's need for it.