#!/usr/bin/python
#jupyter-qtconsole

import numpy as np
from scipy.signal import dlti
from control import *
from scipy import signal
from scipy.signal import freqz, get_window
import matplotlib.pyplot as plt
import sys

#=================================================
#  		CAHIER DES CHARGES
#=================================================
sample_rate=44100.0
Filter_Order = 4
Filter_Type='low' # low / high / bandpass / bandstop
cutoff_1_hz = 0   # 0.0 if low or high pass filter
cutoff_2_hz = 440.0
numtaps = Filter_Order  # Length of the filter (number of coefficients, i.e. the filter order + 1)
#=================================================
#       CALCUL DES COEFFICIENTS
#=================================================
b,a = signal.iirfilter(numtaps, cutoff_2_hz, btype=Filter_Type,
analog=False, ftype='butter', fs=44100,
)
                       
print(b)    
print(a)            
#=================================================
#       TRACES
#=================================================
plt.subplot(211)
f, h = signal.freqz(b,a, fs=sample_rate)
h_dB = 20*np.log10(abs(h))
plt.plot(f,h_dB)

plt.subplot(212)
h_Phase = np.unwrap(np.arctan2(np.imag(h),np.real(h)))
plt.subplot(212)
plt.plot(f,h_Phase)

# MISE EN FORME

plt.subplot(211)
plt.ylabel("Module $H(e^{j\omega})\ en\ dB$")
plt.xlabel("Frequence [Hz]")
plt.xlim([0, 4*cutoff_2_hz])
plt.ylim([-60, 10])
plt.grid()

plt.subplot(212)
plt.ylabel("Phase $H(e^{j\omega})\ en\ rad$")
plt.xlabel("Frequence [Hz]")
plt.xlim([0, 4*cutoff_2_hz])
plt.grid()

plt.show()

#===========================================================