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. As well as plots from the Impulse Response for better understanding for Spectral Division Aproach and for Convolution with Inverse x. 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() Continue reading...
