1
$\begingroup$

I am struggling to simulate the frequency response of a simple op-amp low-pass transfer function in Python. The results I get are not accurate. The transfer function is $H(s)=\frac{1}{1 + s\tau}$, where $\tau=\frac{1}{2\pi f_{cutoff}}$. So I would expect a frequency response with a gain of 0dB until the cutoff frequency, then -20dB/dec, classic. But as shown below, I have gain of $\approx -7.5dB$, and also the cutoff frequency seems to be off. What am I doing wrong?

BTW, in case you are wondering why I filter a sine-wave going through the PA - just to check that the filtered sine ends up on a single bin to have a "clean" result, I don't want any transients, thus the cyclic convolution (take-away from this question Unexpected frequency components after applying bandpass filter in python).

import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.ndimage import convolve1d
from scipy.fftpack import fft
import control
from scipy import signal

# Parameters
f_in = 1e6
M = 64
f_samp = f_in * M
f_pa_cutoff = 10e6
N = int((f_samp / f_in) * 128)

# Signals and Axis
t = np.linspace(0, float((1 / f_samp) * (N - 1)), num=int(N), dtype=float)
x = np.sin(2 * math.pi * f_in * t)

freq_axis = np.fft.fftfreq(N, 1/f_samp)
freq_axis = freq_axis[0:N//2]

# PA transfer function
imp = np.concatenate(([1], np.zeros(len(t) - 1)))

pa_tau = 1 / (2 * math.pi * f_pa_cutoff)
pa_gain = 1
pa_tf_vec = ([0, pa_gain], [pa_tau, 1])

_, sig_opamp_imp_resp, _ = signal.lsim(pa_tf_vec, U=imp, T=t)
y = convolve1d(x, np.array(sig_opamp_imp_resp), mode='wrap')

# Plot FFT
fft_x = 20*np.log10(np.divide(np.abs(fft(x)), N/2))[0:N//2]
fft_y = 20*np.log10(np.divide(np.abs(fft(y)), N/2))[0:N//2]
fft_opamp_filt = 20*np.log10(np.abs(fft(sig_opamp_imp_resp)))[0:N//2]

plt.figure()
plt.semilogx(freq_axis, fft_opamp_filt, label = 'PA TF', color='green', linestyle='dashed')
plt.scatter(freq_axis, fft_x, label = 'PA In', marker = 'x')
plt.scatter(freq_axis, fft_y, label = 'PA Out', marker = 'x')
plt.grid(True)
plt.xlabel('Frequency [Hz]')
plt.ylabel('[dB]')
plt.ylim(-100, 5)
plt.show()

enter image description here

$\endgroup$

0

Your Answer

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