# -*- coding: utf-8 -*-
# Programme de calcul brut de la TFD d'un signal échantillonné
# Dominique Lefebvre pour TangenteX.com
# 6 avril 2016

# importation des librairies
import matplotlib.pyplot as plt
from scipy import sin, pi, exp, arange, dot, log10

# définition du signal à traiter
Fe = 40.0                   # fréquence d'échantillonnage en Hz
dF = 0.1                    # résolution fréquentielle en Hz

# le nombre d'échantillons est déterminé par la résolution fréquentielle
# hypothèse N = L
N = int(Fe/dF)              # nombre d'échantillons

# la durée du signal à analyser (la période d'échantillonnage Te) dépend du
# nombre d'échantillons à acquérir et de la fréquence d'échantillonnage
Te = 1.0/dF                 # durée du signal à échantillonner

# construction du signal échantillonné s
t = arange(0.0,Te,Te/N)     # vecteur temps discrétisé
a0 = 5.0; f0 = 10.0         # amplitude et fréquence de la sinusoïde
a1 = 2.0; f1 = 5.0
s = a0*sin(2.0*pi*f0*t) + a1*sin(2.0*pi*f1*t)

# construction du vecteur fréquence
f = arange(N)*dF

# définition de la matrice W
n = arange(N)
k = n.reshape((N,1))
W = exp(-2j*pi*k*n/N)

# calcul de la transformée de Fourier discrète par simple produit
# matriciel X = M.s
X = dot(W,s)

# calcul du spectre SP normalisé
SP = (2.0/N)*abs(X)

# pour affichage en dB
SP1 = 20*log10(SP/SP.max())

# affichage du spectre
fig = plt.figure(figsize=(14,8))
plt.grid(True)
plt.title(u"Spectre de fréquence d'un signal composite")
plt.xlabel(u'Fréquence', fontsize = 16)
plt.ylabel('Amplitude (en dB)', fontsize = 16)
plt.plot(f[1:N/2],SP1[1:N/2])