#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Feb 11 08:52:51 2019

@author: laurentmillot (contacts _DOT_ lmillot _AT_ free _DOT_ fr)
"""




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



h_Array = np.array([
        -0.00042, 0.000334, -0.00058, -0.00029, 0.000359, -0.00089, -0.0001, 0.00033,
      -0.00135, 0.000255, 0.000123, -0.0019, 0.000852, -0.00043, -0.00239, 0.001684, 
        -0.00148, -0.00257, 0.002624, -0.00313, -0.00215, 0.003408, -0.00535, -0.00084,
        0.003646, -0.00793, 0.001561, 0.002858, -0.01048, 0.005134, 0.000525, -0.0124,
        0.009757, -0.00387, -0.0129, 0.015115, -0.01088, -0.01095, 0.020724, -0.02127,
        -0.00498, 0.025989, -0.03688, 0.008325, 0.030301, -0.06444, 0.041262, 0.033129,
        -0.15409, 0.253095, 0.700781, 0.253095, -0.15409, 0.033129, 0.041262, -0.06444,
        0.030301, 0.008325, -0.03688, 0.025989, -0.00498, -0.02127, 0.020724, -0.01095, 
        -0.01088, 0.015115, -0.0129, -0.00387, 0.009757, -0.0124, 0.000525, 0.005134,
        -0.01048, 0.002858, 0.001561, -0.00793, 0.003646, -0.00084, -0.00535, 0.003408,
        -0.00215, -0.00313, 0.002624, -0.00257, -0.00148, 0.001684, -0.00239, -0.00043,
        0.000852, -0.0019, 0.000123, 0.000255, -0.00135, 0.00033, -0.0001, -0.00089,
        0.000359, -0.00029, -0.00058, 0.000334, -0.00042
        ]) 



     
plt.figure(1)
plt.plot(h_Array, 'k-')
plt.grid(True)
plt.xlabel("time (sample)")
plt.title("Impulse response")

plt.savefig("impulse_Response" + ".pdf", format='pdf')


w1, H_Array = signal.freqz(h_Array, 1, 2**18)
G_Array = abs(H_Array)
GdB_Array = 20 * np.log10( abs(H_Array) )
Phi_Array = np.unwrap(np.angle(H_Array))


plt.figure(2)
plt.plot(w1 / w1[-1], GdB_Array, 'k-')
plt.grid(True)
plt.xlabel("reduced frequency (frequency / Nyquist frequency)")
plt.title("Gain (dB)")

plt.savefig("Gain_in_dB" + ".pdf", format='pdf')


plt.figure(3)
plt.plot(w1 / w1[-1], Phi_Array, 'k-')
plt.grid(True)
plt.xlabel("reduced frequency (frequency / Nyquist frequency)")
plt.title("Phase (rad)")

plt.savefig("Phase" + ".pdf", format='pdf')


w2_Array, groupDelay_Array = signal.group_delay((h_Array, 1))


plt.figure(4)
plt.plot(w2_Array/w2_Array[-1], groupDelay_Array, 'k-')
plt.ylim( (-1, 60))
plt.grid(True)
plt.ylabel("Group delay (samples)")
plt.xlabel("reduced frequency (frequency / Nyquist frequency)")
plt.title("Group delay")

plt.savefig("Group_Delay" + ".pdf", format='pdf')