4
$\begingroup$

I am implementing the exponential/logarithmic sine sweep approach to impulse response measurement (as presented by Farina - http://pcfarina.eng.unipr.it/Public/Papers/134-AES00.PDF).

My query is regarding amplitude correction of the inverse sweep. I have found a lot of differing information about what correction should be applied, and the only method that seems to work as I expect, doesn't have any information on how it was derived.

So far, I have generated the sweep (in Python) as follows:

f1 = 20     # Start frequency
f2 = 20000  # Stop frequency
T = 1       # Duration in seconds
Fs = 48000  # Sample rate

L = T/np.log(f2/f1)      # Constant used in sweep calculation

t = np.arange(0,T,1/Fs)  # Time ramp

sweep = np.sin(2*np.pi*f1*L*(np.exp(t/L) - 1))

Then to generate the inverse sweep, you have to time reverse the sweep and then compensate for the fact that, due to the non-linear variation in frequency, the sweep contains less energy at higher frequencies. In terms of the exact compensation to apply, I have found varying information.

In Farina's paper, he says to generate the inverse sweep, you must time reverse the sweep and then apply an amplitude correction that starts at 0 dB and ends at $-6\log{}_2(\frac{f2}{f1})$ dB.

In their work on "Synchronized Swept-Sine", Novak et al. (https://hal.science/hal-02504286/document) give the amplitude correction as $\frac{f1}{L}\exp{}(-\frac{t}{L})$.

In the answer to a question on calculating the inverse of the exponential sine sweep on here (Calculating the inverse filter for the (exponential) sine sweep Method) jojeck gives the correction as $\frac{1}{\exp(\frac{tR}{T})}$ where $R = \log(\frac{f2}{f1})$.

Finally, in another question on here (Impulse response amplitude : Sine sweep method), the OP, bouaaah, gives the correction as (np.exp(-time/L))/L*f2*T**2 but doesn't say where the equation comes from.

Now, I should say at this point that what I am aiming for is that, when the sweep is deconvolved directly with its inverse, the magnitude of the resulting frequency response is 0 dB in the frequency range of interest, which I think is a reasonable thing to expect, but am open to being corrected.

Here is the code to apply all the above corrections and compare results:

rev_sweep = np.flip(sweep)

env = np.linspace(0,-6*(np.log(f2/f1)/np.log(2)),rev_sweep.size)
env = 10**(env/20)

inv_sweep = rev_sweep*env

N = sweep.size + inv_sweep.size

sweepFFT = fft.fft(sweep, N) * 2 / sweep.size
invFFT = fft.fft(inv_sweep, N) * 2 / inv_sweep.size

H = sweepFFT*invFFT
H = H[0:int(H.size/2 + 1)]
f = np.linspace(0,Fs/2,H.size)

fig, ax = plt.subplots()
fig.set_tight_layout(True)
ax.semilogx(f,20*np.log10(np.abs(H)),label='Farina')
ax.grid(True)

env = (f1/L)*np.exp(-t/L)

inv_sweep = rev_sweep*env

N = sweep.size + inv_sweep.size

sweepFFT = fft.fft(sweep, N) * 2 / sweep.size
invFFT = fft.fft(inv_sweep, N) * 2 / inv_sweep.size

H = sweepFFT*invFFT
H = H[0:int(H.size/2 + 1)]
f = np.linspace(0,Fs/2,H.size)

ax.semilogx(f,20*np.log10(np.abs(H)),label='Novak et al.')

R = np.log(f2/f1)
env = 1/np.exp((t*R)/T)

inv_sweep = rev_sweep*env

N = sweep.size + inv_sweep.size

sweepFFT = fft.fft(sweep, N) * 2 / sweep.size
invFFT = fft.fft(inv_sweep, N) * 2 / inv_sweep.size

H = sweepFFT*invFFT
H = H[0:int(H.size/2 + 1)]
f = np.linspace(0,Fs/2,H.size)

ax.semilogx(f,20*np.log10(np.abs(H)),label='jojeck')

env = np.exp(-t/L)/L*f2*T**2

inv_sweep = rev_sweep*env

N = sweep.size + inv_sweep.size

sweepFFT = fft.fft(sweep, N) * 2 / sweep.size
invFFT = fft.fft(inv_sweep, N) * 2 / inv_sweep.size

H = sweepFFT*invFFT
H = H[0:int(H.size/2 + 1)]
f = np.linspace(0,Fs/2,H.size)

ax.semilogx(f,20*np.log10(np.abs(H)),label='bouaaah')
ax.legend()

plt.show()

And this results in:

Different inverse sweep corrections

As you can see, the only correction that gives 0 dB in the passband is the one in bouaaah's question (Impulse response amplitude : Sine sweep method). It should be noted that Farina and jojeck's approaches give very similar results so it's difficult to see the separate lines.

So my question is twofold:

  1. Firstly, is it reasonable to expect that when you deconvolve the sweep directly with its inverse, you should get 0 dB in the passband?
  2. Secondly, if so, can anyone explain where bouaaah's equation comes from? Have I missed any references etc.?
$\endgroup$
5
  • 2
    $\begingroup$ 🤯 looks like there is a bug in my code! Thank you for spotting it! The resulting IR should be a delta function with an amplitude of 1 (or -3 dBFS RMS across in frequency domain). I will go back and fix my math/code. $\endgroup$
    – jojeck
    Aug 7, 2023 at 13:40
  • 1
    $\begingroup$ This paper gives another correction, which is different from bouaaah's method but is correct as well. You can check it out. $\endgroup$
    – ZR Han
    Aug 16, 2023 at 2:41
  • $\begingroup$ A bit late to the party. I am struggling a bit to understand what is going on here. It is an exciting topic and I’d very much like to see it solved as I am also very interested in it. I tried to replicate the results and got the behaviour you described. On the other hand, the impulse returned by the method has a humongous peak, which scales the estimated/measured system response wildly. I may be wrong or missing something here, but I believe that having your deconvolved signal in the ideal case be centred at $0 \, \textrm{dB}$ is not very helpful. You may have another reason to do (cont.) $\endgroup$
    – ZaellixA
    Jan 9 at 21:42
  • $\begingroup$ (cont.’ed) so but I can’t think of a practical advantage of that. However, I didn’t manage to get a scaling that would result in an impulse with unity (close to unity ‘cause some energy will be “spilt” to the ripples next to the peak) and none of the methods provides that (or at least I didn’t manage to achieve it, I may have screwed up somehow). Would you like to elaborate a bit on why you’d like to achieve the $0 \, \textrm{dB}$ average or what advantages this may provide? $\endgroup$
    – ZaellixA
    Jan 9 at 21:45
  • $\begingroup$ Nevermind, I found the answer in your answer to the question providing the “correct” scaling you cited. $\endgroup$
    – ZaellixA
    Jan 10 at 11:18

0

Your Answer

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

Browse other questions tagged or ask your own question.