1
$\begingroup$

I calculated the transfer function of a system in Frequency Domain by division of Y and X. Since I only deal with exponential sweep sines (ESS) as excitation, I want to reproduce the transfer function in Time Domain by convolution of y with the inverse of x. I get both spectra plotted, but it seems there is a amplitude scaling missing, since the whole spectrum of time convolution is shifted by ~+60dB. When I plot the spectra normalized to 1 kHz, the plots are the same in the frequency range of excitation. What am I missing?

Thank you very much for your help.

Here are the plots of the Transfer function and the normalized transfer function 1 2.

As well as plots from the Impulse Response for better understanding for the Spectral Division Aproach and for Convolution with Inverse filter approach 3 4

import numpy as np import scipy.signal as sig import matplotlib.pyplot as plt # ESS Sweep Signal f1 = 1 f2 = 22000 T = 8 T_sig = 6 T_silence = 2 fs = 48000 t_sig = np.arange(0, T_sig*fs)/fs t = np.arange(0, T*fs)/fs R = np.log(f2/f1) x = np.sin((2*np.pi*f1*T_sig/R)*(np.exp(t_sig*R/T_sig)-1)) x = 0.25*x # Inverse filter with equalized, time reversed x # # Angelo Farina's Apporach to calculate the impulse response # http://pcfarina.eng.unipr.it/Public/Presentations/AES122-Farina.pdf # https://aes2.org/publications/elibrary-page/?id=14106 k = np.exp(t_sig*R/T_sig) x_inv = x[::-1]/k silence_samples = int(T_silence * fs) # Add silence x = np.concatenate([x, np.zeros(silence_samples)]) x_inv = np.concatenate([np.zeros(silence_samples), x_inv]) # Recording #y = record(x, fs) y = azpf_out # response of the system # Time Domain # convolution with inverse filter ir_convolve = sig.convolve(y, x_inv, mode='full') # Crop impulse response, keep the first half and excess the last part ir_convolve_truncated = ir_convolve[int(len(x_inv)):] plt.figure(1) plt.title('Impulse Response (Convolution with inverse x)') x_plot = np.arange(len(ir_convolve_truncated), 2*len(ir_convolve_truncated)) plt.plot(ir_convolve, label='IR') plt.plot(x_plot, ir_convolve_truncated, label='IR truncated') plt.legend() # Plot spectrum of impulse response H = np.fft.fft(ir_convolve_truncated) len_spec = len(H) half_spectrum = H[:len_spec//2 + 1] amplitude_dB = 20 * np.log10(np.abs(half_spectrum)) f = np.linspace(0, fs/2, len(half_spectrum)) plt.figure(2) plt.semilogx(f, amplitude_dB) # Set FFT length to at least 2x original length to prevent circular wraparound N_fft = 2 ** int(np.ceil(np.log2(2 * x.size))) # next power of 2 print(f"Zero-padding factor: {N_fft/x.size:.1f}x") input_fft = np.fft.fft(x, n=N_fft) output_fft = np.fft.fft(y, n=N_fft) # Calculate transfer function (frequency response) H = output_fft / input_fft # Crop IR: Remove convolution artifacts h = np.fft.ifft(H).real h_truncated = h[ : x.size - 1] plt.figure(3) plt.title('Impulse Response (Frequency Domain Divison)') plt.plot(h, label='IR');plt.plot(h_truncated, label='IR truncated') plt.legend() # Plot the truncated impulse response H_truncated = np.fft.fft(h_truncated) len_spec = len(H_truncated) half_spectrum = H_truncated[:len_spec//2 + 1] f = np.linspace(0, fs/2, len(half_spectrum)) amplitude_dB = 20 * np.log10(np.abs(half_spectrum)) plt.figure(2) plt.semilogx(f, amplitude_dB, linestyle='--',alpha=0.5) plt.xlabel('Frequency (Hz)') plt.ylabel('Amplitude (dB)') plt.title('Frequency Response of H(jw) = Y(jw)/X(jw)') plt.legend(['Convolution with inverse X', 'Spectral division']) plt.grid(True) plt.xlim([10, 10**4]) plt.show()
Improve this question
$\endgroup$
3
  • $\begingroup$ Do you have a complex frequency response or only a real frequency response? $\endgroup$ Commented May 22 at 19:51
  • $\begingroup$ @robert bristow-johnson I receive a complex frequency response, but i'm only interested in the magnitude response. $\endgroup$ Commented May 26 at 11:10
  • $\begingroup$ The resulting impulse response is different if you inverse Fourier Transform the magnitude response than if you inverse Fourier Transform the complete complex frequency response. $\endgroup$ Commented May 26 at 16:17

1 Answer 1

Reset to default
1
$\begingroup$

It's a scaling problem with your time-domain approach.

In the frequency domain,

$$H(f) = \frac{Y(f)}{X(f)}$$

is automatically scale-correct (since the same normalization cancels in numerator / denominator), which is why it shows the “right” amplitude.

However in the time domain, when you convolve your recorded sweep with the inverse filter, you don’t get out a unit‐amplitude impulse unless you explicitly normalize your inverse filter so that

sig.convolve(x, x_inv, mode='full')

produces a Dirac of height $1$.

What’s happening in your case is that

delta = sig.convolve(x, x_inv, mode='full')

has a peak of $\approx 900$ ($\approx 60\tt{dB}$) instead of $1$ so when you do

ir_convolve = sig.convolve(y, x_inv, mode='full')

you’re getting out the true impulse response times that same ~900x factor.

You need to normalize your time-domain impulse response by the "self-convolution" peak of your sweep and inverse filter:

# Time Domain # convolution with inverse filter ir_convolve = sig.convolve(y, x_inv, mode='full') delta = sig.convolve(x,x_inv,mode='full') norm = np.max(np.abs(delta)) # self-convolution peak ir_convolve /= norm
Here is what you get before and after normalization (I simulated a generic system since I don't have access to your data):

Before Normalization

enter image description here

After Normalization

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you very much for your answer and explanation, @Jdip! With the normalization I get almost similar results for both approaches. Evaluated at 100 Hz: 6 dB for the normalized time domain approach and 5 dB for spectral division. Where does that difference come from? Does it come from the different deviating energy of both approaches outside the excitation (10 Hz - 22 kHz), as seen in the spectrum? $\endgroup$ Commented May 26 at 11:06
  • $\begingroup$ You're welcome. The difference you see come mainly from the fact that x_inv $\neq$ ifft(1/fft(x)), which is what the frequency domain approach uses. Farina's method includes the factor 1/k = exp(−t·R/T_sig) which whitens the sweep’s exponentially‐varying energy so that every octave carries equal total energy. 1/fft(x) is just a pointwise magnitude/phase inverse of your particular sampled sweep and does not impose a logarithmic energy equalization. $\endgroup$ Commented May 26 at 23:56

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.