#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Nov 19 2021

@author: Stéphane & François
"""

########################################
#                                      #
#          Filtre passe-bas            #
# sur signaux carré, triangle et sinus #
#                                      #
########################################


import scipy.signal as scsig
import scipy.fft as scfft
import numpy as np
import matplotlib.pyplot as plt

"""
Compléter ci-dessous :
La fonction de transfert complexe
Le facteur de qualité
La fréquence de résonance/coupure
Le nom du filtre
"""
# La fonction filtre donne la valeur de la fonction de transfert complexe
# à la fréquence freq, avec le facteur de qualité Q
# et la frequence de résonance ou de coupure f0
# Attention : 1j est le nombre complexe unité
def filtre(freq,Q,f0):
    return 1/(1+1j*freq/f0) # filtre passe-bas d'ordre 1
# caractéristiques du filtre
Q=10 # facteur de qualité du filtre, si besoin
fc=0.5 # fréquence de résonance ou de coupure du filtre
f0=5 # fréquence des signaux en Hz
Nom_Filtre="Filtre passe-bas d'ordre 1"

# fréquence d'échantillonnage 1000 Hz, durée du signal 1 s
N=1000
t=np.linspace(0,1,N)
offset=2

x1=int(np.log10(fc))  # permet de "centrer" le diagramme de Bode sur fc

# signal carré de fréquence f0
sigcar=scsig.square(2*np.pi*f0*t)+offset

# signal triangle de fréquence f0, symétrique
sigtri=scsig.sawtooth(2*np.pi*f0*t,width=0.5)+offset

# signal sinusoïdal de fréquence f0
sigsin=np.sin(2*np.pi*f0*t)+offset

# fft et récupération des fréquences
transfcar=scfft.fft(sigcar)
transftri=scfft.fft(sigtri)
transfsin=scfft.fft(sigsin)
freq=scfft.fftfreq(len(t),1/N)

# Préparation des fonctions de sortie et cas particulier de la fréquence nulle
specsortcar=np.zeros(len(freq),dtype=complex)
specsortcar[0]=0 # fréquence nulle traitée à part

specsorttri=np.zeros(len(freq),dtype=complex)
specsorttri[0]=0 # fréquence nulle traitée à part

specsortsin=np.zeros(len(freq),dtype=complex)
specsortsin[0]=0 # fréquence nulle traitée à part

# passage du spectre dans le filtre 
for i in range(1,len(freq)-1):
    specsortcar[i]=filtre(freq[i],Q,fc)*transfcar[i]
    specsorttri[i]=filtre(freq[i],Q,fc)*transftri[i]
    specsortsin[i]=filtre(freq[i],Q,fc)*transfsin[i]

# fft inverse pour récupérer les signaux temporels de sortie
sigsortcar=scfft.ifft(specsortcar)+offset
sigsorttri=scfft.ifft(specsorttri)+offset
sigsortsin=scfft.ifft(specsortsin)+offset

#Tracé du diagramme de Bode en gain
def GdB(ff):
    return 20*np.log10(np.abs(filtre(10**ff,Q,fc)))      # Gain en dB 

ff=np.linspace(x1-3,x1+4,1000)   # crée un tableau allant de log(f0)-3 à log(f0)+3 contenant 1000 points
gg=np.ones(1000)                 # crée un tableau de longueur 1000 ne contenant que des 1
gg=GdB(ff)                      # on remplit le tableau gg avec les valeurs de GdB(ff)
# maintenant on trace le diagramme de Bode
plt.grid()    
plt.plot(ff, gg,"b")
plt.title("Diagramme de Bode du filtre $f_c=$ " +str(fc)+ " Hz")
plt.ylabel("Gain en dB")
plt.xlabel("log(fréquence)")
name_Bode='Passe-bas Q '+str(Q)+' fc '+str(fc)+' Bode.png'
plt.savefig(name_Bode,dpi=300,format='png')
plt.show()

# Tracé du signal sinusoïdal
fig, axs = plt.subplots(2,2)
fig.suptitle(Nom_Filtre+" $Q$="+str(Q)+" et $f_0$="+str(f0)+' Hz'+ " sinus")   # titre global
fig.subplots_adjust(hspace=0.5)
fig.subplots_adjust(wspace=0.5)
axs[0,1].set_ylim([-1+offset, 1+offset])
axs[0,0].set(xlabel='temps en s', ylabel='amplitude')
axs[0,0].plot(t,sigsin,'r')
axs[0,1].set(xlabel='temps en s', ylabel='amplitude')
axs[0,1].plot(t,sigsortsin,'b')
axs[1,0].set(xlabel='freq en Hz', ylabel='amplitude')
axs[1,0].plot(np.abs(freq[:100]),np.abs(transfsin[:100])*2/N,'r')
axs[1,0].set(xlabel='freq en Hz', ylabel='amplitude')
axs[1,1].plot(np.abs(freq[:100]), np.abs(specsortsin[:100])*2/N,'b')
name_sin='Passe-bas Q '+str(Q)+' f0 '+str(f0)+' Sinus.png'
plt.savefig(name_sin,dpi=300,format='png')

# Tracé du signal triangulaire
fig, axs = plt.subplots(2,2)
fig.suptitle(Nom_Filtre+" Q="+str(Q)+" et $f_0$="+str(f0)+" Hz"+" triangle")   # titre global
fig.subplots_adjust(hspace=0.5)
fig.subplots_adjust(wspace=0.5)
axs[0,0].set(xlabel='temps en s', ylabel='amplitude')
axs[0,0].plot(t,sigtri,'r')
axs[0,1].set(xlabel='temps en s', ylabel='amplitude')
axs[0,1].plot(t,sigsorttri,'b')
axs[1,0].set(xlabel='freq en Hz', ylabel='amplitude')
axs[1,0].plot(np.abs(freq[:100]),np.abs(transftri[:100])*2/N,'r')
axs[1,0].set(xlabel='freq en Hz', ylabel='amplitude')
axs[1,1].plot(np.abs(freq[:100]), np.abs(specsorttri[:100])*2/N,'b')
name_tri='Passe-bas Q '+str(Q)+' f0 '+str(f0)+' Triangle.png'
plt.savefig(name_tri,dpi=300,format='png')

# Tracé du signal carré
fig, axs = plt.subplots(2,2)
fig.suptitle(Nom_Filtre+" Q="+str(Q)+" et $f_0$="+str(f0)+" Hz"+" carre")   # titre global
fig.subplots_adjust(hspace=0.5)
fig.subplots_adjust(wspace=0.5)
axs[0,0].set(xlabel='temps en s', ylabel='amplitude')
axs[0,0].plot(t,sigcar,'r')
axs[0,1].set(xlabel='temps en s', ylabel='amplitude')
axs[0,1].plot(t,sigsortcar,'b')
axs[1,0].set(xlabel='freq en Hz', ylabel='amplitude')
axs[1,0].plot(np.abs(freq[:100]),np.abs(transfcar[:100])*2/N,'r')
axs[1,0].set(xlabel='freq en Hz', ylabel='amplitude')
axs[1,1].plot(np.abs(freq[:100]), np.abs(specsortcar[:100])*2/N,'b')
name_car='Passe-bas Q '+str(Q)+' f0 '+str(f0)+' Carre.png'
plt.savefig(name_car,dpi=300,format='png')