import
numpy as np
import
scipy.signal as signal
import
matplotlib.pyplot as plt
N
=
2
Fs
=
8000
fc
=
3400
Td
=
1
/
Fs
wd
=
2
*
np.pi
*
fc
print
(wd)
wc
=
(
2
/
Td)
*
np.tan(wd
*
Td
/
2
)
print
(
'Order of the filter='
, N)
print
(
'Cut-off frequency (in rad/s)='
, wc)
b, a
=
signal.butter(N, wc,
'low'
, analog
=
'True'
)
z, p
=
signal.bilinear(b, a, fs
=
Fs)
print
(
'Numerator Coefficients:'
, z)
print
(
'Denominator Coefficients:'
, p)
wz, hz
=
signal.freqz(z, p,
512
)
fig
=
plt.figure(figsize
=
(
10
,
8
))
Mag
=
20
*
np.log10(
abs
(hz))
Freq
=
wz
*
Fs
/
(
2
*
np.pi)
sub1
=
plt.subplot(
2
,
1
,
1
)
sub1.plot(Freq, Mag,
'r'
, linewidth
=
2
)
sub1.axis([
1
, Fs
/
2
,
-
60
,
5
])
sub1.set_title(
'Magnitude Response'
, fontsize
=
15
)
sub1.set_xlabel(
'Frequency [Hz]'
, fontsize
=
15
)
sub1.set_ylabel(
'Magnitude [dB]'
, fontsize
=
15
)
sub1.grid()
sub2
=
plt.subplot(
2
,
1
,
2
)
Phase
=
np.unwrap(np.angle(hz))
*
180
/
np.pi
sub2.plot(Freq, Phase,
'g'
, linewidth
=
2
)
sub2.set_ylabel(
'Phase (degree)'
, fontsize
=
15
)
sub2.set_xlabel(r
'Frequency (Hz)'
, fontsize
=
15
)
sub2.set_title(r
'Phase response'
, fontsize
=
15
)
sub2.grid()
plt.subplots_adjust(hspace
=
0.5
)
fig.tight_layout()
plt.show()