\documentclass{TDPhysique}
\partie{Informatique pour tous}
\chapitre{Transformée de Fourier discrète}
\annee{2020-2021}


\begin{document}

\begin{TD}
\begin{center}
\begin{LARGE}
\textcolor{DodgerBlue!50!black}{Autour de la transformée de Fourier discrète}
\end{LARGE}
\end{center}
\end{TD}


\section*{Introduction}

Des renseignements très utiles sur un signal donné\footnote{D'un point de vue mathématique, le signal est une fonction réelle ou complexe de la variable réelle.} peuvent être obtenus en déterminant son spectre fréquentiel. 

La découverte
de l'intérêt de la décomposition spectrale d'un signal est due à Joseph Fourier qui a établi que tout signal périodique peut se décomposer en une somme finie de signaux sinusoïdaux de fréquences et d'amplitudes constantes. L'invention de Fourier, en 1811, à 43 ans, tire son origine de la résolution de l'équation de la chaleur à l'aide d'une fonction périodique, somme de fonctions trigonométriques. 

C'est l'ensemble fini des amplitudes et des fréquences de ces fonctions trigonométriques que l'on nomme le spectre du signal. 

\begin{flushleft}

\includegraphics[width=18cm]{Fig_Intro_Spectre}

\end{flushleft}


L'outil mathématique pour obtenir la décomposition d'un signal périodique et par conséquent le spectre correspondant est le développement en séries de Fourier (DSF). 

Dans le cas d'un signal qui n'est pas périodique, une analyse comparable du signal peut être effectuée pour mettre en évidence les composantes
fréquentielles principales d'un signal. On a recours dans ce cas à un outil mathématique appelé transformée de Fourier (TF). 


Enfin, l'analyse spectrale d'un signal numérique, donc discret, fait intervenir quant à elle un outil mathématique appelé transformée de Fourier discrète (TFD).





\begin{center}
	\begin{tabular}{l|c|c|c|c|c|c|}
	 & (1) & & (2) &  & (3) &  \\
\rule[-2ex]{0pt}{5.5ex} Signal d'entrée & $\longrightarrow$ & Spectre & $\longrightarrow$ & Spectre filtré & $\longrightarrow$ & Signal de sortie \\
\hline
\hline
\rule[-2ex]{0pt}{5.5ex} périodique  & DSF & & filtrage &  & DSF$^{-1}$ &  \\
\hline
\rule[-2ex]{0pt}{5.5ex} non-périodique  & TF & & filtrage &  & TF$^{-1}$ &  \\
\hline
\rule[-2ex]{0pt}{5.5ex} discret  & TFD & & filtrage &  & TFD$^{-1}$ &  \\
\hline
	\end{tabular} 
\end{center}

DSF : Développement en séries de Fourier

TF : Transformée de Fourier

TFD : Transformée de Fourier discrète\\

Si le signal d'entrée est une fonction du temps, on parle alors d'un spectre de fréquences temporelles. C'est le cas de tous les signaux temporels de la physique.

Si le signal d'entrée est une fonction de l'espace, on parle alors d'un spectre de fréquences spatiales. C'est le cas dans les images et leur traitement.

Dans les deux cas, l'opération de filtrage a lieu sur le spectre en fréquence. Le tableau ci-dessus en résume le principe.

En dehors de l'analyse spectrale d'un signal, les outils mathématiques que sont les séries de Fourier et la transformée de Fourier sont également utiles pour d'autres usages comme la résolution d'équations aux dérivées partielles, le calcul de la valeur de certaines intégrales, le calcul de la valeur de la somme de certaines séries, etc. 


 
\section{Présentation de la TFD}

Nous allons travailler sur des signaux numériques qui sont donc par nature discrétisés. Obtenir le spectre, ici en fréquence, d'un signal numérisé en fonction du temps, consiste à chercher la transformée de Fourier discrète (TFD) de la série de valeurs prises par le signal d'entrée. 

\schema{15cm}{Fig_Echantillonnage}


	\subsection*{Définition}
La transformée de Fourier discrète -- TFD -- réalise la correspondance entre deux suites de $N$ termes.
\begin{center}
$\{s_k$, $k\in[0,N-1]\}$ \qquad $\longrightarrow$ \qquad $\{S_m$, $m\in[0,N-1]\}$
avec 
$$S_m=\sum_{k=0}^{N-1}s_k e^{-j\frac{2\pi}{N}km} $$
\end{center}

En posant $w_N=e^{-j\frac{2\pi}{N}}$, racine N-ième de l'unité, la TFD s'écrit 

\begin{equation}
S_m=\sum_{k=0}^{N-1}s_k w_N^{km} 
\label{Eq1}
\end{equation}

$s_k$ est réel alors que $S_m$ est complexe. 
Le module de $S_m$, $|S_m|$ en fonction de $m$ donne le spectre en amplitude du signal $\{s_k\}$ et $\arg{(S_m)}$ sa phase. 


	\subsection*{Mise en forme matricielle de la TFD}
	
La somme précédente a une forme matricielle. Par exemple, pour $N=4$

\begin{eqnarray*}
 	S_0 = s_0(w_4)^{0.0}+s_1(w_4)^{1.0}+s_2(w_4)^{2.0}+s_3(w_4)^{3.0}\\
 	S_1 = s_0(w_4)^{0.1}+s_1(w_4)^{1.1}+s_2(w_4)^{2.1}+s_3(w_4)^{3.1}\\
 	S_2 = s_0(w_4)^{0.2}+s_1(w_4)^{1.2}+s_2(w_4)^{2.2}+s_3(w_4)^{3.2}\\
 	S_3 = s_0(w_4)^{0.3}+s_1(w_4)^{1.3}+s_2(w_4)^{2.3}+s_3(w_4)^{3.3}\\
\end{eqnarray*}
 	

soit 


\[\begin{pmatrix}
   S_0 \\
   S_1 \\
   S_2 \\
   S_3 \\
\end{pmatrix}
=
\begin{pmatrix}
   1 & 1 & 1 & 1 \\
   1 & w & w^2 & w^3\\
   1 & w^2 & w^4 & w^6\\
   1 & w^3 & w^6 & w^9\\
\end{pmatrix}
\begin{pmatrix}
   s_0 \\
   s_1 \\
   s_2 \\
   s_3 \\
\end{pmatrix}\]
	
	
	
soit 

\begin{equation}
\mathcal{S}=W_N.\textit{s} \label{Eq2}
\end{equation}

 où  $\mathcal{S}$ est le vecteur  complexe de la transformée de Fourier discrète  du vecteur signal $\textit{s}$ 

	\quest	Montrer que le vecteur signal $\textit{s}=\{2,3,-1,1\}$ à comme TFD le vecteur spectre $\textit{S}=\{5+0j,3-2j,-3+0j,3+2j\}$.


\[\begin{pmatrix}
   5 \\
   3-2j \\
   -3 \\
   3+2j \\
\end{pmatrix}
=
\begin{pmatrix}
   1 & 1 & 1 & 1 \\
   1 & -j & -1 & j\\
   1 & -1 & 1 & -1\\
   1 & j & -1 & -j\\
\end{pmatrix}
\begin{pmatrix}
   2 \\
   3 \\
   -1 \\
   1 \\
\end{pmatrix}\]

	\quest Le spectre consiste à tracer le module des composantes de $\mathcal{S}$ en fonction de $m$ : $|S_m|(m)$ 

Tracer le spectre du vecteur signal précédent. 

\begin{center}
\schema{8cm}{Fig_ExempleSpectre}
\end{center}


		\quest Montrer que $\frac{S_0}{N}$ est la valeur moyenne du signal $\{s_k\}$.\\

\hrulefill		
		
$S_0$ correspond à l'amplitude de la fréquence nulle et est nécessairement réelle. Cette amplitude \textit{serait} la valeur moyenne du signal ${s_k}$ si on la divisait par $N$. En effet, on a 
$$S_0=\sum_{k=0}^{N-1}s_k \qquad \rightarrow \qquad \frac{S_0}{N}=\frac{1}{N}\sum_{k=0}^{N-1}s_k = \langle s_k \rangle$$
		
\hrulefill		
		
		\quest Le spectre en fréquence est donc le module des coefficients $S_m$ complexes obtenus par la TFD. Montrer que pour un spectre de $N$ valeurs, $S_{N-m}=\overline{S_m}$ et que par conséquence le spectre tracé sur la totalité des points est dédoublé : le spectre attendu et son spectre miroir. \\
		
\hrulefill
		
Pour $1\leq m \leq N-1$, on montre que $S_{N-m}=\overline{S_m}$ ce qui réduit la partie utile du spectre en fréquence, puisque l'on trace le module de $S_m$, $|S_m|(m)$.

En effet, 

\begin{align*}
S_{N-m}	&=\sum_{k=0}^{N-1}s_k e^{ -j\frac{2\pi}{N}k(N-m)}\\  	
		&=\sum_{k=0}^{N-1}s_k e^{ -j\frac{2\pi}{N}kN}.e^{j\frac{2\pi}{N}km}\\
		&=\sum_{k=0}^{N-1}s_k e^{j\frac{2\pi}{N}km}\\
		&=\overline{S_m}
\end{align*}


Finalement, pour tracer le spectre en fréquence, soit le module de $S_m$ en fonction de $m$, nous n'avons besoin que des valeurs $S_m$, pour $1\leq m\leq N/2$.

\hrulefill
	
		
\section{Algorithme naïf d'obtention du vecteur spectre $\mathcal{S}$}
On charge les bibliothèques suivantes : 

\inputpython{TFD_Latex_Algorithmes.py}{1}{2}

La bibliothèque \pyth{numpy} permet de travailler simplement avec des matrices à valeurs complexes. Une annexe donne quelques fonctions utiles.
	
	\subsection*{Fonction $w_N$}	

	\quest Coder la fonction \pyth{w(N)} définie pour l'équation (\ref{Eq1}). Cette fonction prend en argument $N$, le nombre de points du signal.
	
\hrulefill	
	
On construit dans un premier temps la fonction \pyth{w(N)} définie pour l'équation (\ref{Eq1}).



\inputpython{TFD_Latex_Algorithmes.py}{4}{5}

\hrulefill

	\subsection*{TFD à partir de la somme}

	\quest Écrire une fonction \pyth{TFD1(s)}, qui prend en argument un vecteur  signal \pyth{s} et qui renvoie un vecteur  spectre \pyth{S}, à partir de la définition donnée par l'équation (\ref{Eq1}).

\hrulefill

\inputpython{TFD_Latex_Algorithmes.py}{8}{14}

\hrulefill


	\quest Vérifier que la fonction précédente donne $\textit{S}=\{5,3-2j,-3,3+2j\}$  à partir du signal $\textit{s}=\{2,3,-1,1\}$.
	
	On prendra donc ici pour \pyth{s}, \pyth{s = np.array([2,3,-1,1],dtype = complex)}

On pourra arrondir les flottants avec la fonction \pyth{numpy} : \pyth{np.round()}

\hrulefill

\begin{python}
>>> s = np.array([2,3,-1,1],dtype = complex)
>>> s 
array([ 2.+0.j,  3.+0.j, -1.+0.j,  1.+0.j])
>>> TFD1(s)
array([ 5. +0.00000000e+00j,  3. -2.00000000e+00j, -3. -9.79717439e-16j,
        3. +2.00000000e+00j])
>>> np.round(TFD(s)
array([ 5.+0.j,  3.-2.j, -3.-0.j,  3.+2.j])
>>> 
\end{python}

\hrulefill


	\subsection*{TFD à partir de la matrice $W_N$}

		\quest Écrire une fonction \pyth{M_TFD(N)} qui prend en argument le nombre d'éléments du signal -- dimension du vecteur signal -- et qui génère la matrice $W_N$ adéquate définie par (\ref{Eq2}).

\hrulefill

\inputpython{TFD_FFT_Algorithmes.py}{17}{23}

\hrulefill

		\quest Écrire une fonction \pyth{TFD(s)}, qui prend en argument un vecteur  signal et qui renvoie la TFD à partir de son écriture matricielle.

\hrulefill

\inputpython{TFD_Latex_Algorithmes.py}{27}{31}

\hrulefill

Vérifier à nouveau que la fonction créée donne $\textit{S}=\{5,3-2j,-3,3+2j\}$  à partir du signal $\textit{s}=\{2,3,-1,1\}$.

	\subsection*{Complexité}

	\quest Évaluer la complexité de l'algorithme naïf. 

\hrulefill
	
Les boucles imbriquées de la fonction \pyth{MTFD(s)} indique que $$\mathcal{C}(N)=\mathcal{O}(N^2)$$

\hrulefill

	\subsection*{TFD inverse : obtenir le vecteur signal $s$}

Une fois le signal filtré, et par conséquent le vecteur spectre $\mathcal{S}$ modifié, nous  construisons le signal de sortie à partie de la TDF inverse telle que : 

$$s_k = \frac{1}{N}\sum_{m=0}^{N-1}S_m w_N^{-mk}$$

soit, sous forme matricielle, 

\begin{equation}
\textit{s} = \frac{1}{N}\overline{W_N}\mathcal{S}
\label{Eq3}
\end{equation}

où $\overline{W_N}$, la matrice conjuguée de $W_N$ est telle que $$\overline{W_N}.W_N=N.I_N$$ 

$I_N$ est la matrice identité d'ordre $N$ et $\overline{W_N}=(\overline{w_N^{km}})_{0\leq m,k \leq N-1}$. 

		\quest Écrire une fonction \pyth{M_TFDI(N)} qui prend en argument le nombre \pyth{N}  d'éléments du vecteur spectre et qui génère la matrice $\overline{W_N}$, matrice conjuguée de $W_N$.

\hrulefill

\inputpython{TFD_Latex_Algorithmes.py}{37}{39}

\hrulefill

		\quest Écrire une fonction \pyth{TFDI(S)} qui prend en argument un vecteur spectre $\mathcal{S}$ et que renvoie le vecteur signal \textit{s}.

Tester cette fonction à partir de l'exemple précédent. 

\hrulefill

\inputpython{TFD_Latex_Algorithmes.py}{42}{46}
	
\hrulefill

	\subsection*{Filtrage}
Le principe du filtrage est relativement simple. 

On peut modifier\footnote{... atténuer, amplifier, déphaser etc.} les composantes (complexes) du vecteur spectre $\mathcal{S}$ associées aux fréquences du signal. 


Pour que ce filtrage soit parlant pour un physicien, nous allons travailler sur un signal physique dont le fondamental est de fréquence $f_s$ échantillonné à la fréquence $f_e$ sur une durée totale d'acquisition $t_A$.

		\subsection*{Présentation d'un signal physique}\label{signalPhysique}
		
	\begin{center}
	\begin{tabular}{|c|c|}
	\hline 
	\rule[-1.5ex]{0pt}{4.5ex} \textsc{Signal fonction du temps} & \textsc{Spectre fonction de la fréquence} \\ 
	\hline 
	\hline
	\multicolumn{2}{|c|}{Avec Fourier, tout est inversé} \\ 
	 
	\hline 
	\hline
	\rule[-1ex]{0pt}{2.5ex} \qquad Durée totale d'acquisition du signal \qquad& Fréquence d'échantillonnage du spectre \\ 

	\rule[-3ex]{0pt}{6.5ex} $t_A$ & $\displaystyle f_A= \frac{1}{t_A}$ \\ 
\hline
\hline
	\multicolumn{2}{|c|}{\rule[-1ex]{0pt}{3.5ex} \'Echelle temporelle la plus grande $\leftrightarrow$ \'Echelle fréquentielle la plus petite	} \\ 
\hline 
\hline 
	Période du signal & Fréquence du signal\\
	\rule[-3ex]{0pt}{6.5ex} $ \displaystyle T_s=\frac{1}{f_s}$ & $f_s$ \\ 
	\hline 
	Période d'échantillonnage du signal & \qquad Fréquence d'échantillonnage du signal \qquad \\
	\rule[-3ex]{0pt}{6.5ex} $T_e = \displaystyle\frac{T_s}{N_s}$ & $f_e=N_s.f_s$\\ 

\hline
\hline
	\multicolumn{2}{|c|}{\rule[-1ex]{0pt}{3.5ex} \'Echelle temporelle la plus petite $\leftrightarrow$ \'Echelle fréquentielle la plus grande	} \\ 
\hline 
\hline 
	
	\multicolumn{2}{|c|}{$N_s$ est le nombre de points échantillonnés par période} \\ 	\multicolumn{2}{|c|}{$N_s>2$ est le critère de Shannon} \\ 	
	\hline 
	\rule[-1ex]{0pt}{2.5ex} $t_A=N.T_e$ & $f_e=N.f_A$ \\ 
	\hline
	\hline
	
	\multicolumn{2}{|c|}{$N$ est le nombre total de points échantillonnés du signal} \\ 	
	\multicolumn{2}{|c|}{$N=2^n$ doit être une puissance de 2 : $1024 = 2^{10}$} \\ 	
	\hline 
	\end{tabular} 
	\end{center}

\begin{center}

\schema{15cm}{Fig_ExempleSignal}
\end{center}

\begin{center}
\schema{16cm}{Fig_ExempleSpectreCalcule}
\end{center}


	\subsection*{Construction et tracé du signal \pyth{signal} }

	\quest On souhaite construire un signal \pyth{s} de $N=2^{10}=1024$ échantillons, somme de trois sinusoïdes de même amplitude, mais de fréquences $f_s=100$ Hz, $2.f_s$ et $5.f_s$, et qui présente 100 échantillons par période $T_s=\frac{1}{f_s}$ du fondamental du signal.   \\
	
Définir la fonction \pyth{fn(x)} qui prend en entrée $x$ la fréquence du fondamental et renvoie le signal demandé.

\hrulefill


\begin{python}
def fn(x):
    return  np.sin(x) + np.sin(2*x) + np.sin(5*x)
\end{python}
	
\hrulefill	
	
	\quest 	Compléter les lignes du code suivant : 

\begin{python}
N = 2**10 # le nombre d echantillons est une puissance de 2
fs = A COMPLETER #(en Hz) frequence du signal
Ts = 1/fs # (en s) du signal
Ns = A COMPLETER # nombre de points par periode du signal
fe = A COMPLETER fe est fonction de Ns et fs # frequence d echantillonage. Ns>2 Shannon largement verifiee
Te= A COMPLETER # periode d echantillonage 
tA = A COMPLETER # tA duree d acquisition/totale du signal 
\end{python}

\hrulefill

\begin{python}
N = 2**10 # nombre de points lors de l acquisition puissance de 2
fs = 100 #(en Hz) du signal
Ts = 1/fs # (en s) du signal
Ns = 100  # nombre de points pour echantillonage
fe = Ns*fs # frequence de l echantillonage Ns>2 Shannon largement verifiee
Te=1/fe # periode d echantillonage 
#duree d acquisition tA
tA = N*Te #
\end{python}

\hrulefill



	\quest A l'aide de la fonction \pyth{np.arange(debut,fin,pas)} présentée en annexe, générer une liste \pyth{numpy} des temps \pyth{t} sur le durée totale du signal avec le pas adéquate. 
	
\hrulefill	
	
\begin{python}	
#on peut generer la liste des temps
t = np.arange(0,tA,Te) #N points car tA = N*Te
\end{python}

\hrulefill

	\quest Générer le signal \pyth{signal} à l'aide la fonction \pyth{fn(x)} puis tracer \pyth{signal}. 
	
	Pour le tracé, compléter le code suivant : 
	
\begin{python}
plt.plot(A COMPLETER)
plt.xlabel('temps [s]')
plt.ylabel('signal numerique lisse')
plt.show()
\end{python}  

\hrulefill

\begin{python}
signal = fn(2*np.pi*fs*t)
plt.subplot(4,1,1)
plt.plot(t,signal)
plt.xlabel('temps [s]')
plt.ylabel('signal numerique lisse')
plt.show()
\end{python}

\hrulefill

	\subsection*{Calcul et tracé du spectre \pyth{spectre} de \pyth{signal}}



	\quest Calculer le spectre en fréquence du signal \pyth{signal}.

\hrulefill

\begin{python}
#Calcul du spectre
spectre = TFD(signal)
\end{python}

\hrulefill

	\quest Pour tracer le module des composantes complexes du spectre en fonction de \pyth{f}, il faut générer la liste \pyth{numpy} de fréquences. A partir de la fonction \pyth{numpy}, \pyth{np.arange(debut,fin,pas)}, générer la liste des fréquences \pyth{f} associées à ce signal. 
	
	Tracer le module des composantes complexes du spectre en fonction de \pyth{f}. On indique que le module d'un complexe \pyth{A=0.1+2.j} et donné par \pyth{abs(A)}.
	
\hrulefill	
	
\begin{python}
#Pour le spectre en frequence
fA=1/tA # On genere les valeurs de frequences
f = np.arange(0,fe,fA)
#Trace du spectre sur N//2
plt.subplot(4,1,2)
plt.plot(f[:N//2],abs(spectre[:N//2]))
plt.xlabel('frequence [Hz]')
plt.ylabel('spectre calcule')
plt.show()
\end{python}



Comme nous l'avons montré, le tracé du spectre se réduit à un tracé sur la moitié des valeurs calculées, soit ici de $0$ à $\frac{f_e}{2}$.

\begin{center}
\schema{16cm}{Fig_ExempleSpectreCalcule}
\end{center}

\hrulefill

	\subsection*{Filtrage}
	
		\quest Réaliser un filtre passe-bas numérique consiste, dans une conception extrêmement simple, à mettre à 0 l'ensemble des valeurs du vecteur spectre supérieures à $1.5f_s$ si on a choisi une fréquence de coupure telle que $f_c=1.5f_s$. 
		
		Écrire une fonction \pyth{filtrePBas(spectre)} qui prend en argument un vecteur spectre complexe et qui renvoie le spectre filtré à $1.5f_s$. L'appliquer au spectre précédent et tracer le module des composantes de ce spectre filtré.

\hrulefill

\inputpython{TFD_Latex_Algorithmes.py}{120}{125}
	
\begin{python}
spectre_filtre = filtrePBas(spectre)
\end{python}

\begin{python}
#Trace du spectre filtre sur N//2
plt.subplot(4,1,3)
plt.plot(f[:N//2],abs(spectre_filtre[:N//2]))
plt.xlabel('frequence [Hz]')
plt.ylabel('spectre calcule')
plt.show()
\end{python}

		\quest Une fois le vecteur spectre modifié, on lui applique la transformée de Fourier discrète inverse définit précédemment \pyth{TDFI(spectre)}. Pour terminer, tracer la partie réelle du signal filtré.
		
\hrulefill

\begin{python}
signal_filtre = TFDI(spectre_filtre) 
\end{python}


\hrulefill

\begin{python}
plt.subplot(4,1,4)
plt.plot(t,signal_filtre).real)
plt.show()
\end{python}

\begin{center}
\includegraphics[width=15cm]{Fig_Filtrage_passe_bas.pdf}
\caption{Filtrage passe-bas}
\end{center}

		\quest Réaliser un filtre passe-bande numérique consiste, dans une conception extrêmement simple à mettre à 0 l'ensemble des valeurs du vecteur spectre inférieures et supérieures à deux valeurs choisies. 


\hrulefill

\inputpython{TFD_Latex_Algorithmes.py}{129}{134}
	
\begin{python}
spectre_filtre = filtrePBande(spectre)
\end{python}

Une fois le vecteur spectre modifié, on lui applique la transformée de Fourier discrète inverse et on extrait la partie réel du signal calculé.

\begin{python}
signal_filtre = (TFDI(spectre_filtre)).real 
\end{python}


\hrulefill

\begin{center}
\includegraphics[width=15cm]{Fig_Filtrage_passe_bande.pdf}
\caption{Filtrage passe-bande autour de $2.f_s$}
\end{center}



\newpage{}

\section{Algorithme de TFD rapide : FFT de James Cooley et John Tukey (1965)}

	\subsection*{Complexité de l'algorithme naïf de TFD : un inconvénient rédhibitoire}

		\quest Quelle est la complexité d'une TFD ? 

\hrulefill

On multiplie une matrice $(N\times N)$ par une matrice colonne de taille $N$. On a donc une complexité en $\mathcal{O}(N^2)$. 

Pour un vecteur de taille 1000, on a donc $10^6$ opérations. Puis pour un vecteur de taille 3000, le temps de calcul est par $3^2=9$

Pour des petites valeurs de $N$, cela ne pose pas de problème. En revanche, le temps de calcul devient prohibitif pour des applications en temps réel de la TFD, par exemple, pour les sons audibles par une oreille humaine qui ont des fréquences s'étendant jusqu'à environ 20000 Hz. 

Pour analyser correctement un son, il faudra donc l'échantillonner à 40000 Hz, nous avons environ 2 milliards d'opérations. Si c'est faisable en quelques secondes, cela interdit le traitement du son nécessaire en temps réel, comme la simple écoute d'un son numérisé. 

Finalement, l'algorithme naïf est trop lent. 

Dans la suite nous étudions l'algorithme de TFD rapide de Cooley et Tukey beaucoup plus rapide. 


\hrulefill

 	\subsection*{Histoire de la naissance de la TFD rapide, FFT : Fast Fourier Transform}
	

	Cet algorithme\footnote{L'histoire de la naissance de cet algorithme est tirée de l'encyclopédie en ligne  Wikipedia.}, y compris sa version récursive, a été inventé vers 1805 par Carl Friedrich Gauss, qui l'a utilisé pour interpoler les trajectoires des astéroïdes Pallas et Junon. Son travail est resté inaperçu du fait d'avoir été publié uniquement à titre posthume et en néo-latin. Il semble que Gauss n'ai pas analysé le temps de calcul asymptotique. 
	
	Diverses formes limitées ont également été redécouvertes à plusieurs reprises au cours du XIX\up{e} et au début du XX\up{e} siècle. 
	
	Les FFT sont devenues populaires après que James Cooley d'IBM et John Tukey de l'université de Princeton aient publiés un article en 1965 réinventant l'algorithme et décrivant comment l'exécuter facilement sur un ordinateur. 
	
	Tukey aurait eu l'idée lors d'une réunion du comité consultatif scientifique du président Kennedy sur les moyens de détecter les essais d'armes nucléaires en Union soviétique en utilisant les données des sismomètres situés à l'extérieur du pays.
	
	La difficulté venait du fait que l'analyse de cette grande quantité de données (nombre de capteurs et durée d'acquisition) nécessitait des algorithmes rapides pour calculer la TFD. 
	
	Cette tâche était essentielle pour la ratification du projet d'interdiction des essais nucléaires afin que toute violation puisse être détectée sans qu'il soit nécessaire de visiter les installations soviétiques. 
	
	Un autre participant à cette réunion, Richard Garwin d'IBM, qui avait saisi le puissance de cette méthode, a mis Tukey en contact avec Cooley en s'assurant toutefois que Cooley ne connaissait pas le but initial.  Il a alors été dit à Cooley que ce travail était nécessaire pour déterminer les périodicités des orientations de spin dans un cristal 3-D d' hélium-3. 
	
	La publication de l'article de Cooley et Tukey, suivi du développement simultané de convertisseurs analogique-numérique capables d'échantillonner à des fréquences allant jusqu'à 300 kHz a rendu l'algorithme de la FFT très opératoire. 

	\subsection*{Principe de l'algorithme récursif de la FFT : \og diviser pour régner \fg{}  }

		
L'algorithme de la FFT divise, à chaque étage récursif, une TFD de taille $N$ en deux TFD entrelacées de taille $N/2$. 

Cette idée provient d'une écriture récursive de la TFD pour laquelle on sépare la somme initiale en deux sommes : une somme sur les indices pairs $k=2p$, et une somme sur les indices impairs $k=2p+1$. Cela impose à $N$ d'être une puissance de 2 : $N=2^q$. \\

		\quest Monter que $$S_m=\underbrace{\sum_{p=0}^{N/2-1}s_{2p} e^{-j\frac{2\pi}{N/2}pm}}_{\text{TDF des index pairs } s_{2p}}+e^{-j\frac{2\pi}{N}m}\underbrace{\sum_{p=0}^{N/2-1}s_{2p+1} e^{-j\frac{2\pi}{N/2}pm}}_{\text{TDF des index impairs } s_{2p+1}}  $$

soit $$S_m=E_m+e^{-j\frac{2\pi}{N}m}O_m  $$

si on note la TFD des signaux d'index pair $E_m$ (Even-index pour index pairs) et la TFD des entrées d'index impair $O_m$ (Odd-index pour index impairs) telles que 


 $$E_m=\underbrace{\sum_{p=0}^{P-1}e_{p} e^{-j\frac{2\pi}{P}pm}}_{\text{TDF des index pairs } s_{2p}=e_p}$$
 
 et 
 
 $$O_m=\underbrace{\sum_{p=0}^{P-1}o_{p} e^{-j\frac{2\pi}{P}pm}}_{\text{TDF des index impairs } s_{2p+1}=o_p}  $$ 


\hrulefill

Rappelons la définition de la transformée de Fourier discrète (TFD), 

$$S_m=\sum_{k=0}^{N-1}s_k e^{-j\frac{2\pi}{N}km}=\sum_{k=0}^{N-1}s_k w_N^{km}$$

La première séparation peut donc s'écrire à partir des indices paires $(2p)$ et impaires $(2p+1)$ :

$$S_m=\sum_{p=0}^{N/2-1}s_{2p} e^{-j\frac{2\pi}{N}2pm}+\sum_{p=0}^{N/2-1}s_{2p+1} e^{-j\frac{2\pi}{N}(2p+1)m}  $$


On peut factoriser un multiplicateur commun de la seconde somme. 

$\sum_{p=0}^{N/2-1}s_{2p+1} e^{-j\frac{2\pi}{N}(2p+1)m}=e^{-j\frac{2\pi}{N}m}\sum_{p=0}^{N/2-1}s_{2p+1} e^{-j\frac{2\pi}{N}2pm}  $



$${\rm d'où} \qquad S_m=\underbrace{\sum_{p=0}^{N/2-1}s_{2p} e^{-j\frac{2\pi}{N/2}pm}}_{\text{TDF des index pairs } s_{2p}}+e^{-j\frac{2\pi}{N}m}\underbrace{\sum_{p=0}^{N/2-1}s_{2p+1} e^{-j\frac{2\pi}{N/2}pm}}_{\text{TDF des index impairs } s_{2p+1}}  $$

Si on pose $e_p=s_{2p}$, $o_p=s_{2p+1}$ et $P=\frac{N}{2}$, il est apparait que les deux sommes sont la TFD de la partie à index pair et la TFD de la partie à index impair de la somme initiale.

$$S_m=\underbrace{\sum_{p=0}^{P-1}e_{p} e^{-j\frac{2\pi}{P}pm}}_{\text{TDF des index pairs } s_{2p}=e_p}+e^{-j\frac{2\pi}{N}m}\underbrace{\sum_{p=0}^{P-1}o_{p} e^{-j\frac{2\pi}{P}pm}}_{\text{TDF des index impairs } s_{2p+1}=o_p}  $$ 


\begin{center}
\begin{tabular}{|l|cccccccc|}
 \hline 
signal $\{s_k\}$& \{$s_0$,&$s_1$,&$s_2$,&$s_3$,&$s_4$,&$s_5$,&$s_6$,&$s_7$\} \\
indice $k$ & 0&1&2&3&4&5&6&7 \\  
 \hline 
division en 2 signaux & \{$s_0$,&$s_2$,&$s_4$,&$s_6$\}&\{$s_1$,&$s_3$,&$s_5$,&$s_7$\} \\
nouvel indice $p$ & 0&1&2&3&0&1&2&3 \\ 

nouveaux signaux $\{e_p\}$ et $\{o_p\}$  & \{$e_0$,&$e_1$,&$e_2$,&$e_3$\},&\{$o_1$,&$o_2$,&$o_3$,&$o_4$\} \\
 \hline 
 \end{tabular}  
\end{center}


On note la TFD des signaux d'index pair $E_m$ (Even-index pour index pairs) et la TFD des entrées d'index impair $O_m$ (Odd-index pour index impairs) et on obtient : 
 
$$S_m=E_m+e^{-j\frac{2\pi}{N}m}O_m  $$


\hrulefill

		\quest La particularité de la définition de la TFD, nous permet de ne faire que la moitié des calculs, soit $N/2$ au lieu de $N$. En effet, montrer que 
		
		$$S_{m+\frac{N}{2}}	=E_m-e^{-j\frac{2\pi}{N}m}O_m$$ 
		
		
\hrulefill
		
		La particularité de la définition de la TFD, nous permet de ne faire que la moitié des calculs, soit $N/2$ au lieu de $N$. Explication.

Calculons 

$$S_{m+\frac{N}{2}}=\sum_{p=0}^{N/2-1}s_{2p} e^{-j\frac{2\pi}{N}2p({m+\frac{N}{2}})}+\sum_{p=0}^{N/2-1}s_{2p+1} e^{-j\frac{2\pi}{N}(2p+1)({m+\frac{N}{2}})}  $$


Soit, en procédant comme précédemment :

\begin{align*}
S_{m+\frac{N}{2}}	&=\sum_{p=0}^{N/2-1}s_{2p} e^{-j\frac{2\pi}{N}2p({m+\frac{N}{2}})}+e^{-j\frac{2\pi}{N}({m+\frac{N}{2}})}\sum_{p=0}^{N/2-1}s_{2p+1} e^{-j\frac{2\pi}{N}2p({m+\frac{N}{2}})} \\
					&=\sum_{p=0}^{N/2-1}s_{2p} e^{-j\frac{2\pi}{N}2pm}e^{-j2\pi p}+e^{-j\frac{2\pi}{N}m}e^{-j\pi}\sum_{p=0}^{N/2-1}s_{2p+1} e^{-j\frac{2\pi}{N/2}pm}e^{-j2\pi p}\\
					&=\sum_{p=0}^{N/2-1}s_{2p} e^{-j\frac{2\pi}{N}2pm}e^{-j2\pi p}-e^{-j\frac{2\pi}{N}m}\sum_{p=0}^{N/2-1}s_{2p+1} e^{-j\frac{2\pi}{N/2}pm}\\
					&=E_m-e^{-j\frac{2\pi}{N}m}O_m 
\end{align*}


\hrulefill

Finalement, le seul calcul des valeurs $E_m$ et $O_m$, sur la moitié des valeurs, soit $\frac{N}{2}$, permet d'obtenir la totalité du spectre $S_m$. 



\begin{center}
\[\left(
   \begin{array}{l c l}
      S_m  & = & E_m+e^{-j\frac{2\pi}{N}m}O_m \qquad \text{ pour } 0\leq m \leq \frac{N}{2}-1\\
      S_{m+\frac{N}{2}}   & = & E_m-e^{-j\frac{2\pi}{N}m}O_m \qquad \text{ pour } \frac{N}{2}\leq m \leq N-1\\
   \end{array}
\right.
\]
\end{center}


Ce résultat exprime récursivement la TFD de longueur $N$ en deux TFD de taille $\frac{N}{2}$ ce qui est au cœur de la transformée de Fourier discrète rapide dite FFT.

\begin{center}
\includegraphics[scale=1]{Fig_FFT_Recursif}
\caption{Schéma de principe d'un étage récursif}
\end{center}
		
		
		\subsection*{Complexité}

 	\quest Quelle est la complexité de ce nouvel algorithme ?
 	
 
\hrulefill

Partons des deux formulations vues précédemment et rappelons-nous que pour l'algorithme naïf, pour un vecteur $N$, on a $\mathcal{O}(N^2)$.

	$$S_m  =  E_m+e^{-j\frac{2\pi}{N}m}O_m$$
			\paragraph{Pour le premier étage récursif}
			
\begin{enumerate}
\item Complexité de $E_m$ qui met en jeu $N/2$ éléments est en $\mathcal{O}_1\left(\left(\frac{N}{2}\right)^2\right)$ 
\item De même pour $O_m$
\item Ce à quoi il faut ajouter $N$ multiplications pour  $e^{-j\frac{2\pi}{N}m} \times O_m $  
\end{enumerate}

La complexité de ce premier étage récursif est donc : $$\mathcal{C}_1(N)=2.\mathcal{O}_1\left(\left(\frac{N}{2}\right)^2\right)+\mathcal{O}_1(N)$$

Soit $$\mathcal{C}_1(N)=2\textcolor{red}{\left(\frac{N}{2}\right)^2}+N=\frac{N^2}{4}+N$$ 
			\paragraph{Pour le second étage récursif}

Pour ce second étage, les complexités de $E_m$ et $O_m$ sont réduites puisqu'il n'y a que $P=\frac{N}{2}$ valeurs. 

$$\textcolor{red}{\left(\frac{N}{2}\right)^2}\leftarrow \textcolor{blue}{\left[2\left(\frac{N}{4}\right)^2+\frac{N}{2}\right]}$$

La complexité de s'écrit alors 
\begin{align*}
\mathcal{C}_2(N)	&=2\textcolor{blue}{\left[2\left(\frac{N}{4}\right)^2+\frac{N}{2}\right]}+N \\
					&=4\left(\frac{N}{4}\right)^2+N+N \\
					&=\frac{N^2}{2^2}+2N
\end{align*}


\paragraph{Pour le troisième étage récursif}

Pour ce troisième étage, les complexités de $E_m$ et $O_m$ sont réduites de la même façon : $$\textcolor{blue}{\left(\frac{N}{4}\right)^2}\leftarrow \textcolor{purple}{\left\{2\left(\frac{N}{8}\right)^2+\frac{N}{4}\right\}}$$

 
\begin{align*}
\mathcal{C}_3(N)	&=2\textcolor{blue}{\Bigg[2\frac{}{}}\textcolor{purple}{\left\{2\left(\frac{N}{8}\right)^2+\frac{N}{4}\right\}}\textcolor{blue}{\left.+\frac{N}{2}\right]}+N	\\
				&= 4\left[2\left(\frac{N}{8}\right)^2+\frac{N}{4}\right]+N+N\\
				&= 8\left(\frac{N}{8}\right)^2+N+N+N\\
				&= \frac{N^2}{8}+3N\\
				&= \frac{N^2}{2^3}+3N\\
\end{align*}

\paragraph{Pour le n-ième étage récursif}
 
\begin{align*}
\mathcal{C}_n(N)	&= \frac{N^2}{2^n}+nN\\
\end{align*}

Comme $N=2^n$, et que $\log_2(N)=n.\log_2(2)$, soit $n=\log_2(N)$, 

on obtient finalement 


$$\mathcal{C}_n(N)	= N+N.\log_2(N)\simeq \mathcal{O}(N.\log_2(N)) $$


\hrulefill

La FFT ramène le calcul de la transformée de Fourier discrète de $\mathcal{O}(N^2)$ à $\mathcal{O}(N\log(N))$. Cette réduction de complexité a ouvert la résolution de nombreux problèmes. 

Pour poursuivre l'exemple du son, une centaine de milliers
d'opérations seront mobilisées pour la FFT, au lieu de 2 milliards donc pour le TFD naïve : 10000 fois plus rapide.\\





		\subsection*{Pseudo-code de l'algorithme récursif de la FFT}
		
\begin{python}
1. Fonction FFT(s)
2. 		N = longueur de s # s est une puissance de 2 
3.		si N = 1
4.			alors retourner s
5.		sinon 
6.			e = liste des s a indice pair 
7.			o = liste des s a indice impair
8.			E = la FFT de la liste e
9.			O = la FFT de la liste o
10.			pour m de 0 a N/2-1 faire
11.				S[m] 	 = E[m] + exp(-j*2*\pi /N*m)*O[m]
12.				S[m+N/2] = E[m] - exp(-j*2*\pi /N*m)*O[m]
13.		retourner S
\end{python}

\begin{enumerate}
	\item La ligne \pyth{2.} calcule la longueur de la liste (type \pyth{array}) de \pyth{N} qui doit être une puissance d'un entier $n$ : $N=2^n$.
	\item Les lignes \pyth{3.} et \pyth{4.} renvoient \pyth{s} car pour un seul élément, la transformée de Fourier est l'élément lui-même. 
	\item Les lignes \pyth{6.}  et \pyth{7.} font la division de la liste initiales à partir de la parités des indices des grandeurs du signal dans la liste \pyth{s}.
	\item Les lignes \pyth{8.}  et \pyth{9.} sont proprement récursives.
	\item Enfin, les lignes \pyth{11.}  et \pyth{12.} construisent le vecteur spectre \pyth{S} complexe, à partir de la reformulation donnée plus haut.
	\item La dernière ligne renvoie le spectre complexe.
\end{enumerate}


Un algorithme qui divise les objets à traiter, qui traite les sous-objets puis les recombine, est appelé un algorithme \og diviser pour régner \fg.


		\subsection*{Code en Python}		
Tous les listes \pyth{numpy} que nous considérerons
auront une taille en puissance de 2 : $N=2^n$. Cette condition est assez restrictive même s'il y a des moyens de la lever.

Comme nous travaillerons avec des listes \pyth{numpy}. Il sera inutile de faire des vecteurs colonnes.

			\paragraph{Diviser}

Nous devons scinder en deux moitiés les $s_k$ du signal en deux listes  constituées respectivement des éléments du signal d'indices pairs $s_{2p}=e_p$ et des éléments d'indices impairs $s_{2p+1}=o_p$.

	\quest Définir une fonction \pyth{diviser(s)} qui prend en argument un signal \pyth{s} (liste \pyth{numpy}) et renvoie \pyth{e} et \pyth{s} deux listes  constituées respectivement des éléments du signal d'indices pairs $s_{2p}=e_p$ et des éléments d'indices impairs $s_{2p+1}=o_p$

\hrulefill

Par exemple.

\begin{python}
def diviser(s):
	'''Prend un signal (liste np.array) en entree
    Renvoie deux listes np.array : e indice pair, o indice impair'''
    N = len(s)
    e = []
    o = []
    for i in range(N):
        if i%2 == 0:
            e.append(s[i])
        else:
            o.append(s[i])
    return np.array(e) , np.array(o)

\end{python}

ou bien 

\begin{python}
def diviser(s):
	'''Prend un signal (liste np.array) en entree
    Renvoie deux listes np.array : e indice pair, o indice impair'''
    e = s[::2]
    o = s[1::2]
    return e,o
\end{python}


ou encore

\begin{python}
def diviser(s):
    '''Prend un signal (liste np.array) en entree
    Renvoie deux listes np.array : e indice pair, o indice impair'''
    e = np.array([s[i] for i in range(len(s)) if i % 2 == 0])
    o = np.array([s[i] for i in range(len(s)) if i % 2 == 1])
    return e,o
\end{python}

\hrulefill

		\paragraph{Algorithme de la FFT}
	
\quest Définir la fonction \pyth{FFT(s)} qui prend en entrée un signal \pyth{s} et renvoie son spectre \pyth{S}. 

\hrulefill
			
Effectuer récursivement les FFT des deux moitiés en divisant en deux moitiés les deux moitiés précédentes, jusqu'à n'obtenir que des "moitiés" à un seul élément. On est alors ramené au cas où $N=1$. Il s'agit alors, après avoir mis en mémoire, après avoir empilés les différents calculs non exécutés en mémoire, de les calculer effectivement, à l'envers, à partir du cas où $N=1$.

 

\begin{python}
def FFT(s):
    N = len(s) #N est une puissance de 2
    S = np.zeros((N),dtype=complex)
    if N == 1:
        return np.array([s[0]+0.j])
    else:
        e,o = diviser(s)
        E = FFT(e)
        O = FFT(o)
        for m in range(N//2):
            S[m] = E[m] + np.exp((-2j*np.pi*m)/(N))*O[m]
            S[m+N//2] = E[m] - np.exp((-2j*np.pi*m)/(N))*O[m]
    return S
\end{python}	

\hrulefill

			\paragraph{Test}

	\quest Vérifier que la FFT donne $\textit{S}=\{5,3-2j,-3,3+2j\}$  à partir du signal $\textit{s}=\{2,3,-1,1\}$.



	\subsection*{FFT inverse}
	
		
Le travail pour la FFT peut être refait rapidement pour concevoir une FFT inverse.


A partir de $$s_k=\frac{1}{N}\sum_{m=0}^{N-1}S_m e^{j\frac{2\pi}{N}mk}$$

Si on note la TFDI des spectres d'index pair $e_k$ (Even-index pour index pairs) et la TFDI des entrées d'index impair $o_k$ (Odd-index pour index impairs) et on obtient : 
 
$$s_k=\frac{1}{2}(e_k+e^{j\frac{2\pi}{N}k}o_k)  $$

De même, on peut s'épargner le calculs de la totalité des composantes du signal en calculant

$$s_{k+\frac{N}{2}}=\frac{1}{2}(e_k-e^{j\frac{2\pi}{N}k}o_k)  $$


Finalement
\begin{center}
\[\left(
   \begin{array}{l c l}
      s_k=\frac{1}{2}(e_k+e^{j\frac{2\pi}{N}k}o_k) \qquad \text{ pour } 0\leq k \leq \frac{N}{2}-1\\
      s_{k+\frac{N}{2}}=\frac{1}{2}(e_k-e^{j\frac{2\pi}{N}k}o_k) \qquad \text{ pour } \frac{N}{2}\leq k \leq N-1\\
   \end{array}
\right.
\]
\end{center}


\hrulefill

Rappelons la définition de la transformée de Fourier discrète inverse (TFDI), 

$$s_k=\frac{1}{N}\sum_{m=0}^{N-1}S_m e^{j\frac{2\pi}{N}mk}$$

La première séparation peut donc s'écrire à partir des indices paires $(2p)$ et impaires $(2p+1)$ :

$$s_k=\frac{1}{N}\sum_{p=0}^{N/2-1}s_{2p} e^{j\frac{2\pi}{N}2pk}+\frac{1}{N}\sum_{p=0}^{N/2-1}s_{2p+1} e^{j\frac{2\pi}{N}(2p+1)k}  $$


On peut factoriser un multiplicateur commun de la seconde somme. 

$\frac{1}{N}\sum_{p=0}^{N/2-1}s_{2p+1} e^{j\frac{2\pi}{N}(2p+1)k}=e^{j\frac{2\pi}{N}k}\frac{1}{N}\sum_{p=0}^{N/2-1}s_{2p+1} e^{j\frac{2\pi}{N}2pk}  $



$${\rm d'où} \qquad s_k=\frac{1}{2}.\underbrace{\frac{2}{N}\sum_{p=0}^{N/2-1}s_{2p} e^{j\frac{2\pi}{N}2pk}}_{\text{TDFI des index pairs}}+e^{j\frac{2\pi}{N}k}\frac{1}{2}.\underbrace{\frac{2}{N}\sum_{p=0}^{N/2-1}s_{2p+1} e^{j\frac{2\pi}{N}2pk}}_{\text{TDFI des index impairs}}  
$$



Si on pose $E_p=S_{2p}$, $O_p=S_{2p+1}$ et $P=\frac{N}{2}$, il est apparait que les deux sommes sont la TFDI de la partie à index pair et la TFDI de la partie à index impair de la somme initiale.

$$s_k=\frac{1}{2}\left(\underbrace{\frac{1}{P}\sum_{p=0}^{P-1}E_{p} e^{j\frac{2\pi}{P}pk}}_{\text{TDFI des index pairs}\; E_p=S_{2p}}+e^{j\frac{2\pi}{N}k}\underbrace{\frac{1}{P}\sum_{p=0}^{P-1}O_p e^{j\frac{2\pi}{P}pk}}_{\text{TDFI des index impairs}\; O_p=S_{2p+1}} \right) $$


On note la TFDI des spectres d'index pair $e_k$ (Even-index pour index pairs) et la TFDI des entrées d'index impair $o_k$ (Odd-index pour index impairs) et on obtient : 
 
$$s_k=\frac{1}{2}(e_k+e^{j\frac{2\pi}{N}k}o_k)  $$

De même, on peut s'épargner le calculs de la totalité des composantes du signal en calculant

$$s_{k+\frac{N}{2}}=\frac{1}{2}(e_k-e^{j\frac{2\pi}{N}k}o_k)  $$


Finalement
\begin{center}
\[\left(
   \begin{array}{l c l}
      s_k=\frac{1}{2}(e_k+e^{j\frac{2\pi}{N}k}o_k) \qquad \text{ pour } 0\leq k \leq \frac{N}{2}-1\\
      s_{k+\frac{N}{2}}=\frac{1}{2}(e_k-e^{j\frac{2\pi}{N}k}o_k) \qquad \text{ pour } \frac{N}{2}\leq k \leq N-1\\
   \end{array}
\right.
\]
\end{center}

\hrulefill


		\quest Définir la fonction \pyth{FFTI(S)} qui prend en entrée un spectre  \pyth{S} et renvoie son signal complexe \pyth{s}. 
		
\hrulefill		
		
\begin{python}
def FFTI(S):
    N = len(S) #N est une puissance de 2
    s = np.zeros((N),dtype=complex)
    if N == 1:
        return np.array([S[0]+0.j])
    else:
        E,O = diviser(S)
        e = FFTI(E)
        o = FFTI(O)
        for k in range(N//2):
            s[k] = 1/2*(e[k] + np.exp((2j*np.pi*k)/(N))*o[k])
            s[k+N//2] = 1/2*(e[k] - np.exp((2j*np.pi*k)/(N))*o[k])
    return s
\end{python}
		
\hrulefill


\section{Application à l'analyse d'un son}

Les sons que vous allez traiter sont des signaux de codes DTMF (Dual Tone Multiple Frequency) utilisés pour distinguer par leur doublet de fréquences les touches numérotées d'un clavier téléphonique (DTMF sur wikipedia).

L'appuie sur une touche génère donc un son composé de deux fréquences, selon le tableau suivant :

\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline 
 & 1 209 Hz  & 1 336 Hz  & 1 477 Hz & 1 633 Hz \\ 
\hline 
697 Hz & 1 & 2 & 3 & A \\ 
\hline 
770 Hz & 4 & 5 & 6 & B \\ 
\hline 
852 Hz & 7 & 8 & 9 & C \\ 
\hline 
941 Hz & * & 0 & \verb'#' & D \\ 
\hline 
\end{tabular} 

\includegraphics[width=5cm]{Clavier_DTMF.png}
\end{center}


		\subsection*{Travailler avec un fichier son sous Python}

Deux modules sont pratiques pour lire et extraire les données importantes d'un fichier son au format \pyth{.wav}.


\begin{python}
import soundfile as sf
import sounddevice as sd
\end{python}
			\subsection*{Lire et extraire des données du fichier \pyth{son.wav}}

\begin{python}
s, fe = sf.read(r'son.wav')
#s est le signal
#fe la frequence d echantillonnage du signal
#Le nombre d echantillon est donne par s.shape[0] pour la voie 0 
N = s.shape[0]
\end{python}
			\subsection*{Écouter le fichier \pyth{son.wav}}
\begin{python}
#Pour ecouter le son
sd.play(s,fe)
\end{python}		
	

			\subsection*{Analyser le son \pyth{son1.wav}}

				\subparagraph{L'échantillonnage est-il en $2^n$ ?\\}

Le nombre d'échantillons $N\neq 2^n$. Ici, \pyth{N=22051}. On va  ajouter des '0' pour avoir $N=2^n$. Initialement $p$ est tel que $N=2^p$ avec $p$ non entier. On a donc $\log_2{(N)}=p.\underbrace{\log_2{(2)}}_{=1}$. Dans notre exemple, $p=14.42855646219518$
		
		On prend pour $n$ la valeur entière immédiatement supérieure : $$n={\rm int}(p)+1={\rm int}(\log_2(N))+1$$
 
 soit finalement, $n=15$. 
 		
		Ajoutons ensuite au signal les zéros manquants. 

	\quest Construire une fonction \pyth{puissance2(liste)} 
	\begin{enumerate}
		\item qui prend en entrée une liste  \pyth{liste}
		\item et renvoie 
		\begin{enumerate}
			\item le nombre de points en $2^n$ avec $n$ l'entier supérieur à $p$ de la nouvelle liste
			\item ainsi que la nouvelle liste des valeurs du signal complétée des '0' nécessaires.	
		\end{enumerate} 
	\end{enumerate} 	
		
\hrulefill 		

\begin{python}
def puissance2(liste):
    N = len(liste) 
    N_nouv = 2**(int(np.log2(N))+1)
    N_zeros=np.zeros(N_nouv-N)
    #Ajouter des N_zeros 0 pour avoir N=2^n
    return N_nouv,np.concatenate((liste,N_zeros))

N_nouv,s_nouv = puissance2(s)
\end{python}

\hrulefill 		


				\subparagraph{De quelle touche s'agit-il ?}
				 
			\subsection*{Analyser le son \pyth{son3.wav}}
\quest En prenant les mêmes précautions que précédemment, indiquer de quelle(s) touche(s) il s'agit.

			\subsection*{Analyser le son \pyth{son4.wav}}
\quest Déduire de l'analyse spectrale le numéro de téléphone tapé sur le clavier DTMF?

\hrulefill

\begin{enumerate}
	\item Les premières difficultés
	\begin{enumerate}
		\item Le nombre d'échantillons $N\neq 2^n$. Ici, \pyth{N=55459}.
		\item A l'écoute, le signal est une succession de 10 sons. 
	\end{enumerate} 
	\item Les solutions envisagées
	\begin{enumerate}
		\item Ajouter des '0' pour avoir $N=2^n$ en se servant de \pyth{puissance(liste)}
		
Faisons la FFT du signal complété des zéros adéquates. 

\begin{center}
\includegraphics[scale=0.8]{Fig_son4.pdf}
\caption{Signal complété de '0' pour avoir $N=2^n$. }
\end{center}


		\item Travailler en isolant les différents signaux. On prendra des tranches du signal initial -- 10 tranches donc -- de largeur le nombre d'échantillons initial divisé par dix : $N/10$.
		
\begin{python}
#on travaille sur une tranche (le signal est coupe en 10 tranches)
tranche=int(len(s)/10)
N=tranche
k = 0
s=s[k*tranche:(k+1)*tranche]

\end{python}			
	\end{enumerate}
\end{enumerate}

\begin{center}
\includegraphics[scale=0.7]{Fig_son4_tranche1.pdf}
\caption{Pour la première tranche : touche \pyth{0 (941 Hz, 1336 Hz)}}
\end{center}


Nous avons les touches : 

\begin{itemize}
	\item 0 (941 Hz, 1336 Hz)
  	\item 3 (694 Hz, 1476 Hz)
  	\item 2 (698 Hz, 1330 Hz) 
  	\item 0 (941 Hz, 1336 Hz)
 	\item 6 (771 Hz, 1476 Hz) 	
 	\item 0 (941 Hz, 1336 Hz)
 	\item 5 (770 Hz, 1336 Hz)
 	\item 0 (941 Hz, 1336 Hz)
\end{itemize} 

\begin{center}
\includegraphics[scale=0.7]{Fig_son4_tranche2.pdf}
\caption{Pour la seconde tranche : touche \pyth{3 (694 Hz, 1476 Hz)}}
\end{center}


\newpage{}
\section{Annexe : \pyth{numpy}}


	\subsection*{Créer une matrice de zéros complexes}
	
\begin{center}
\pyth{M = np.zeros((n,m),dtype = complex)}
\end{center}

\begin{python}
>>> np.zeros((3,2),dtype = complex)
array([[ 0.+0.j,  0.+0.j],
       [ 0.+0.j,  0.+0.j],
       [ 0.+0.j,  0.+0.j]])
\end{python}

	\subsection*{Accéder à une ligne, une colonne, une valeur d'une matrice}
\begin{python}
	
>>> A=np.array([[1.0-2.j,3.0+1.j],[0.1+0.2j,0.0+0.0j],
       [ 0.+0.j,  0.+0.j]])
>>> A
array([[1. -2.j , 3. +1.j ],
       [0.1+0.2j, 0. +0.j ],
       [0. +0.j , 0. +0.j ]])
\end{python}

\pyth{A[1,:]} renvoie les valeurs de toutes les colonnes de la ligne d'indice 1 (la seconde ligne), donc la seconde ligne.
\begin{python}
>>> A[1,:]
array([0.1+0.2j, 0. +0.j ])
\end{python}

\pyth{A[:,0]} renvoie les valeurs de toutes les lignes de la colonne d'indice 0 (la première colonne), donc la première colonne.
\begin{python}
>>>A[:,0]
array([1. -2.j , 0.1+0.2j, 0. +0.j ])
\end{python}

\pyth{A[0,1]} renvoie la valeur de la lignes d'indice 0 et de la colonne d'indice 1.
\begin{python}
>>> A[0,1]
(3+1j)
\end{python}


	\subsection*{Produit matricielle}
Si \pyth{A} est une matrice $(n \times m)$ et que \pyth{B} une matrice de taille $(m\times p)$, le produit matriciel \pyth{AB} s'écrit avec \pyth{numpy}.

\begin{center}
\pyth{np.dot(A,B)}
\end{center}

\begin{python}
>>> A=np.array([[1.0-2.j,3.0+1.j],[0.1+0.2j,0.0+0.0j]])
>>> A
array([[1. -2.j , 3. +1.j ],
       [0.1+0.2j, 0. +0.j ]])
       
>>> B = np.array([[1.5-2.5j],[1.0+0.7j]])
>>> B
array([[1.5-2.5j],
       [1. +0.7j]])
       
>>> C = np.dot(A,B)
>>> C
array([[-1.2 -2.4j ],
       [ 0.65+0.05j]])
\end{python}

	\subsection*{Matrice conjuguée}
	
Si $\mathcal{A}$ est une matrice, sa matrice conjuguée $\overline{\mathcal{A}}$ est donnée par la méthode \pyth{conjugate()}

\begin{center}
\pyth{A.conjugate()}
\end{center}

	\subsection*{Ajouter des éléments à une liste \pyth{numpy}}


Attention, ici, \pyth{A} et \pyth{B} sont des listes \pyth{numpy} et non des tableaux à une ligne. 

\begin{center}
\pyth{np.concatenate((A,B))}
\end{center}
	
\begin{python}

>>> A=np.array([1.0-2.j,3.0+1.j,1.5-2.5j,0.0+0.j])
>>> A
array([1. -2.j , 3. +1.j , 1.5-2.5j, 0. +0.j ])
>>> B = np.array([1.5-2.5j,1.0+0.7j])
>>> B
array([1.5-2.5j, 1. +0.7j])
>>> C = np.concatenate((A,B))
>>> C
array([1. -2.j , 3. +1.j , 1.5-2.5j, 0. +0.j , 1.5-2.5j, 1. +0.7j])
\end{python}


	\subsection*{Créer une liste \pyth{numpy}}

\begin{center}
\pyth{np.arange(debut,fin,pas)}
\end{center}

\begin{python}
>>> L=np.arange(0,2,0.3)
>>> L
array([0. , 0.3, 0.6, 0.9, 1.2, 1.5, 1.8])
\end{python}

\end{document}

\thispagestyle{empty}