
\chapter{Le cas discret}

\section{Transformée de Fourier Discrète, applications}
\label{s:TFD}\setcounter{exo}{0}

\subsection{La dimension 1}

La transformée de Fourier discrète est un moyen de calculer les coefficients de Fourier d'une fonction $a$-périodique
$u$ à partir de ses $N$ échantillons $u(\frac{ka}N)$, $k=0, ... N\!-\!1$.  Cela n'est possible que si la fonction
présente un nombre de fréquences inférieur ou égal à $N$. Pour coller à la pratique numérique, nous supposerons
toujours dans ce chapitre que $N$ est pair. En général, $N$ est une puissance de 2. Tous les résultats énoncés (sauf
ceux du paragraphe \ref{TFR} relatifs à la transformée de Fourier rapide) s'adaptent toutefois sans difficulté au cas
$N$ impair.  Soit $u(x)$ une fonction réelle ou complexe  de p{é}riode $a$, et $N$ un entier pair. On cherche un
polyn{\^o}me trigonom{é}trique de la forme
\begin{equation}
\label{poly} P(x)=\sum_{n=-\frac N2}^{\frac N2-1}\tilde{u}_n \exp\left(\frac{2i\pi n x}{a}\right),
\end{equation}
qui soit {é}gal {à} $u$ aux points $\frac{ka}{N}$ pour $k=0, ..., N\!-\!1$. On dira dans la suite que $P$ est de degré
$\frac N2$.  Le but est donc d'interpoler les échantillons $u(\frac{ka}{N})=u_k$.
\newline\newline
{\small Pourquoi choisir un polynôme trigonométrique ? La raison est physique : tous les dispositifs d'acquisition de
signaux (sons) ou images ont une {\it bande passante}, \ccaadd un intervalle de fréquences captées par le dispositif
d'enregistrement ; les autres fréquences sont perdues ou tellement atténuées qu'on les néglige : on suppose donc ici
que la "bande passante" est contenue dans l'intervalle $[-\frac N2, \frac N2-1]$.  Il n'y a par contre aucune raison de
supposer que le signal ou l'image soit périodique, et en plus d'une période qui coïncide avec la fenêtre d'observation
$[0, a]$ comme nous sommes en train de le suppposer pour $P$. Cette hypothèse est donc imposée à la donnée par une
périodisation abusive, et provoque une distorsion près des bords de l'image que l'on peut voir et évaluer : c'est le
phénomène de Gibbs. Si  la fonction $u$ dont on possède les $N$ échantillons n'avait pas une bande de fréquence
contenue dans $[-\frac N 2, \frac N2-1]$, son interpolation par un polynôme trigonométrique de degré $\frac N2$
provoque une autre distorsion, très grave, que nous allons évaluer précisément dans ce chapitre : l'aliasage.}
\newline
On va commencer par calculer les coefficients de $P$.

\dexer  On pose $\omega_N=\exp\left(\frac{2i\pi}{N}\right)$, racine $N$-ième de l'unité. Montrer que
$\sum_{k=0}^{N\!-\!1}\omega_N^k=0$, puis que $\sum_{k=0}^{N\!-\!1}\omega_N^{kl}=0$ pour $l\neq 0$ modulo $N$ et
finalement que pour tout $k_0$, $\sum_{k=k_0}^{k_0+N\!-\!1}\omega_N^{kl}=0$ pour tout $l\neq 0$ modulo $N$.

\fexer
\begin{Def} \index{A}{transformée de Fourier!discrète}
 On pose
$u_k=u(\frac{ka}{N})$ et, pour $n=-\frac N2, ...,\frac N2-1$,
\begin{equation}
\label{defTFD} \tilde u_n = \frac{1}{N}\sum_{l=0}^{N\!-\!1}u_l\omega_N^{-nl}.
\end{equation}
Les $N$ coefficients $\tilde u_n$ sont appelés {\rm transformée de Fourier discrète} (TFD) des $N$ échantillons $u_k$.
On appelle transformée de Fourier discrète inverse l'application de $\CC^N$ dans lui même définie par
\begin{equation}\label{TFDinv}
u_k=\sum_{n=-\frac N2}^{\frac N2-1} \tilde u_n \omega_N^{kn}, \; k=0, ..., N\!-\!1.
\end{equation}
\end{Def}
\begin{pro} \index{A}{transformée de Fourier!discrète}
\label{coeff} Les coefficients $(\tilde u_n)$ d{é}finis par (\ref{defTFD}) sont les uniques coefficients tels que le
polyn{\^o}me trigonom{é}trique (\ref{poly}) v{é}rifie $P\left(\frac{ka}{N}\right)=u_k$,  pour tout $k=0, ..., N\!-\!1$.
En d'autres termes, la transformée de Fourier discrète composée avec son inverse donne bien l'identité.
\end{pro}
\paragraph{D{é}monstration}
Pour $k=0, ..., N\!-\!1$,
\begin{eqnarray*}
P(\frac{ka}{N})&=&\sum_{n=-\frac N2}^{\frac N2-1}\tilde u_n \omega_N^{nk}\\
&=&\frac{1}{N}\sum_{n=-\frac N2 }^{\frac N2-1}(\sum_{l=0}^{N\!-\!1}u_l
\omega_N^{-nl})\omega_N^{nk}\\
&=& \frac 1N \sum_{l=0}^{N\!-\!1}u_l(\sum_{n=-\frac N2}^{\frac N2-1}\omega_N^{nk-nl})\\
&=&\frac{1}{N}\sum_{l=0}^{N\!-\!1}N\delta(k-l)u_l = u_k,\\
\end{eqnarray*} o{\`u} on a noté $\delta$ la fonction
d{é}finie sur les entiers, valant 1 en 0, et 0 ailleurs. L'unicité provient du fait que toute application linéaire
surjective de $\complex^N$ dans lui-même   est aussi injective.\cqfd


\begin{cor}\index{A}{polynôme trigonométrique}
\label{commut-poly} Si $u$ est un polynôme trigonométrique $u(x)=\sum_{n=-\frac N2}^{\frac N2-1}\tilde{u}_n
\exp\left(\frac{2i\pi n x}{a}\right),$ les coefficients $\tilde u_n$ sont obtenus à partir des échantillons de $u$ par
la formule (\ref{defTFD}). Ce sont les coefficients de Fourier de $u$.
\end{cor}

\dexer  On note $u$ un vecteur de $\CC^N$ et $\hbox{TFD}(u)=\tilde u$ sa transformée de Fourier discrète. Vérifier que
$\sqrt{N}\;\hbox{TFD}$ est unitaire et calculer $\hbox{TFD}^{-1}$ grâce à la formule  $\hbox{TFD}^{-1}=
N\;^t\overline{\hbox{TFD}}.$

\fexer

   On rappelle
d'autre part que si $u \in L^2(0,a)$, les coefficients de la s{é}rie de Fourier de $u$ sont d{é}finis, pour $n \in \Z$,
par
\begin{equation}
\label{cn} c_n(u)=\frac{1}{a} \int_0^a u(x) \exp \left( \frac{-2i\pi nx}{a}\right).
\end{equation}


Les coefficients $\tilde u_n$ de la transformée de Fourier discrète sont approchés par les termes de la TFD de $(u_k)$
au sens suivant:

\begin{pro} Soit $u$ continue et $a$-périodique. Alors
les $\tilde u_n$ sont des approximations des $c_n(u)$ par la formule des trapèzes, pour
$n=\frac{-N}{2},...,\frac{N}{2}-1$.
\end{pro}

\paragraph{D{é}monstration}
Il suffit d'{é}crire l'approximation de l'int{é}grale (\ref{cn}) par la m{é}thode des trapèzes  en tenant compte du
fait que $u(a)=u(0)$ pour une fonction $a$-périodique.  \cqfd


\begin{figure}
\begin{center}
\psfrag{a}{$u$} \psfrag{b}{$u_k$} \psfrag{c}{$\hat u_k$} \psfrag{TFD}{TFD} \psfrag{Serie de Fourier}{Série de Fourier}
\psfrag{echantillonage}{échantillonnage}
\includegraphics[width=6cm]{Figures/diagTFD.eps}
\caption[]{La TFD après échantillonnage calcule bien les coefficients de Fourier {\bf si la fonction $u$ est un polynôme
trigonométrique} (corollaire \ref{commut-poly})}
\end{center}
\end{figure}

\begin{pro}
\label{sym1} On suppose que les échantillons $u_k$ sont réels. Alors
 $\tilde u_0$ et $\tilde u_{-\frac N2}$ sont r{é}els, et pour $k=1...\frac N2-1$,
$\tilde u_k=\overline{\tilde u_{-k}}$.
\end{pro}
\paragraph{D{é}monstration}
$\tilde u_0=\frac 1N\sum_k u_k$, et $\tilde u_{-\frac N2}=\frac 1N\sum (-1)^k u_k$ ; ces deux coefficients sont donc
r{é}els. D'autre part
$$\tilde u_{-n}=\frac{1}{N}\sum_{k=0}^{N\!-\!1}u_k
\omega_N^{kn}=\frac{1}{N}\sum_{k=0}^{N\!-\!1}u_k \overline{\omega_N^{-nk}}=\overline{\tilde u_n}.$$ \cqfd

  Remarquons  le
r{\^o}le particulier jou{é} par le terme $\tilde u_{-\frac N2}$, qui n'a pas de terme conjugu{é} lui correspondant.
\begin{pro}\label{coeffnulNsurdeux}
si $u$ est un polynôme trigonométrique réel dont les fréquences sont parmi $-\frac N2, ..., \frac N2 -1$, le terme
$\tilde u_{-\frac N2}$ est nul.
\end{pro}

\paragraph{Démonstration}
En effet, en regroupant les termes conjugués, on a, pour le polynôme trigonométrique P dont les coefficients sont les
$\tilde u_n$ :
$$
P(x)=\tilde u_0+\sum_{n=1}^{\frac N2-1} (\tilde u_n e^\frac{2in\pi x}a + \tilde u_{-n}e^\frac{-2in\pi x}a) +\tilde
u_{-\frac N2}e^\frac{-iN\pi x}a.
$$
Tous les termes de la somme sont réels sauf le dernier, qui ne l'est que si $\tilde u_{-\frac N2}=0$.
\enddemo

  La Figure \ref{aaah} montre un
exemple de signal (repr{é}sentant le son A) et le module de sa TFD.


\begin{figure}
\begin{center}

\includegraphics[width=10cm]{Figures/aaah.eps}

\end{center}
\caption[Transformée de Fourier Discrète du son a]{Haut: un signal correspondant à la voyelle "Ah" (le signal
représente la pression de l'air en fonction du temps); bas: module de la TFD (coefficients $|\tilde u|$, voir le
texte). On remarque que le module du spectre est symétrique, et qu'il existe trois pics importants correspondant aux
fréquences dominantes.} \label{aaah}
\end{figure}

\subsection{La dimension 2}


On consid{è}re un r{é}el $a$, une fonction $u$ de $\real^2$ dans $\real$, telle que $u(x+a,y+a)=u(x,y)$. On fixe à
nouveau un entier $N$, et l'on pose $u_{k,l}=u\left(\frac{ka}{N},\frac{la}{N}\right)$. On d{é}finit la TFD des
$u_{k,l}$ comme la suite des coefficients, pour $m,n \in \{-\frac N 2,...,\frac N 2 -1\}$,
\begin{equation}
\label{TFD2} \tilde{u}_{m,n}=\frac{1}{N^2}\sum_{k=0}^{N\!-\!1}\sum_{l=0}^{N\!-\!1}u_{k,l}\omega_N^{-mk}\omega_N^{-nl }.
\end{equation}

\dexer Montrer que la transformation ainsi d{é}finie est ``s{é}parable'', c'est-à-dire que le passage des $u_{k,l}$ aux
$\tilde u_{m,n}$ s'effectue par deux TFDs {à} une dimension successives.

\fexer
%\paragraph{Remarque}: la suite des coefficients $u_{k,l}$ repr{é}sente souvent
%en pratique une image num{é}rique. Rappelons qu'une telle image se
%pr{é}sente comme une matrice de nombres correspondant {à} la quantit{é} de
%lumi{è}re emmagasin{é}e par les capteurs d'un appareil optique, comme
%expliqu{é} à la fin du chapitre \ref{intlebesgue}. Nous attirons l'attention du lecteur sur le
%fait que dans ce chapitre nous avons consid{é}r{é} les valeurs $u_{k,l}$ comme
%les valeurs ponctuelles d'une certaine fonction $u$ de $L^2$, alors que pour
%illustrer le concept de martingale il s'agissait des moyennes locales de cette fonction de
%d{é}part, ce qui correspond au mécanisme des capteurs des appareils
%photos numériques, les CCD. Ceci ne pose pas de probl{è}me particulier, car la moyenne locale
%peut {\^e}tre vu comme une valeur ponctuelle de la fonction apr{è}s convolution
%avec une fonction {à} support born{é}e.

De m{\^e}me qu'en dimension 1, nous avons la propri{é}t{é} d'interpolation suivante:
\begin{pro}
Soient les coefficients $\tilde u_{m,n}$ définis, pour $m,n = -\frac N2,\; ...,\;\frac N2-1$, par (\ref{TFD2}).
Considérons le polynôme trigonométrique
$$P(x,y)=\sum_{m,n=-\frac N2}^{\frac N2-1}\tilde u_{m,n}\exp(\frac{2i\pi mx}{a})\exp(\frac{2i\pi
ny}{a}).$$ Les coefficients $\tilde u_{m,n}$ sont les seuls nombres complexes tels que, pour tout $k, l \in \{0, ...,
N\!-\!1\}$, on ait $P\left(\frac{ka}{N},\frac{la}{N}\right)=u\left(\frac{ka}{N},\frac{la}{N}\right)$. Par conséquent,
la transformée discrète inverse de $u_{k,l}\to \tilde u_{m,n}$ est donnée par le calcul du polynôme aux échantillons
$(\frac {ka}N, \frac {la}N)$, $0\leq k,l\leq N\!-\!1$ :
$$u(k,l)=P(\frac{ka}N, \frac {la}N)=\sum_{-\frac N2}^{\frac N2-1}\sum_{-\frac N2}^{\frac
N2-1}\tilde u_{m,n}\omega_N^{km+ln}.$$
\end{pro}

\dexer Montrer la proposition précédente. Le calcul est exactement le m{\^e}me qu'en dimension 1. De m{\^e}me qu'en
dimension 1, nous pouvons identifier un certain nombre de sym{é}tries  des $\tilde u_{m,n}$ si l'image est à valeurs
réelles.
 On suppose {à} nouveau que $N$ est pair. Montrer également la proposition suivante.
 \fexer


\begin{pro}\label{sym} Supposons que les échantillons $u_{k,l}$ soient réels.
Alors les coefficients  $\tilde u_{0,0}$, $\tilde u_{0,-\frac N2}$, $\tilde u_{-\frac N2,0}$, et $\tilde u_{-\frac N2,
-\frac N2}$ sont r{é}els; de plus
$$\forall m, n \in \{-\frac N2 +1, ...,\frac N2-1\} \hskip 1cm\tilde u_{m,n}=\overline{\tilde
u_{-m,-n}}$$
 \end{pro}

\dexer A nouveau, comme en dimension 1, les coefficients $(\tilde u_{m,n})$ correspondent aux fr{é}quences de l'image
$u$, ordonn{é}es des n{é}gatives aux positives. Plus pr{é}cis{é}ment, si $u \in L^1$ et que l'on d{é}finit les
coefficients de la s{é}rie de Fourier de $u$ par
$$c_{m,n}=\frac{1}{4 \pi ^2} \int_{[0,\;2\pi]^2}u(x,y) \exp \left( \frac{-2i\pi mx}{a}\right)\exp
\left( \frac{-2i\pi n y}{a}\right),$$ alors, pour $m, n = -\frac N2, ...,\; \frac N2-1$ les $\tilde u_{m, n}$ sont des
approximations des $c_{m,n}$ par la m{é}thode des trapèzes. \fexer
 La figure \ref{image-TFD}
présente une image et le logarithme du module de sa transform{é}e de Fourier discr{è}te (le logarithme est utilis{é}
car le module des TFD des images usuelles d{é}croit tr{è}s vite lorsque l'on s'{é}loigne des basses fr{é}quences).

%\begin{figure}
%\begin{center}
%\includegraphics[width=8cm]{Figures/aero.ps}
%\includegraphics[width=8cm]{Figures/aerofft2.ps}
%\end{center}
%\caption{Haut: une image num{é}rique de taille $256\times256$;
%bas: le logarithme du module de sa TFD, dans le sens de la
%d{é}finition des $\tilde u$ (voir le texte)} \label{image-TFD}
%\end{figure}

\begin{figure}
\begin{center}

\includegraphics[width=7cm]{Figures/aero.ps}
\includegraphics[width=7cm]{Figures/aerofft2.ps}

\end{center}
\caption[Transformée de Fourier Discrète d'une image]{Gauche: une image numérique de taille $256\times256$; droite:
application d'un changement de contraste logarithmique $s\to C+D\log s$ au module de sa TFD. (Les constantes $C$ et $D$
sont choisies pour que l'image résultante ait ses valeurs entre 0 et 255). Le spectre décroît rapidement aux hautes
fréquences (rappelons que l'image étant bornée, son spectre est dans $L^2$). C'est la grande vitesse de décroissance du
spectre qui rend nécessaire l'utilisation du logarithme pour la visualisation. La symétrie centrale du module de la TFD
est visible. Les lignes blanches horizontale et verticale correspondent aux bords verticaux et horizontaux de l'image,
respectivement. En effet, l'image étant implicitement périodisée elle présente un contraste brusque entre bord haut et
bord bas d'une part, coté droit et coté gauche d'autre part. Remarquer également les lignes obliques qui correspondent
aux bords obliques de l'image (voir en particulier les dalles sur le sol).} \label{image-TFD}
\end{figure}

\subsection{Le ph{é}nom{è}ne du repliement de spectre ou aliasage}
\index{A}{repliement de spectre}


Le but de ce paragraphe est de calculer les perturbations auxquelles est exposée la transform{é}e de Fourier discrète
d'un signal  lorsque celui-ci est sous-{é}chantillonn{é}. On vient de voir que la transformée de Fourier discrète
calculait exactement les coefficients de Fourier  d'un polynôme trigonométrique de degré $\frac N2$,
$P(x)=\sum_{n=-\frac N2}^{\frac N2-1}\tilde{u}_n \exp\left(\frac{2i\pi n x}{a}\right)$, dont on connaissait $N$
échantillons $u(\frac {ka}N)$, $N=0,..., N\!-\!1$. Dans cette section, on considère une fonction $u\in L^2(0,a)$ et sa
série de Fourier
$$u(x)= \sum_{n\in \relative} c_n (u)e^{\frac {2in\pi x}a}.$$
Dans toute la suite, on supposera que $\sum_{n\in \relative}|c_n(u)|< +\infty$, ce qui implique que $u$ est continue et
$a$-périodique.  Cette hypothèse n'est pas irréaliste. En effet, étant donné un signal $v$ régulier  ($C^2$ par
exemple) sur $[0, a/2]$, on peut le rendre pair en posant
 $u(-x)=\tilde v(x)$ pour $x\in [-a/2,0]$, $u(x)=v(x)$ sur $[0, a]$.  On voit que la
$a$-périodisée
 de cette extension $u$ reste Lipschitz et $C^2$ par morceaux et on peut en déduire (exercice
!) que
 la série des coefficients de Fourier de $u$ est convergente.
  On suppose également, ce qui est réaliste,
 qu'un signal $u$ n'est en fin de compte connu que par ses
échantillons sur $[0, a]$, $u(0), ..., u(\frac {N\!-\!1}{N}a)$.



\begin{thf}
\label{Shannondiscret} Soit $u$ définie sur $[0, a]$, vérifiant
 $\sum_n|c_n(u)|<+\infty$.
  \index{A}{Shannon !théorème discret de}
 Alors la transformée de Fourier discrète des $N$ échantillons de $u$ sur $[0,a]$ est
la $N$-périodisée de la suite des coefficients de Fourier de $u$ :
\begin{equation}
\label{aliasagediscret} \tilde u_n= \sum_{q=-\infty}^{+\infty} c_{n+qN}(u),\; \; n=-\frac N2,...,\frac N2-1.
\end{equation}
\end{thf}

\paragraph{Démonstration}
On rappelle la notation $\omega_N=e^{\frac{2i\pi}N}$ et $(\omega_N)^N=1$. Comme
$$u(x)=\sum_{m\in \relative} c_m(u)e^\frac{2im\pi x}{a}, $$ on a
$$u(\frac {ka}N)=\sum_{m\in \relative}c_m(u)\omega_N^{mk}.$$
On pose pour $m\in \relative$, $m=qN+n$, $-\frac N2\leq n\leq \frac N2-1$. En regroupant les termes de la série de
Fourier on obtient
$$u(\frac{ka}N)=\sum_{n=-\frac N2}^{\frac N2-1}
\left(\sum_{q=-\infty}^{+\infty} c_{n+qN}(u)\right)\omega_N^{nk}, \; k=0,..., N\!-\!1.$$ Mais on a aussi (formule
d'inversion de la transformée de Fourier discrète):
$$u(\frac{ka}N)=\sum_{n=-\frac N2}^{\frac N2-1}\tilde u_n\omega_N^{nk}, \; k=0,..., N\!-\!1.$$
Ces deux dernières formules définissent toutes deux la transformée de Fourier discrète et par identification on obtient
la formule de "repliement de spectre" (\ref{aliasagediscret}).

\enddemo
On a le résultat analogue en deux dimensions:


\begin{thf}
\label{Shannondiscret2D} Soit $u$ définie sur $[0, a]^2$, vérifiant
 $\sum_n|c_{m,n}(u)|<+\infty$.
  \index{A}{Shannon !théorème discret de}
 Alors la transformée de Fourier discrète de $u$ est
la $(N,N)$-périodisée de la suite des coefficients de Fourier de $u$ :
\begin{equation}
\label{aliasagediscret2D} \tilde u_{m,n}= \sum_{p,q=-\infty}^{+\infty} c_{m+pN,n+qN}(u),\; \; m, n=-\frac N2,...,\frac
N2-1.
\end{equation}
\end{thf}
\dexer Montrer le théorème \ref{Shannondiscret2D}.
 \fexer

Ce théorème va nous permettre d'interpréter les effets de moiré visibles dans beaucoup d'images digitales ou  de films
digitalisés (DVD). Ces effets de moiré sont dûs à un "repliement de spectre", ou "aliasage".  Le repliement de spectre
provient d'un sous-échantillonnage abusif.   Le terme aliasage se réfère à la présence des coefficients parasites
$c_{m+pN,n+qN}$, pour $(p,q)\neq (0,0)$ dans le calcul du coefficient de la fréquence $(m,n)$, $\tilde u_{m,n}$. Quand
la transformée de Fourier discrète fait correctement son travail, qui est de retrouver le coefficient $c_{m,n}$ de la
fréquence $(m,n)$ de $u$, on doit avoir $\tilde u_{m,n}=c_{m,n}$. Les coefficients $c_{m+pN, n+qN}$ qui s'y ajoutent
dans (\ref{aliasagediscret2D}) sont des répliques, ou "alias" de coefficients correspondant aux fréquences plus grandes
$(m+pN, n+qN)$, $(p,q)\neq (0,0)$. D'où le terme d'aliasage. Dans une expérience sur ordinateur, les signaux et images
sont toujours déjà échantillonnés.  Voyons donc maintenant comment interpréter l'expérience pratique suivante: un
signal est donné par ses échantillons, et on le sous-échantillonne.


\begin{Def} \index{A}{sous-échantillonnage d'un signal}
Soit un signal échantillonné $(u_k)$, $k=0,...,N\!-\!1$, et soit $p$ un entier divisant $N$. On d{é}finit l'op{é}rateur
``sous-{é}chantillonnage d'ordre $p$'' comme suit:
$$ S_p: \RR^{N} \longrightarrow \RR^{ N/p}$$
$$ (u_k)_{k=0,\;...,\; N\!-\!1} \longrightarrow (v_k)=(u_{kp})_{k=0, ..., N/p}.$$
Le signal $(v_k)$ est dit sous-{é}chantillonn{é} d'un facteur $p$.
\end{Def}
 Nous commen{\c c}ons par le cas, technologiquement classique, où $p=2$.

\begin{cor}
\label{sousech2} Soit $(v_k)=S_2((u_k))$ (on suppose que $\frac N2$ est pair). Alors $(\tilde v _n)$, la transform{é}e
de Fourier discr{è}te de $(v_k)$, s'écrit, pour $n=-\frac N4,...,\frac N4-1$,
\begin{equation}
\tilde v_n= \tilde u_n+\tilde u_{n-\frac N2}+\tilde u_{n+\frac N2},
\end{equation}
le deuxième terme étant par ailleurs nul si $n<0$ et le troisième étant nul si $n\geq 0$.
\end{cor}

\paragraph{D{é}monstration}
Appliquons le théorème \ref{Shannondiscret} à l'unique polynôme trigonométrique $P$ à $N$ coefficients qui a pour
échantillons les $u_k$. Alors par définition de la transformée de Fourier discrète, $\tilde u_n=c_n(P)$. On a donc pour
$\frac N4\leq n\leq \frac N4-1$,

$$\tilde v_n=\sum_{q\in \relative}c_{n+q\frac N2}(P)=\tilde
u_n+\tilde u_{n-\frac N2}+\tilde u_{n+\frac N2}.$$ Remarquons que si $n\geq 0$ cela donne $\tilde v_n=\tilde u_n+\tilde
u_{n-\frac N2}$, l'autre coefficient étant nul.  De même,  si $n< 0$, on obtient $\tilde v_n=\tilde u_n+\tilde
u_{n+\frac N2}$.

%\begin{eqnarray*}
%\hat u_k + \hat u_{\frac N2 + k} &=& \sum_{l=0}^{N\!-\!1} u_{l} \omega_{2N}^{kl} +
%u_{l}\omega_{2N}^{(\frac N2 + k)l} \\
%&=&\sum_{l=0}^{N\!-\!1} u_l \omega_{2N}^{kl} + u_l (-1)^l \omega_{2N}^{kl}\\
%&=&\sum_{n=0}^{\frac N2 -1} (u_{2n}) \omega_{N}^{k2n}=\sum_{n=0}^{\frac N2 -1}
%v_n \omega_{\frac N2}^{kn}=\hat v_k.
%\end{eqnarray*}
\enddemo
Cette proposition indique que le spectre du signal sous-échantillonné d'un facteur deux s'obtient en superposant à
lui-m\^eme le spectre du signal original avec un décalage de $\frac N2$. On dit qu'il y a repliement de spectre. Ainsi,
le spectre du signal sous-échantillonné contient généralement des informations non présentes dans le spectre du signal
de départ, ce qui se traduit sur le signal sous-échantillonné par l'apparition de structures périodiques n'ayant pas de
lien direct avec le contenu du signal. Ceci est particulièrement frappant dans le cas des signaux bi-dimensionnels,
pour lesquels on a un résultat identique à celui de la proposition \ref{sousech2}.

\begin{cor}
\label{sousech2-2D} Soit $(v_{k,l})=S_2((u_{k,l}))$ la sous-échantillonnée d'un facteur 2 d'une image discrète. (on
suppose que $\frac N2$ est pair). Alors $(\tilde v _{m,n})$, la transform{é}e de Fourier discr{è}te de $(v_{k,l})$,
s'écrit, pour $m,n=-\frac N4,...,\frac N4-1$,
\begin{equation}\label{aliasing2Dparsousechantillonnage}
\tilde v_{m,n}= \sum_{(\varepsilon_1,\varepsilon_2)\in\{0,1,-1\}}\tilde u_{m+\varepsilon_1\frac N2,
n+\varepsilon_2\frac N2}.
\end{equation}
\end{cor}
\dexer Montrer le corollaire \ref{sousech2-2D}.  Parmi les neuf termes de la somme de droite de
(\ref{aliasing2Dparsousechantillonnage}), il y a en général quatre termes non nuls: lesquels? Faire un dessin. \fexer

 Pour illustrer ce résultat, nous montrons deux exemples d'images sous-échantillonnées aux figures
\ref{aliassynth} (image synthétique) et \ref{mouche}, exemple o\`u l'apparition de structures périodiques est due à la
superposition, lors du sous-échantillonnage, des hautes fréquences de l'image. La manipulation numérique à faire pour
créer des effets de moiré dans une image est aussi simple que son interprétation est subtile : il suffit de prendre "un
point sur deux" de l'image. L'interprétation de l'opération se fait en Fourier : on a créé de basses fréquences
parasites en $c_n$ qui correspondent au "repliement" de hautes fréquences $c_{n+\frac N2}$. D'où l'apparition de
sinusoïdes qui n'ont rien à voir avec le signal original et qui créent des effets de moiré.


\begin{figure}
\label{aliassynth}
\begin{center}
\includegraphics[width=5cm]{Figures/alias.eps}
\end{center}
\caption[Repliement de spectre 1]{Exemple de repliement avec une image synthétique. En haut à gauche: image originale,
à droite son spectre. En bas à gauche: l'image sous-échantillonnée d'un facteur deux dans chaque direction, à droite le
spectre correspondant. Le spectre de l'image sous-échantillonnée est obtenu en périodisant le spectre de l'image
originale avec pour période le carré visible en surimpression. }
\end{figure}

\begin{figure}
\label{mouche}
\begin{center}
\subfigure[Image originale]{\includegraphics[width=5cm]{Figures/mouche.ps}}
%\hspace{1cm}
\subfigure[Sa TFD, non nulle en dehors du carré visible en surimpression]
{\includegraphics[width=5cm]{Figures/mouchetfd2c.eps}}

%\vspace{1cm}
\subfigure[Image sous-{é}chantillonn{é}e d'un facteur 2] {\includegraphics[width=5cm]{Figures/mouchesub.ps}}
%\hspace{1cm}
\subfigure[La TFD correspondante, sur laquelle il y a
repliement]{\includegraphics[width=5cm]{Figures/mouchesubtfd2.ps}}

\end{center}
\caption[Repliement de spectre 2]{Sous-{é}chantillonnage et repliement: le cas d'une image mal {é}chantillonn{é}e. Pour
les images (a), (b), (c), (d), le principe est le m\^eme que dans la figure \ref{aliassynth}, mais le détail de la
transformation du spectre est plus difficile à suivre ! les effets du repliement ({\it aliasing} en anglais) sont
particulièrement visibles sur les yeux de la mouche, image (c), qui présentent des oscillations à basse fréquence. Les
structures quasi-périodiques de l'image originale sont visibles sous formes de taches et de filaments sur le spectre
(b). Le repliement est d\^u à la présence de ces structures aux hautes fréquences: la TFD de l'image originale n'est
pas nulle en dehors du carré visible en surimpression figure (b). Ce type d'effet de moiré est visible dans de nombreux
DVD commerciaux.}
\end{figure}

\begin{figure}
\label{filtrage}
\begin{center}

\subfigure[Image obtenue par TFD inverse de b] {\includegraphics[width=4cm]{Figures/mouchefilt.ps}}
%\hspace{1cm}
\subfigure[Image obtenue en mettant à zéro les hautes fréquences de
\ref{mouche}-a]{\includegraphics[width=4cm]{Figures/mouchefilttfd.ps}}

\subfigure[Sous-échantillonnage: le repliement a disparu] {\includegraphics[width=4cm]{Figures/mouchefiltsub.ps}}
%\hspace{1cm}
\subfigure[TFD de c]{\includegraphics[width=4cm]{Figures/mouchefiltsubtfd.ps}}


\end{center}
\caption[Filtrage passe-bas d'une image]{Une solution possible pour éviter les effets de repliement illustrés sur la
figure \ref{mouche}. L'image (a) est l'image dont le spectre est le m\^eme que celui de l'image \ref{mouche}-(a) à
l'intérieur du carré, et est nul à l'extérieur (filtrage passe-bas). L'image (c) est l'image sous-échantillonnée
correspondante. On observe que l'effet de repliement a disparu.}
\end{figure}

Le résultat de la proposition \ref{sousech2} se généralise dans le cas d'un sous-échantillonnage d'ordre plus élevé,
comme le montre la proposition suivante:

\begin{pro}
\label{sousech} Soit $(v_k)=S_p((u_k))$ (on suppose que $N=pM$, pour un certain entier $M$). Alors $(\tilde v_k)$, la
transform{é}e de Fourier discr{è}te de $(v_k)$, s'écrit, pour $k=1...M-1$,
\begin{equation}
\tilde v_k=  \sum_{a=-p+1}^{p-1}\tilde u_{k+\frac{aN}{p}}.
\end{equation}
\end{pro}

\paragraph{Démonstration}
Appliquer de nouveau le théorème \ref{Shannondiscret} à l'unique polynôme trigonométrique à $N$ coefficients qui a pour
échantillons les $u_k$.  Ce polynôme vérifie $c_n(P)=\tilde u_n.$
\enddemo

On peut comparer les propositions \ref{sousech2} et \ref{sousech} et le corollaire \ref{sousech2-2D} aux théorèmes
\ref{Shannondiscret} et \ref{Shannondiscret2D}. Ces théorèmes  nous donnent  les conditions générales de Shannon et
Whittaker pour qu'un signal ou une image soient correctement échantillonnés : ces conditions sont que le spectre soit
borné (nombre fini $N$ ou $N^2$ de coefficients de Fourier) et que l'on dispose d'au moins $N$ (ou $N^2$) échantillons.
Les propositions \ref{sousech2} et \ref{sousech}
 sont plus pratiques : elles
ne donnent aucune hypothèse sur le signal ou l'image qui ont été échantillonnés et
 ont l'avantage de s'appliquer donc à un signal ou une image discrets, quelconques, qu'ils soient ou non
 issus d'un bon échantillonnage.

 \subsection{Manipulations propres et  simples et leur visualisation}

\begin{cor} {\bf (Représentation d'une image à partir de ses coefficients de Fourier)}
Les $\tilde u_{m,n}$ étant donnés pour $m,n=-\frac N2, \dots, \frac N2-1$, on désire représenter une image qui a ces coefficients de Fourier. Pour  cela, on fixe le domaine de l'image comme étant $[0,b]^2$ et on  pose donc
$$u(\bx)=\sum_{\bm\in\{-\frac N2,\dots,\frac N2-1\}^2} \tilde u_{\bm}e^{\frac{2i\pi \bm.\bx}{b}}.$$ On peut donc construire une image ayant des coefficients de Fourier donnés sur un domaine carré de largeur arbitraire $b$, et l'image reconstruite est donc $b$-périodique. Ses échantillons $u(\frac {\bk b}N)$, $\bk\in \{0,\dots, N-1\}^2$ redonnent par DFT les $\tilde u_\bm$, et inversement ces échantillons sont obtenus par DFT inverse à partir des $\tilde u_\bm$. Il  en résulte que, les coefficients de Fourier $\tilde u_\bm$ étant donnés l'on doit toujours fixer le domaine de reconstruction en donnant $b$. Les images reconstruites à partir de différents $b$ se déduisent l'une de l'autre par simple homothétie.
\end{cor}

 On part d'une image digitale donnée par ses échantillons $\bu=(u(k,l))$,  $k,l\in [0,N-1]^2$. Donc l'image  qui  apparaît  sur l'écran est de taille $N\times N$. Si on procède par exemple à un sous-échantillonnage de deux:  $\bv:=S_2\bu$ est une image avec quatre fois moins d'échantillons, et donc de dimensions deux fois plus petites sur l'écran. Pour visualiser $\bv$  et  $\bu$ à la même taille, on procède par {\it zero-padding}. Soit $v$ l'unique polynôme trigonométrique $N$-périodique et de degré $N/2$  passant  par les échantillons pairs de $\bu$. On a donc $\bv=S_2v$. Les coefficients $\tilde v_\bm$ de $v$ sont obtenus par la formule d'aliasage \ref{???} à partir de ceux de $u$. Ces coefficients sont en nombre $(\frac N2)^2$. Donc pour  comparer visuellement l'image sous-échantillonnée $\bv$ à $\bu$, il faut zoomer $\bv$ pour en montrer autant d'échantillons. Pour cela on procède par zéro-padding, en posant $\tilde v_{\bm}=0$ pour
  $\bm\in \{-\frac N2,\dots,\frac N2-1\}^2\setminus \{-\frac N4,\dots,\frac N4-1\}^2$. Alors on voit autant d'échantillons ($N^2$) pour $v$ que pour $u$ et on peut donc les comparer sur un écran à la même  taille.  On notera $ZP_M$ l'opération qui consiste à border de zéros un spectre de manière à passer d'une largeur $N$ à une largeur $M+N$, et donc à $(M+N)^2$ échantillons.

  Pour avoir  une notation computationnellement  complète,  il faut préciser donc quand on fait une DFT le nombre d'échantillons  $N^2$. On notera donc $DFT_N$ pour la DFT à $N$ échantillons si c'est un signal et $DFT_N$ également si c'est la DFT à $N^2$ échantillons d'une image. Avec ce formalisme simple et qui  correspond à ce qui est programmé, les opérations d'intérêt à mener sont:

  \begin{itemize}
  \item Partant de $N^2$ échantillons $\bu=u(\bk)$, faire $DFT_N^{-1}DFT(\bu)$  (simple vérification: on doit retrouver $\bu$);
  \item Faire $\bv=DFT_N^{-1}ZP_{\frac N2}DFT_{\frac N2}S_2\bu$:  met en évidence les effets néfastes de l'aliasage;
  \item Faire $\bw=DFT_{N}^{-1}ZP_{\frac N2}(DFT_{N}\bu)_{|\{-\frac N4,\dots\frac N4-1\}})$. En d'autres termes on fait les choses correctement: on annule les coefficients de Fourier plus grands que $N/4$, ce qui  autorise un sous-échantillonnage de 2.
  \item Comparer $\bv$ et $\bw$
  \end{itemize}

\subsection{L'équation de la chaleur}

Consider a digital image $\bu(\bk)$ with $\bk=(k,l)\in [0,N-1]^2$ and its DFT interpolate
\begin{equation}
u(\bx) =\!\!\! \sum_{\bm=(m,n)\in[-N/2,N/2-1]^2}  \tilde
u_\bm e^{\frac{2i\pi\bm.\bx}N},
\end{equation}
where $\tilde u_\bm$ with $\bm = (m,n)$ are the Discrete Fourier
Transform (DFT) coefficients of the $N^2$ samples $u(\bk)$.  The scale space convolves
$u$ with a Gaussian $G_\sigma(\bx)=\frac
1{2\pi\sigma^2}e^{-\frac{\bm^2}{2\sigma^2}}.$ A direct
 calculation shows that the result of the convolution of $u$
with $G_\sigma$ is
\begin{equation}
v(\bx):= (G_\sigma\ast u)(\bx)=\sum_{\bm\in[-\frac N2,\frac
N2-1]^2} \tilde {u}_{\bm} \hat
{G_\sigma}\left(\frac{2\bm\pi}N\right)e^{\frac{2i\pi\bm.\bx}N}
\end{equation}
 with
 $\hat G_\sigma(\xi)=e^{\frac{-\sigma^2\xi^2}2}$.
 Indeed, $$G_\sigma\ast e^{\frac{2i\pi\bm.\bx}N}=\int_{\R^2}G_\sigma(\by)e^{\frac{2i\pi\bm.(\bx-\by)}N}d\by=
 e^{\frac{2i\pi\bm.\bx}N}
 \int_{\R^2}G_\sigma(\by)e^{-\frac{2i\pi\bm.\by}N}d\by= e^{\frac{2i\pi\bm.\bx}N}\hat{G_\sigma}(\frac{2\pi\bm}{N}).$$
 The convolved digital image is
\begin{equation}
\label{eqn:DFT:v} v(\bk)=\sum_{\bm \in[-\frac N2,\frac N2-1]^2}
 \tilde v_\bm e^{\frac{2i\pi\bm.\bk}N}.
\end{equation}

Thus, the algorithm giving the {\it exact} solution $u(\sigma,\bx)$ of the heat equation starting from an initial trigonometric polynomial  $u_0(\bx)=u(0,\bx)$ will be a simple formula: Denote by  $\bu_0$ the $N^2$ samples  of $u_0$ and $\tilde u_\bm =DFT_N(\bu_0)$ its  DFT coefficients.  Call $\bu(\sigma)$ the $N^2$ samples of  $u(\sigma,\bx)$.  Then
$$\bu(\sigma)=DFT_N^{-1}((\hat{G_\sigma}(\frac{2\bm.\bx}N)\tilde u_\bm)_{\bm\in\{-\frac N2,\dots\frac N2-1\}}).$$
We denote  this discrete filter by $\bG_\sigma$: $\bu(\sigma)=\bG_\sigma \bu_0$.

Finally we can make the experiment checking on what  kind of Gaussian convolution is enough to make
a subsampling.
Trying $t=2, 4, 8$ and taking an initial digital image $\bu_0$ with $N^2$ samples, we just have to compare for several values  of $\bc$ the digital images
 $$\bu_1 =\bS_1 \bG_{t\bc}\bu \hbox{ and } \bu_2=\bS_{\frac 1t}\bS_t\bG_{t\bc}\bu_0.$$
 Notice that $\bS_{\frac 1t}$  is an oversampling, obtained therefore by zero padding.   More precisely,
 $$\bS_{\frac 1t}=: DFT_N^{-1}ZP_{N(1-\frac 1t)}.$$
 Thus it is not true in general that $\bS_{\frac 1t}S_t$ is the identity.
 The difference between $\bu_1$ and  $\bu_2$ can be visualized, but can also  be
 reliably computed as the root mean square error per pixel (RMSE). If the RMSE is lower than 1 or 2, the visual  difference is tiny, and if the RMSE is lower than, say, 0.1,  it becomes outright negligible.




%\paragraph{D{é}monstration}
%\begin{eqnarray*}
%\sum_{a=0}^{p-1}\hat u_{k+\frac{aN}{p}} &=& \sum_{a=0}^{p-1}\sum_{l=0}^{N\!-\!1} u_{l}
%\omega_N^{\left(k+\frac{aN}{p}\right)l} \\
% &=&\sum_{l=0}^{N\!-\!1} u_l \sum_{a=0}^{p-1}\exp\left(\frac{2i\pi
% l}{p}\right)\omega_{N}^{kl} \\
%&=&p\sum_{l'=0}^{M-1}u_{pl'}\omega_{M}^{kl'} = \hat v_k
%\end{eqnarray*}




 \subsection{La transform{é}e de Fourier rapide}
\label{TFR} \index{A}{transformée de Fourier!rapide}

Comme nous l'avons vu plus haut, le calcul des coefficients de Fourier $\tilde u_n$ revient {à} l'{é}valuation d'un
certain polyn{\^o}me aux racines $N$-i{è}mes de l'unit{é}. Dans le cas g{é}n{é}ral, l'{é}valuation classique (ex.
m{é}thode de H{\"o}rner) d'un polyn{\^o}me de degr{é} $N\!-\!1$ en un point prend $\mathcal{O}(N)$ op{é}rations. Donc
si l'on r{é}p{è}te cela pour les $N$ racines de l'unit{é} on devra effectuer $\mathcal{O}(N^2)$ op{é}rations.
L'algorithme de la Transform{é}e de Fourier Rapide (TFR) permet de r{é}soudre le probl{è}me en $\mathcal{O}(N\log N)$
op{é}rations. Appelons ``calcul d'ordre $N$'' l'{é}valuation d'un polyn{\^o}me de degr{é} $N\!-\!1$ aux racines
$N$-i{è}mes de l'unit{é}. Et soit $T(N)$ le nombre d'op{é}rations (additions et multiplications) demand{é}es par ce
calcul.

On se place dans le cas $N=2^n$ et soit un polyn{\^o}me
$$P(X)=\sum_{k=0}^{N\!-\!1}a_kX^k.$$
On pose
$$Q(X)=\sum_{k=0}^{\frac N2-1}a_{2k}X^k,$$
$$R(X)=\sum_{k=0}^{\frac N2-1}a_{2k+1}X^k.$$
Alors \begin{equation}\label{etoile}
P(\omega_N^k)=Q\left(\left(\omega_N^k\right)^2\right)+\omega_N^kR\left(\left(\omega_N^k\right )^2
\right).\end{equation}
Or, si $N$ est pair les $\left(\omega_N^k\right)^2$ sont exactement les racines d'ordre
$\frac{N}{2}$ de l'unit{é}. Il suffit donc d'{é}valuer les deux polyn{\^o}mes $Q$ et $R$ aux racines d'ordre
$\frac{N}{2}$ de l'unit{é} ce qui est un probl{è}me d'ordre $\frac{N}{2}$. On a donc, en tenant compte des additions et
multiplications demand{é}es par (\ref{etoile}),
$$T(N)=2T\left(\frac{N}{2}\right)+2N.$$
On en tire ais{é}ment $T(N)=\mathcal{O}(N\log(N))$. \cqfd

\dexer Démontrer la relation $T(N)=\mathcal{O}(N\log(N))$.

\fexer
\begin{rem}
\label{TFDMatlab} Les programmes usuels de calcul numérique ne calculent pas les coefficients $\tilde u_n$, mais les
coefficients $\hu_n$, définis par la formule suivante, pour $n=0,\;...,\;N\!-\!1$:
$$
 \hu_n= \left\{ \begin{array}{ll}
                                  \tilde u_n & \mbox{si $n=0,...,\frac{N}{2}-1$} \\
                                  \tilde u_{n-N} & \mbox{si $n=\frac{N}{2}, ..., N$} \\
                                 \end{array}\right.
                                 $$
\end{rem}


\subsection{L'utilisation de la transform{é}e de Fourier discr{è}te pour
d{é}finir zoom, translations et rotations des images} \label{zoom}
\paragraph{Le zoom}

Nous pr{é}sentons une m{é}thode d'interpolation reposant sur une extension de la TFD d'un signal ou d'une image. Nous
d{é}taillons la m{é}thode, dite du ``prolongement par des zéros (``zéro-padding''), en une dimension, le principe se
g{é}n{é}ralisant sans mal pour une image. Comme pr{é}c{é}demment, consid{è}rons des échantillons $u_k$, $k$ variant de
0 {à} $N\!-\!1$, et $\tilde{u}_n = \frac{1}{N}\sum_{k=0}^{N-1}u_{k}\omega_N^{-kn}$. On suppose que N
est pair et que l'on veut zoomer d'un facteur 2, c'est {à} dire que l'on veut construire un signal de taille deux fois
plus grande (avec deux fois plus d'échantillons) que le signal de d{é}part. On d{é}finit un nouveau signal $v$, de
taille $2N$ comme {é}tant la TFD inverse de $\tilde {v}$, donn{é} par

\begin{equation}
\label{TFDzoom}
 \tilde{v}_n =  \tilde u_n \; \hbox { si } -\frac N2\leq n\leq \frac N2-1,\;\; \tilde v_n=0
\hbox{
 si  } n\in[-N, -\frac N2-1]\cup  [\frac N2, N-1].
\end{equation}

\begin{pro} \index{A}{zoom discret par Fourier} \label{zeropadding}
Le signal $v$ dont la TFD est donn{é}e par la formule (\ref{TFDzoom}) v{é}rifie  $v_{2k}=u_k$, pour $k=0, ...,
N\!-\!1$.
\end{pro}

\paragraph{D{é}monstration}
On a $$v_{2k}=\sum_{-N}^{N\!-\!1}\tilde v_n\omega_{2N}^{2nk}=\sum_{-\frac N2}^{\frac N2-1}\tilde u_n\omega_N^{nk}
=u_k.$$ En effet, $\omega_{2N}^{2nk}=\omega_N^{nk}.$ \cqfd
 \begin{rem} Ce résultat est  évident sans démonstration : en effet, on peut considérer
l'unique
 polynôme trigonométrique de degré $\frac N2$ passant par les échantillons $u_k$. Les
échantillons
 $v_k$ s'interprètent immédiatement comme des échantillons de ce même polynôme.
 \end{rem}

\begin{rem}
On remarquera que les nouveaux  échantillons obtenus par cette méthode peuvent \^etre complexes, m\^eme lorsque le
signal original est réel (ceci étant d\^u au terme d'aliasage $u_{-\frac N 2}$). Aussi, on ne retient que  la partie
réelle à des effets de visualisation.
\end{rem}
\begin{figure}
\begin{center}

\psfrag{a}{a} \psfrag{b}{b}
\includegraphics[width=10cm]{Figures/aerozoom2.eps}

\includegraphics[width=10cm]{Figures/aerodup.ps}
\end{center}
\caption[Zoom par prolongement du spectre 1]{Zoom sur une partie de l'image \ref{image-TFD}. Haut: zoom par TFD
(zéro-padding), bas: zoom par r{é}plication des pixels. Le zoom par TFD est obtenu en prolongeant par des zéros le
spectre de l'image initiale. Celui par réplication des pixels en rempla\c{c}ant chaque pixel par quatre pixels de la
m\^eme valeur. Remarquons tout d'abord la plus grande régularité du zoom par TFD, qui supprime les effets de ``blocs''
très visibles sur le zoom par réplication. En contrepartie, le phénomène de Gibbs (voir paragraphe \ref{Gibbs}) est
très visible sur le zoom par TFD, puisque l'on a mis à zéro brutalement des coefficients de la TFD. Ce phénomène est
particulièrement visible le long des frontières de l'image, qui correspondent à des discontinuités puisque l'image est
périodisée (par exemple zone a). On observe aussi le phénomène de Gibbs le long des contours des objets (par exemple
zone b).} \label{fftzoom}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=8cm]{Figures/mains.ps}
\includegraphics[width=8cm]{Figures/dalles2.ps}
\end{center}
\caption[Zoom par prolongement du spectre 2]{détails apr{è}s zoom, {à} gauche par TFD, {à} droite par r{é}plication des
pixels.} \label{details}
\end{figure}
La m{é}thode se g{é}n{é}ralise aux cas des images. Nous consid{é}rons une image num{é}rique $(u_{k,l})$, et nous
d{é}finissons une image zoom{é}e $\left(v_{i,j}\right)_{i,j=0,\;...,\; 2N\!-\!1}$ comme {é}tant la transform{é}e de
Fourier discr{è}te inverse de $\tilde v_{i,j}$ définie pour $i,j=-N,\;...,\;N\!-\!1$ par

\begin{equation}
\label{TFDzoom2}
 \tilde{v}_{m,n} = \tilde u_{m,n} \hbox{  si  } -\frac N2\leq m, n\leq \frac N2-1, \;\; \tilde
v_{m,n}=0
 \hbox {  sinon}.
\end{equation}

\dexer  Ecrire l'équivalent de la proposition \ref{zeropadding} et la montrer en dimension 2. \fexer
La figure
\ref{fftzoom} montre la partie réelle d'une partie de l'image \ref{image-TFD} zoom{é}e par zéro-padding, et la compare
celle obtenue par r{é}plication des pixels (chaque pixel est remplacé par quatre pixels de la m\^eme valeur). On
remarque que le zoom par TFD produit une image bien plus r{é}guli{è}re, et {é}vite l'effet ``marche d'escalier''
visible sur l'image zoom{é}e par r{é}plication. La figure \ref{details} illustre ce point sur un d{é}tail. Une autre
remarque concerne l'effet de Gibbs (cf. paragraphe \ref{Gibbs}). Ce ph{é}nom{è}ne produit des rebonds le long de la
frontière du domaine de l'image. En effet, et comme nous l'avons déjà mentionné, le calcul des coefficients de Fourier
de l'image (dont les coefficients de la TFD sont une approximation) suppose l'image p{é}riodique, ce qui fait
apparaître des discontinuit{é}s le long des fronti{è}res de son domaine de d{é}finition. Le phénomène de Gibbs est
également visible le long des discontinuités dans l'image, les contours. Le phénomène est mis en évidence sur la figure
\ref{fftzoom}. Expliquons pourquoi le ph{é}nom{è}ne appara{\^\i}t dans le cas du zoom: une nouvelle image $v_{k,l}$ de
taille $2N\times2N$ est obtenue en utilisant les valeurs prises par le polyn{\^o}me $P(x)$ entre les points dont on
dispose au d{é}part. Cette utilisation de $P$ fait appara{\^\i}tre des oscillations du polynôme trigonométrique qui
{é}taient invisibles sur les échantillons de d{é}part  $(u_{k,l})$. Comme nous l'avons d{é}j{à} {é}voqu{é}, les
oscillations aux fronti{è}res du domaine de l'image peuvent {\^e}tre supprim{é}es par utilisation de la transform{é}e
en cosinus. En revanche, le probl{è}me subsistera le long des discontinuit{é}s pr{é}sentes {à} l'int{é}rieur de
l'image, à moins de faire un ``filtre anti-aliasage''. Ce filtre consiste à atténuer les fréquences près du bord,  pour
que leur transition à  zéro soit moins brutale.

%\begin{figure}
%\begin{center}
%\includegraphics[width=9cm]{Figures/detail.ps}
%\end{center}
%\caption{Zoom sur une partie de l'image utilisant sa TFD (voir le
%texte). On distingue nettement les oscillations dues au
%ph{é}nom{è}ne de Gibbs le long des contours, ainsi que le long
%de la fronti{è}re inf{é}rieure de l'image.} \label{figgibbs}
%\end{figure}


\paragraph{La translation}

La m{é}thode pr{é}sent{é}e au paragraphe pr{é}c{é}dent permet de d{é}finir une translation d'une quantit{é} 1/2 (ou
$a/(2N)$ pour revenir {à} notre d{é}finition premi{è}re du signal $u$), en ne gardant que les points d'indice impair du
signal zoom{é} $v$. Plus g{é}n{é}ralement, nous pouvons d{é}finir une translation d'un signal d'une quantit{é}
$0<\alpha<1$. Comme d'habitude, l'opération de translation sur la fonction $u$ dont nous connaissons les échantillons
$u_k$ se fait sous l'hypothèse que celle-ci est un polynôme trigonométrique. En d'autres termes, on translate le
polynôme d'interpolation, la "vraie" fonction $u$ étant inconnue en dehors des échantillons. Le polynôme
d'interpolation est
$$P(x)=\sum_{-\frac N2}^{\frac N 2-1}\tilde u_n e^{\frac {2i\pi nx}a}.$$
En translatant de $\alpha$, on obtient
$$\tau_\alpha P(x)=P(x-\alpha)=\sum_{-\frac N2}^{\frac N2-1} \tilde u_ne^{-\frac{2i\pi n\alpha}a}
 e^{\frac{2i\pi nx}a}.$$ On a donc :


\begin{pro}\label{translation}% \index{translat\'{e}e!discr\`{e}te par Fourier}
 La TFD $(\tilde v_n)$ de $P(x-\alpha)$ s'obtient à partir de la TFD  de $P(x)$, $\tilde u_n$, par
$$\tilde v_n=\tilde u_n e^{-\frac{2i\pi n\alpha}a}.$$
\end{pro}


Cette m{é}thode de translation se g{é}n{é}ralise sans mal au cas des images, en remarquant qu'une translation {à} deux
dimensions peut se d{é}composer en deux translations, une selon les lignes et une selon les colonnes.

\dexer Donner les formules de la translation par TFD en dimension 2. \fexer


\paragraph{La rotation}
Décrivons maintenant une m{é}thode pour impl{é}menter une rotation discr{è}te, due {à} L. Yaroslavsky. En bref, cette
méthode réduit une rotation à des translations en ligne ou  en colonne  de l'image. Commen\c{c}ons par remarquer que
\begin{equation}\small
\label{decomprot} R(-\theta):=\left(\begin{array}{cc}
           \!\! \cos(\theta) &\!\!  \sin(\theta)\!\!  \\
           \!\!  -\sin(\theta) & \!\! \cos(\theta)\!\!
          \end{array} \right)
 \!\! =\!\!  \left(\begin{array}{cc}
            1 & \tan(\frac{\theta}{2}) \\
            0 & 1
          \end{array} \right)\!\!
   \left(\begin{array}{cc}
            \!\! 1 &\!\!  0 \\
            \!\! -\sin(\theta) &\!\!  1
          \end{array} \right)\!\!
   \left(\begin{array}{cc}
            1 & \tan(\frac{\theta}{2}) \\
            0 & 1
          \end{array} \right)\!\! :=T(\theta)S(\theta)T(\theta)
\end{equation}

(sauf si $\theta = \pi$ auquel cas il suffit de retourner l'image).

Une rotation d'angle $\theta$ de l'image discrète $u(i,j)$ consiste à calculer $u(R(-\theta)(i,j))$ que l'on notera
$(R(-\theta)u)(i,j)$. Mais on  a $R(-\theta)u= T(\theta)S(\theta)T(\theta)u$. Donc il  suffit d'expliquer comment
calculer $T(\theta) u$ et $S(\theta)u$. Or ces deux opérations ont la même structure, à savoir une translation ligne
par ligne ou une translation colonne par colonne. Traitons par exemple le cas de $T(\theta)$. On a
$(T(\theta)u)(i,j)=u(i+j\tan(\frac\theta 2),j)$. Donc partant de la matrice $u_{i,j}$, on translate sa premi{è}re ligne
de $\tan(\frac{\theta}{2})$, la deuxi{è}me de $2\tan(\frac{\theta}{2})$, etc. Appliquer $S(\theta)$  revient à faire
une op{é}ration similaire sur les colonnes. Enfin on réapplique $T(\theta)$ et on fait donc à nouveau une translation
 sur les lignes. Or comme on vient de le voir ces translations ligne à ligne ou colonne à colonne se font en temps $N\log N$  en utilisant
   la TFD {à} une dimension.
\begin{figure}
\begin{center}

\includegraphics[width=8cm]{Figures/aerot.ps}

\end{center}
\caption[Rotation par Transformée de Fourier 1]{Rotation de $\pi/4$ par TFD. La rotation est implémentée en remarquant
qu'elle peut se décomposer en trois transformations consistant en des translations selon les lignes ou les colonnes de
l'image (formule \ref{decomprot}). Chacune de ces transformations est ensuite effectuée gr\^ace à une TFD sur la ligne
ou colonne considérée, en utilisant la méthode présentée au paragraphe précédent.} \label{fftrot}
\end{figure}
\begin{figure}
\begin{center}

\includegraphics[width=10cm]{Figures/bilin.ps}
\includegraphics[width=10cm]{Figures/fftrot.ps}

\end{center}
\caption[Rotation par Transformée de Fourier 2]{Bas: apr{è}s douze rotations successives de $\pi/4$ par TFD; haut:
m{\^e}me exp{é}rience en utilisant une interpolation bilin{é}aire (la valeur en un nouveau point $(x,y)$ est obtenue
par combinaisons lin{é}aires
 des valeurs aux quatre points {à} coordonn{é}es enti{è}res de l'image
originale les plus proches de $(x,y)$).} \label{rotiter}
\end{figure}

La figure \ref{fftrot} montre une image apr{è}s une rotation de $\pi/4$ par la m{é}thode d{é}crite ci-dessus. Puis,
pour illustrer la stabilit{é} de la m{é}thode, nous montrons figure \ref{rotiter} le r{é}sultat de l'application
successive de douze rotations de $\pi/4$, et, {à} titre de comparaison, le r{é}sultat de ces douze rotations
successives impl{é}ment{é}s par interpolation bilin{é}aire (les valeurs aux nouveaux points sont des combinaisons
lin{é}aires
 des quatre points {à} coordonn{é}es enti{è}res les plus proches). Cette figure
illustre clairement la sup{é}riorit{é} de la m{é}thode par FFT dans le cas de rotations multiples.




\begin{rem}
Du fait que l'on manipule des fonctions périodiques, une translation conduirait à faire sortir une partie de l'image
par un bord pour la faire entrer par l'autre. Ce qui conduirait à l'apparition, sur les bords de l'image d'un certain
nombre de détails qui sont en fait mal placés.
%Les schémas suivants le montrent bien~:
%\begin{center}
%\begin{tabular}{cc}
%\includegraphics{schemas/rotation.0.eps}&
%\includegraphics{schemas/rotation.1.eps}\\
%{\bf image originale}&
%{\bf prmière translation}\\
%\includegraphics{schemas/rotation.2.eps}&
%\includegraphics{schemas/rotation.3.eps}\\
%{\bf deuxième translation}&
%{\bf troisième translation}
%\end{tabular}
%\end{center}
On  se débarrasse facilement de ce problème en insérant l'image dans un cadre deux fois plus grand\dots
\end{rem}
\begin{rem}
La méthode de rotation n'est pas parfaite. En effet, l'image $u$ continue associée à $u(i,j)$  est dans l'interpolation
  supposée implicitement $N$-périodique, ce qui revient à dire qu'elle est de la forme (pour une image
carrée)
$$u(x,y)=\sum_{k,l=0}^{N\!-\!1}c_{i,j}e^{ 2i\frac{\pi}{N}(kx+ly) }.$$
Mais, si on lui applique une "translation" suivant l'axe des $x$ de valeur $\lambda y$, la formule devient
$$u_1(x,y)=\sum_{k,l=0}^{N\!-\!1}c_{i,j}e^{2i\frac{\pi}{N}(kx+(l-\lambda k)y)}.$$
La fonction $u_1$ n'est pas (pour $\lambda\notin\ZZ$) $N$-périodique en $y$. Or, après la première translation on ne
dispose plus que des échantillons du signal $u_1$ sur une grille carrée $N\times N$. D'après la théorie de Shannon un
tel ensemble de données ne permet pas de capturer toute l'information sur $u_1$ (à la seconde étape on effectue des
translations suivant $y$ qui est justement l'axe qui pose problème). On rencontre encore ce problème à la troisième
translation. Le seul moyen d'avoir une rotation exacte serait d'évaluer $u$ aux points de l'image de
$[0,\;N\!-\!1]\times[0,\;N\!-\!1]$ par une rotation d'angle $-\theta$, mais cette méthode est en $N^4$ ce qui la rend
inopérante\dots
\end{rem}

\subsection{Importances relatives de la phase et du module de la TFD pour une image}

Nous nous int{é}ressons {à} la pertinence visuelle des caract{é}ristiques de la transform{é}e de Fourier discr{è}te
dans le cas des images, et plus particuli{è}rement {à} la phase et au module de la TFD, au moyen de deux exemples. Tout
d'abord nous montrons, figure \ref{phase}, deux images A et B, ainsi que les images obtenues en échangeant les phases
de leurs TFD. Nous remarquons gr\^ace {à} cette exp{é}rience qu'une part tr{è}s importante de l'information
g{é}om{é}trique d'une image est contenue dans la phase de sa TFD. Rappelons que si l'on translate une fonction, les
coefficients de sa s{é}rie de Fourier sont multipli{é}s par des exponentielles complexes de module 1, et que par
cons{é}quent la phase de la TFD contient en en sens des informations sur le placement des constituants de l'image.

Dans la figure \ref{module}, nous montrons deux images de textures, qui visuellement semblent invariantes par
translation, ainsi que les deux images obtenues {à} partir de ces textures en ne conservant que le module de leur TFD,
et en tirant au hasard les phases (selon une loi uniforme). On voit cette fois que le module de la TFD contient
l'information. Cette propri{é}t{é} est caract{é}ristique des textures homog{è}nes du type pr{é}sent{é} figure
\ref{module}, et l'on peut m{\^e}me donner une d{é}finition des ``microtextures'' comme images caract{é}ris{é}es
uniquement par le module de leur transform{é}e de Fourier.
%Pour conclure, nous illustrons figure
%\ref{synth} une possibilit{é} de synth{è}se de textures par
%modification al{é}atoire des phases de la TFD de l'image binaire
%(c'est {à} dire en noir et blanc) d'une forme simple. L'image de
%droite est simplement obtenue en calculant la TFD de l'image de
%droite, en rempla{\c c}ant ses phases par des nombres
%al{é}atoires uniform{é}ment r{é}partis entre 0 et $2\pi$,
%puis en calculant la TFD inverse de l'image obtenue. Se rappelant
%de l'effet d'une translation sur les phases de la TFD, nous
%pouvons consid{é}rer que la méthode consiste en une
%superposition de la forme de d{é}part centr{é}e en des
%positions al{é}atoires, suivie d'un seuillage de l'image.

\begin{figure}
\begin{center}
\subfigure[Image A]{\includegraphics[width=6cm]{Figures/house.ps}}
%\hspace{1cm}
\subfigure[Image A]{\includegraphics[width=6cm]{Figures/lena.ps}}

%\vspace{1cm}
\subfigure[Module de la TF de A et phase de B]{\includegraphics[width=6cm]{Figures/lenmais.ps}}
%\hspace{1cm}
\subfigure[Module de la TF de B et phase de A]{\includegraphics[width=6cm]{Figures/maislen.ps}}
\end{center}
\caption[Phase de la Transformée de Fourier Discrète et information géométrique]{Haut: les deux images de d{é}part;
bas: les deux images apr{è}s {é}change des phases de leurs TFD. L'information g{é}om{é}trique est contenue dans la
phase ! Les formes sont principalement codées dans l'argument des coefficients de Fourier de l'image. Bien que les
images (a) et (c) d'une part, et (b) et (d) d'autre part, aient des modules complétement différents, on y distingue les
m\^emes formes géométriques. Remarquons également que les directions horizontales et verticales très présentes sur
l'image (a) apparaissent sous forme de texture dans l'image (c). Cette remarque est précisée par l'expérience de la
figure \ref{module}.} \label{phase}
\end{figure}
\begin{figure}
\begin{center}

\includegraphics[width=10cm]{Figures/text.ps}

\end{center}
\caption[Module de la Transformée de Fourier Discrète et texture]{Haut: deux images de textures; bas: les deux images
apr{è}s remplacement des phases de leurs TFD par des phases al{é}atoires. Une information essentielle sur la texture
est donc présente dans le module des coefficients de Fourier de l'image. Pour la texture de gauche, il semble que la
plupart de l'information soit contenue dans le module de la TFD. A droite, quelques aspects de la texture sont perdus.
Nous renvoyons le lecteur intéressé à un article (en anglais) sur une méthode de synthèse de texture par modification
de la phase: [Van Wijk]} \label{module}
\end{figure}


%\begin{figure}
%\begin{center}

%\includegraphics[width=7cm]{Figures/spot1.ps}

%\end{center}
%\caption{Gauche: images de d{é}part binaires; droites: textures
%obtenues en remplaçant les phases des TFD des images de d{é}part
%par des valeurs al{é}atoires suivant une loi uniforme}
%\label{synth}
%\end{figure}%

%\begin{figure}
%\begin{center}
%
%\includegraphics[width=7cm]{Figures/spot2.ps}
%
%\end{center}
%\caption{M{\^e}me principe qu'{à} la figure \ref{synth}}
%\label{synth2}
%\end{figure}
%
%
%\section{Lien avec la théorie de Shannon}
%\begin{thf} \index{A}{Shannon!théorème pour les polynômes
%trigonométriques} {\rm (de Shannon pour les polyn\^omes trigonométriques)} Soit un signal trigonométrique
%$$f(t)=\sum_{n=-N}^Nc_ne^{2i\pi\lambda_nt}.$$
%On a encore la formule de Shannon
%$$\forall a\in\left]0,\frac{1}{2\lambda_c}\right[,~\forall t\in\RR,
%~f(t)=\sum_{n=-\infty}^{+\infty}f(na)\frac{\sin\frac{\pi}{a}(t-na)} {\frac{\pi}{a}(t-na)},$$ avec
%$\lambda_c=\max\left\{|\lambda_n|\right\}$. La convergence est ponctuelle.
%\end{thf}
%\begin{rem}
%Ce théorème complète le théorème de Shannon pour un signal qui est ni périodique ni dans $L^2$.
%%la nécessité d'un tel théorème (en plus du th\'orème de Shannon
%%classique) vient du fait que la transformée de
%%Fourier d'une fonction périodique, qui est un peigne de Dirac, n'est en
%%général pas dans $\mbox{{\rm L}}_2$.
%\end{rem}
%{\bf Démonstration}
%
%il suffit de démontrer le résultat dans le cas d'une seule onde. Soit donc
%$$f(t)=e^{2i\pi\lambda t},~~~~\lambda\in\RR.$$
%Soit $g$ périodique de période $\frac{1}{a}$ et égale à $f$ sur $\left(-\frac{1}{2a},\frac{1}{2a}\right)$. les
%coefficients de Fourier de $f$
%  sont
%$$c_n=\frac{a\sin\frac{\pi}{a}(\lambda-na)}{\pi(\lambda-na)}.$$
%Donc
%$$g(t)=\sum_{-\infty}^{+\infty}\frac{\sin\frac{\pi}{a}(\lambda-na)}
%{\frac{\pi}{a}(\lambda-na)}e^{2i\pi
%  nat}.$$
%
%Comme $f$ est $C^1$ sur $\left]-\frac{1}{2a},\frac{1}{2a}\right[$, Cette égalité est ponctuelle pour
%$t\in\left]-\frac{1}{2a},\frac{1}{2a}\right[$ (principe de localisation). D'o\`u
%$$\forall\lambda\in\RR,~~e^{2i\pi\lambda t}=
%\sum_{n=-\infty}^{+\infty}e^{2i\pi
%  nat}\frac{\sin\frac{\pi}{a}(\lambda-na)}
%  {\frac{\pi}{a}(\lambda-na)},~~|t|<\frac{1}{2a}. $$
%
%En intervertissant $\lambda$ et $t$, on obtient
%$$\forall t\in\RR,~~e^{2i\pi\lambda t}=
%\sum_{n=-\infty}^{+\infty}e^{2i\pi
%  na\lambda}\frac{\sin\frac{\pi}{a}(t-na)}{\frac{\pi}{a}(t-na)},~~
%  |\lambda|<\frac{1}{2a}. $$
