import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# filter design
b, a = signal.butter(2, 0.01, btype='low')

# calculate the FRF
F, H = signal.freqz(b, a, worN=4096, fs=48e3)

# plot the FRF
fig, ax = plt.subplots(2)

ax[0].semilogx(F, 20*np.log10(np.abs(H)))
ax[0].set(ylim=[-50, 10])
ax[0].set(xlabel='Frequency [Hz]')
ax[0].set(ylabel='Magnitude [dB]')
ax[0].grid()

ax[1].semilogx(F, np.angle(H))
ax[1].set(ylim=[-np.pi, 0])
ax[1].set(xlabel='Frequency [Hz]')
ax[1].set(ylabel='Phase [rad]')
ax[1].grid()

plt.show()


np.savetxt('IIR_coefficients_from_python.txt', np.append(b, a[1:]))