import math
import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fft,rfft 

#Considérons comme exemple un signal périodique (de période 1) :
    
def signal(t):
    return 0.1+np.cos(2*np.pi*t)+0.5*np.cos(2*2*np.pi*t)+0.25*np.sin(3*2*np.pi*t)+0.01*np.sin(4*2*np.pi*t)
    
# Ce signal a normalement des pics à 0 (amplitude 0.1), 1 (amplitude 1), 2 (amplitude 0.5), 3 (amplitude 0.25) et 4 (amplitude 0.01).    

#Le signal est échantillonné sur une durée T qui doit être beaucoup plus grande que sa période. On fixe cette durée à une valeur qui ne coïncide pas avec une période du signal, comme c'est le cas dans la numérisation des signaux physiques. Dans le cas présent, le choix de l'origine du temps est sans importance ; on choisit l'intervalle [0,T]. La fréquence d'échantillonnage doit être supérieure au double de la plus grande fréquence du signal, ici 4 => on doit donc choisir fe > 8

    
T=100.345821
fe = 10
t=np.linspace(0,T,int(T*fe))   # T*fe=N points régulièrement répartis
echantillons = signal(t)

#On calcule la transformée de Fourier discrète sur [0,fe] puis sa norme (spectre) :
tfd = fft(echantillons)
N=len(echantillons)
spectre = np.absolute(tfd)*2/N

#La fonction rfft calcule la même chose mais un intervalle réduit [0,fe/2]:
rtfd=rfft(echantillons)
spectrebis=np.absolute(rtfd)*2/N
            
#L'échelle des fréquences est construite sachant que l'espacement fréquentiel de deux points de la TFD est l'inverse de la durée T : 

freq=np.linspace(0,fe,N)

#L'échelle des fréquences pour rfft est construite sachant que l'espacement fréquentiel de deux points de la TFD est l'inverse de la durée T mais que l'intervalle est réduit à [0,fe/2] :

freqbis=np.linspace(0,fe/2,(N//2+1))
         
#On trace le spectre pour des fréquences allant de zéro à fe/2 :

plt.figure(figsize=(10,4))
plt.plot(freq,spectre,'r')
plt.xlabel('f')
plt.ylabel('A')
plt.axis([-0.1,fe/2,0,spectre.max()])
plt.title("Calcul du spectre avec fft")
plt.grid()
plt.show()

#On peut aussi tracer le spectre avec rfft sur l'intervalle deja reduit [0,fe/2] :

plt.figure(figsize=(10,4))
plt.plot(freqbis,spectrebis,'r')
plt.xlabel('f')
plt.ylabel('A')
plt.title("calcul du spectre avec rfft")
plt.grid()

plt.show()                        

#L'harmonique de rang 4 est quasiment invisible. Une représentation en décibel permet de voir ces composantes faibles :

spectre_db = 20*np.log10(spectre/spectre.max())
plt.figure(figsize=(10,4))
plt.plot(freq,spectre_db,'r')
plt.xlabel('f')
plt.ylabel('AdB')
plt.axis([-0.1,fe/2,spectre_db.min(),spectre_db.max()])
plt.grid()
plt.title("Calcul avec fft en décibels")
plt.show()