Transformée
de Fourier Rapide
FFT = Fast Fourier Transform
SPECTRE
(suggest update)
Ce programme en langage Turbo-Pascal permet de visualiser le spectre d'une
fonction temporelle ou spatiale.
PROGRAM Spectre; (* Transformee integrale de Fourier *)
(* calcule et trace la F.F.T. d'une
suite regulierement echantillonee ayant 2^N points *)
USES crt, graph; (* USES pour le
Turbo-Pascal *)
CONST
kDIM = 2048; (* kDIM = 2^11 =
taille maximale du vecteur *)
kTAB = 8200; (* il faut kTAB
>= kDIM+2 *)
DoublePi = 6.28318530717959; (*
Do Not Modify *)
kpixelXmax = 639; (* il y a 640
pixels horizontaux (//X) en VGA *)
kpixelYmax = 479; (* il y a 480
pixels verticaux (//Y) en VGA *)
TYPE
VAR
(* Exercice: determiner le pas du
spectre (dw) en fonction de Pas et de kDIM *)
FUNCTION
F1(x:real):real; (* x peut etre le temps t *)
BEGIN F1:= sin(x) + random/2; (*
random/2 varie entre 0 et 1/2 *)
END; (* sinus avec du bruit aleatoire
*)
FUNCTION
F2(t:real):real;
BEGIN IF (x>DoublePi) and
(x<2*DoublePi) THEN F2:= sin(t) ELSE F2:=0;
END; (* une seule periode de sinus
entre 2Pi et 4Pi *)
FUNCTION
F3(t:real):real;
BEGIN F3:= sin(t) + sin(1.1*t);
END; (* battements de 2 frequences
voisines *)
FUNCTION
F4(x:real):real;
BEGIN IF ODD(TRUNC(x/DoublePi)) THEN
F4:=-1 ELSE F4:=1;
END; (* fonction creneau symetrique *)
FUNCTION
F5(x:real):real; (* Exercice: essayer plusieurs fonctions
*)
BEGIN F5:= sin(x) - 0.3*sin(2.3*x);
END; (* sinus(x) avec du bruit
periodique *)
PROCEDURE ENTREE; (*
but: determiner la Fonction et son Pas *)
(* parametre d'entree: aucun *)
(* parametres de sortie: Choix, Pas,
VEC *)
(* Exercice: comment remplir VEC[i] en
lisant un fichier *)
VAR i:integer;
BEGIN
ClrScr; (* clear screen *) writeln;
writeln('Le site
"http://www.geocities.com/facultedegenie"
propose ...');
writeln('le programme SPECTRE
(version 2a) ...');
writeln('de Transformee de Fourier
rapide (F.F.T.) Fast Fourier Transform'); writeln;
writeln('Le premier ecran
representera le debut (640points de 2048) de la
fonction et son spectre');
writeln('Le 2e ecran sera le debut
du "nouveau" spectre et de la fonction
filtree (deduite)');
writeln; writeln('N.B. Le signe de
la FFT n''est pas significatif !');
writeln; write('- Choisir une
fonction de 1 à 5 (ou 0 pour quitter) :');
readln(Choix);
Pas:= DoublePi/50;
CASE Choix OF
1: begin Pas:= DoublePi/50; for
i:= 0 to kDIM do VEC[i]:= F1(i*Pas);
end;
2: begin Pas:= DoublePi/100;
for i:= 0 to kDIM do VEC[i]:= F2(i*Pas);
end;
3: begin Pas:= DoublePi/20; for
i:= 0 to kDIM do VEC[i]:= F3(i*Pas);
end;
4: begin Pas:= DoublePi/50; for
i:= 0 to kDIM do VEC[i]:= F4(i*Pas);
end;
5: begin Pas:= DoublePi/50; for
i:= 0 to kDIM do VEC[i]:= F5(i*Pas);
end;
ELSE Choix:=0;
end; (* of case*)
END;
PROCEDURE FILTRER;
(* but: le filtre modifie le spectre
VEC[i] obtenu apres FFT *)
(* parametre d'entree: VEC spectre issu
de FFT *)
(* parametre de sortie: nouveau spectre
VEC *)
(* Exercice: changer le filtre en
modifiant cette procedure *)
VAR i:integer;
BEGIN (* exemple ci-dessous: filtre
passe bas *)
for i:=imax+11 to kDIM do
begin VEC[i]:=VEC[i]/10;
end;
END;
(******* NE PAS MODIFIER la procedure
FOURIER1 ci-dessous *******)
PROCEDURE FOURIER1(VAR
data:tVECTEUR; n1,isign:integer);
(* parametres d'entree: data =
tableau pour la procedure FOURIER_RAPIDE *)
(* n1 = la demi dimension de ce
tableau *)
(* isign = +1 pour la TF
directe, ou -1 pour l'inverse *)
(* parametres de sortie: data = tableau
pour FOURIER_RAPIDE *)
(* Attention: tVECTEUR doit avoir au
moins 2*n+2 elements *)
(* L'utilisateur n'a pas besoin d'en
savoir plus pour le moment *)
VAR
i1, J1, n, mmax, m, J, istep, i
:integer;
wtemp, wr, wpr, wpi, wi, theta,
tempr, tempi :real;
BEGIN
n:= 2*n1; J:= 1;
FOR i1:= 1 to n1 DO
BEGIN
i:= 2*i1-1;
IF J>i THEN
begin (* permutations *)
tempr:= data[J]; data[J]:=
data[i]; data[i]:= tempr;
tempi:= data[J+1];
data[J+1]:= data[i+1]; data[i+1]:= tempi;
end;
m:= n DIV 2;
WHILE (m >= 2) AND (m <
J) DO
begin
J:= J-m; m:= m DIV 2;
end;
J:= J+m;
END;
mmax:= 2;
WHILE n>mmax DO
begin
istep:= 2*mmax;
theta:= DoublePi/(isign*mmax);
wpr:= -2.0*sqr(sin(0.5*theta));
wpi:= sin(theta); wr:= 1.0;
wi:= 0.0;
FOR i1:= 1 to (mmax DIV 2) DO
BEGIN
m:= 2*i1-1;
FOR J1:= 0 to ((n-m) DIV
istep) DO
begin
i:= m+J1*istep; J:=
i+mmax;
tempr:=
wr*data[J]-wi*data[J+1];
tempi:=
wr*data[J+1]+wi*data[J];
data[J]:=
data[i]-tempr; data[J+1]:=
data[i+1]-tempi;
data[i]:=
data[i]+tempr; data[i+1]:=
data[i+1]+tempi;
end; (* of for J1...*)
wtemp:= wr; wr:=
wr*wpr-wi*wpi+wr; wi:= wi*wpr+wtemp*wpi+wi;
END; (* of for i1...*)
mmax:= istep;
end; (* of while...*)
END;
(******* NE PAS MODIFIER la procedure
FOURIER_RAPIDE ci-dessous *******)
PROCEDURE
FOURIER_RAPIDE(VAR data:tVECTEUR; n,isign:integer);
(* Parametres d'entree: data =
tableau pour FOURIER_RAPIDE *)
(* n = la demi dimension utile
de ce tableau *)
(* isign = +1 pour la TF
directe, ou -1 pour l'inverse *)
(* Parametres de sortie: data = tableau
pour FOURIER_RAPIDE *)
(* Attention: FOURIER_RAPIDE utilise
FOURIER1 *)
(* L'utilisateur n'a pas besoin d'en
savoir plus pour le moment *)
VAR
wr, wi, wpr, wpi, wtemp, theta,
c1, c2, h1r, h1i, h2r, h2i, wrs, wis :real;
i, i1, i2, i3, i4 :integer;
BEGIN
theta:= DoublePi/(2.0*n);
wr:= 1.0; wi:= 0.0; c1:= 0.5;
IF isign>=0 THEN
BEGIN
isign:= 1; (* FFT directe *)
c2:= -0.5;
FOURIER1(data,n,1);
data[2*n+1]:= data[1];
data[2*n+2]:= data[2]
END
ELSE
BEGIN
isign:= -1; (* FFT inverse *)
c2:= 0.5; theta:= -theta;
data[2*n+1]:= data[2];
data[2*n+2]:= 0.0; data[2]:= 0.0
END;
wpr:= -2.0*sqr(sin(0.5*theta));
wpi:= sin(theta);
FOR i:=1 to (n DIV 2)+1 DO
BEGIN
i1:= i+i-1; i2:= i1+1; i3:=
n+n+3-i2; i4:= i3+1;
wrs:= wr; wis:= wi;
h1r:= c1*(data[i1]+data[i3]);
h1i:= c1*(data[i2]-data[i4]);
h2r:= -c2*(data[i2]+data[i4]);
h2i:= c2*(data[i1]-data[i3]);
data[i1]:= h1r+wrs*h2r-wis*h2i;
data[i2]:= h1i+wrs*h2i+wis*h2r;
data[i3]:= h1r-wrs*h2r+wis*h2i;
data[i4]:= -h1i+wrs*h2i+wis*h2r;
wtemp:= wr; wr:=
wr*wpr-wi*wpi+wr; wi:= wi*wpr+wtemp*wpi+wi
END; (* of for i...*)
IF isign>=0 THEN
data[2]:=data[2*n+1]
ELSE FOURIER1(data,n,-1);
END;
PROCEDURE
TRACER1(pYmin,pYmax:integer);
(* but: tracer le tableau VEC[i] sur
l'ecran graphique *)
(* parametres d'entree: les coordonnees
du cadre en pixels *)
(* parametres de sortie: affichage *)
(* attention: le mode graphique doit
etre actif (initGraf) *)
(* N.B. Au lieu
de tracer ici VEC[i] vous pouvez l'enregister dans un
fichier (ecrire une procedure) puis le tracer en
utilisant: Excel, MatCad, MatLab, Mathematica, ou ... *)
VAR
i, ix, iy, jx, jy :integer;
maxA :real;
BEGIN
maxA:= VEC[0]; imax:= 0; (*
initialisations *)
for i:=1 to kDIM do
begin
if abs(VEC[i]) > maxA then
begin maxA:= abs(VEC[i]);
imax:=i;
end;
end; (* fin de la recherche de maxA
en valeur absolue *)
jy:= round( (pYmin+pYmax)/2 + (
(pYmax-pYmin)/(2*maxA) )*VEC[0] );
jx:=0;
for i:=1 to kpixelXmax do
begin
ix:= i; (* le + simple *)
iy:= round( (pYmin+pYmax)/2 + (
(pYmax-pYmin)/(2*maxA) )*VEC[i] );
line(jx,jy,ix,iy); (*
putpixel(ix,iy,c); *)
jx:=ix; jy:=iy; (* garde les
dernieres coordonnees *)
end;
END;
PROCEDURE
TRACER_TOUT(isign:integer);
(* but: tracer le tableau VEC[i],
ensuite faire la TF, puis tracer *)
(* parametres d'entree: isign=1 si FFT
directe, isign=-1 si FFT inverse *)
(* parametres de sortie: isign=1 si FFT
directe, isign=-1 si FFT inverse *)
VAR gd, gm :integer;
BEGIN (* gm:= 0; *)
gd:= Detect; (* particulier au
Turbo-Pascal *)
initGraph(gd,gm,''); (* particulier
au Turbo: initialise le mode graphique *)
setcolor(2); (* particulier au
Turbo: taper Ctrl+F1 pour l'aide *)
TRACER1(0, kpixelYmax div 2); (*
tracer dans la moitie du haut *)
FOURIER_RAPIDE(VEC, kDIM div 2,
isign); (* TF directe ou inverse *)
setcolor(3); (* particulier au
Turbo-Pascal *)
TRACER1(kpixelYmax div 2,
kpixelYmax); (* tracer dans la moitie du bas *)
readln; closegraph; (* particulier
au Turbo *)
END;
PROCEDURE SORTIE;
BEGIN writeln; writlen('Fin du programme SPECTRE venant
du site
"http://www.geocities.com/facultedegenie" page
TOP10'); writeln;
END;
BEGIN (* Programme
Principal *)
REPEAT
ENTREE;
if Choix>0 then
TRACER_TOUT(1); (* F.F.T. directe *)
if Choix>0 then FILTRER; (*
eventuellement *)
if Choix>0 then
TRACER_TOUT(-1); (* F.F.T. inverse *)
UNTIL Choix=0;
SORTIE;
END.
|
(i)
Faculté de Génie, Liban
