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:
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:
- Firstly, is it reasonable to expect that when you deconvolve the sweep directly with its inverse, you should get 0 dB in the passband?
- Secondly, if so, can anyone explain where bouaaah's equation comes from? Have I missed any references etc.?
