5
$\begingroup$

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?

$\endgroup$
6
  • $\begingroup$ This was posted on the DSP subreddit a few days ago; it evaluates a few methods to do what you want: tarabah.me/roomresponse1 $\endgroup$
    – MBaz
    Dec 28, 2016 at 15:35
  • $\begingroup$ do you want this for the general deconvolution of extracting $h$ out of $b$ and $a$? (that's the 2-channel FFT.) or do you want this for swept-frequency sinusoid? $\endgroup$ Dec 29, 2016 at 0:23
  • $\begingroup$ @robertbristow-johnson A general solution would be interesting. (I'm playing right now with ifft(fft(b) / fft(a)) as suggested by your comment, and it's quite good btw!) $\endgroup$
    – Basj
    Dec 29, 2016 at 0:26
  • $\begingroup$ looks like you've succeeded. :-) $\endgroup$ Jan 3, 2017 at 6:57
  • 1
    $\begingroup$ Here is another solution (Least Means Square estimate of the channel given a known input and received signal) which is specifically solving inverse convolution: dsp.stackexchange.com/questions/31318/… Just swap input and output to get the channel instead of the channel compensator! $\endgroup$ Nov 28, 2019 at 13:46

3 Answers 3

1
$\begingroup$

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]$.

$\endgroup$
12
  • $\begingroup$ Thanks. Why do you call this "Two-channel FFT method", and why "Spectrum analyzer", is there a link with something else? $\endgroup$
    – Basj
    Dec 29, 2016 at 0:47
  • $\begingroup$ there is some nice math for this specifically when $$ x[n] = e^{j \tfrac12 \beta n^2} $$ turns out that $X[k]$ is self-similar with different parameters. it takes both a cosine sweep and a sine sweep (both with synchronized recording of the response) and you do the complex math in the mind of the computer. but there is no division by zero if $X[k]$ is an exponential of something. $\endgroup$ Dec 29, 2016 at 0:48
  • $\begingroup$ i dunno. it's just what the guys i would bump into at AES events would call that. it's one method (maybe the most common) of three methods, each have been used multiple times. the other two are the linear sweep method (in which $x[n]$ is the sweep signal above) and the maximum-length sequence method (in which $x[n]= (-1)^{a[n]} = \pm1$ and $a[n] \in \{0,1\}$ is a special pseudo-random bit sequence). both are wide-banded and division by zero is not likely to happen $\endgroup$ Dec 29, 2016 at 0:54
  • 1
    $\begingroup$ the cool thing about 2-channel FFT is that you can use music as a driving signal and the spectrum analyzers have some alg which will deal with the denominator if it gets too close to zero. maybe $$ H[k] = \frac{Y[k]}{X[k]} = \frac{Y[k] \, X^*[k]}{|X[k]|^2} \approx \frac{Y[k] \, X^*[k]}{|X[k]|^2 + \epsilon} $$ for a very tiny $$ \epsilon \ll \ll 1$$ but the cool thing is that you don't have to make the audience listen to some nasty chirp or noise signal and still get an idea of what the room is because the music is usually sufficiently broad-banded. $\endgroup$ Dec 29, 2016 at 0:57
  • $\begingroup$ they call it "2-channel FFT" simply because it had 2 input channels and you could run simultaneous FFT on each channel and you could look at the magnitude of the FFT of either channel or the two on top of each other (like different colors) or the FFT of one divided by the other (which tells you something about the LTI system including and in-between the loudspeakers and microphone. $\endgroup$ Dec 29, 2016 at 1:04
1
$\begingroup$

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

  1. 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} $$

$\endgroup$
7
  • $\begingroup$ Thanks @robert-bristow-johnson. This looks great, I'll need to read it a few times to get the idea (I'm an amateur DSP-ist). In short, is it another method than dsp.stackexchange.com/a/36561/5648 (a better one? if so, in which sense? better SNR?) or just a motivation for using a sweep? Maybe it would be interesting to add a few words to explain the general idea of the method for future reference. Thanks again for your great answers. $\endgroup$
    – Basj
    Jan 3, 2017 at 9:09
  • $\begingroup$ yes, it is an alternative to the 2-channel FFT technique and to the Maximum Length Sequences technique (and to the obvious repeated impulse train technique). note that this sweep is 4 times longer than your sweep. it is two sweeps (a cosine sweep and a sine sweep) and it is double-sided (theoretically sweeps from $-\infty$ through zero to $+\infty$). but this is a very analytical method and you will never divide by zero. it is presented in continuous time because the math is cleaner. $\endgroup$ Jan 3, 2017 at 18:21
  • $\begingroup$ BTW, the best way to make the cosine sweep and sine sweep exact in phase (where it is important), is to make the sweep go from 0 to +Nyquist (that's the second half of the sweep) and then mirror image and copy the sound file for the first half. the cosine sweep does not flip the polarity in the flipping and the sine sweep does flip polarity when you flip it. $\endgroup$ Jan 3, 2017 at 18:25
  • $\begingroup$ Ok thanks. Should it give better results than the method I used? (better in which meaning?) $\endgroup$
    – Basj
    Jan 3, 2017 at 18:25
  • $\begingroup$ it's mathematically complete and it doesn't divide by zero. but the 2-channel FFT method is okay if you deal wisely with the possibility of division by zero. (you don't have to do it this way if you don't want. you might have difficulty aligning $y_\text{re}[n]$ and $y_\text{im}[n]$. maybe not. are you able to play back a recording into the speakers and sample at the very same time with the very same device? $\endgroup$ Jan 3, 2017 at 18:30
0
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Thanks @Jazzmaniac. I added some math at the end of my question. Is my formulation $a * h = b$ correct? PS: yes I'd be very interested if you can add some math in the answer, definitely! PS2: typo: then convolving then its $\endgroup$
    – Basj
    Dec 28, 2016 at 14:54

Your Answer

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

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