Impulse response amplitude : Sine sweep method

Question

I'm trying to implement Farina's method in python to measure impulse response with an exponential sine sweep.

I convolved the output of the system with the inverse filter of the exponential sweep.

Here is my code for the generation of the exponential sweep and the inverse filter:

def logsweep(fmin,fmax,duration,phi0=0):

    global fs

    time = np.arange(0,duration,1/fs)
    omegaMin = fmin * 2 * np.pi
    omegaMax = fmax * 2 * np.pi

    phi = (omegaMin*duration)/np.log(omegaMax/omegaMin) * (np.exp(time/duration * np.log(omegaMax/omegaMin)) - 1) + phi0 # Phase instantanée

    return np.sin(phi)

def inverseFilter(fmin,fmax,duration,phi0=0):

    global fs
    sweep = logsweep(fmin,fmax,duration,phi0)
    time = np.arange(0,duration,1/fs)
 
    L = 1/fmin*(duration*fmin/(np.log(fmax/fmin)))
    normalisation = (np.exp(-time/L))/L*fmax*duration**2
    sweep_reverse = sweep[::-1] * normalisation
 
 
    return sweep_reverse

When I convolved the output of my system with the inverse fitler, the resulting IR has a good shape but is wrong in amplitude.

Here is my code :

    print("Génération du filtre inverse...")
    invFilter = inverseFilter(fmin,fmax,duration)

    print("Calcul de la IR...")
    farinaTempIR = np.convolve(measure,invFilter,mode="full") / fs
    farinaTempIR = farinaTempIR[sweep.size-1:]
    timeFarinaTempIR = np.arange(0,farinaTempIR.size,1) / fs
    plt.figure("IR")
    plt.plot(timeFarinaTempIR,farinaTempIR,'--',label = "Farina Temp")
    plt.legend()

    print("Calcul de la TF...")
    farinaTempTF = np.fft.rfft(farinaTempIR)
    freqFarinaTempTF = np.fft.rfftfreq(farinaTempIR.size,1/fs)
    plt.figure("Fonction de transfert")
    plt.semilogx(freqFarinaTempTF,20*np.log10(np.abs(farinaTempTF)),'--',label = "Farina Temp")
    plt.legend()

I suppose my inverse filter is correct because, when I did the same operation in frequency domain I got the expected result.

Here is my code in frequency domain:

    print("Génération du filtre inverse...")
    invFilter = inverseFilter(fmin,fmax,duration)

    print("Calcul de la FFT du signal de sortie...")
    measureFFT = np.fft.rfft(measure) / measure.size * 2
    freqMeasure= np.fft.rfftfreq(measure.size,1/fs)

    print("Calcul de la FFT du filtre inverse...")
    invFFT = np.fft.rfft(invFilter) / invFilter.size * 2
    freqInv = np.fft.rfftfreq(invFilter.size,1/fs)

    print("Calcul de la TF...")
    farinaFreqTF = measureFFT * invFFT
    plt.figure("Fonction de transfert")
    plt.semilogx(freqInv,20*np.log10(np.abs(farinaFreqTF)),'--',label = "Farina Freq")
    plt.legend()

    print("Calcul de la IR...")
    farinaFreqIR = np.fft.irfft(farinaFreqTF)
    timeFarinaFreqIR = np.arange(0,farinaFreqIR.size,1) / fs
    plt.figure("IR")
    plt.plot(timeFarinaFreqIR,farinaFreqIR,'--',label = "Farina Freq")
    plt.legend()

Here is the plots of my code:

I suppose that I'm missing something about the discrete convolution operator. Should I divide by the sampling frequency? By the sum of the inverse filter? I'm not really confortable with this operator.

Answer 1

All that you have missed is scaling of the sweep and inverse sweep prior to convolution because you have scaled their responses in the frequency domain.

The convolution theorem states that the convolution of two signals in the time domain is equivalent to the inverse Fourier transform of the product of the Fourier transforms of the two signals:

$$r(s) = \{g*h\}(s) = \mathcal{F}^{-1}\{G \cdot H\}$$

(Wikipedia: Convolution theorem)

We can check that this is correct (without any scaling at this stage):

sweep = logsweep(fmin,fmax,duration)

inv_sweep = inverseFilter(fmin,fmax,duration)

irTimeDomain = np.convolve(sweep,inv_sweep)

N = sweep.size + inv_sweep.size - 1

sweepFFT = np.fft.fft(sweep, N)
invFFT = np.fft.fft(inv_sweep, N)

irFreqDomain = np.real(np.fft.ifft(sweepFFT*invFFT))

fig, ax = plt.subplots(3,1)
fig.set_tight_layout(True)
ax[0].plot(irTimeDomain)
ax[0].grid(True)
ax[1].plot(irFreqDomain)
ax[1].grid(True)
ax[2].plot(irTimeDomain-irFreqDomain)
ax[2].grid(True)

We have to make the FFT of each signal N samples long because the convolution in the time domain results in a signal that is sweep.size + inv_sweep.size - 1 long. This gives the following (top time domain, middle frequency domain, bottom difference):

As you can see, the difference between the two is basically precision.

Now, returning to your question, to get equivalence between your time and frequency domain implementations, you simply need to scale the time domain sweep and inverse by (2 / their length) prior to convolution:

irTimeDomain = np.convolve(sweep * 2 / sweep.size, inv_sweep * 2 / 
inv_sweep.size)

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

irFreqDomain = np.real(np.fft.ifft(sweepFFT*invFFT))

fig, ax = plt.subplots(3,1)
fig.set_tight_layout(True)
ax[0].plot(irTimeDomain)
ax[0].grid(True)
ax[1].plot(irFreqDomain)
ax[1].grid(True)
ax[2].plot(irTimeDomain-irFreqDomain)
ax[2].grid(True)

Which, we can again confirm, results in the same thing:

And if we take the Fourier transform of the resulting impulse responses, we get a 0 dB response across the frequency range of interest:

H = np.fft.fft(irFreqDomain)
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)))
ax.grid(True)

Now, may ask you a question? Where did you get the scaling equation for your inverse sweep? I haven't seen it anywhere in the literature and none of the other scaling equations result in a 0 dB frequency response. See my own question here exponential-logarithmic-sine-sweep-inverse-filter-amplitude-correction.