{
Philippe Langevin, Octobre 1996.

Ce programme est une implantation en Turbo-Pascal
sous MSDOS de l'algorithme de transformée de Fourier
rapide, ce qui explique la gestion anormale de la mémoire.
}

const maxlog=9;
      n= 1 shl maxlog;
      d= n shr 1;

{La procédure Fourier transforme des séquences
de longueur n qui ÂBSOLUMENT etre une puissance de 2.
(64,128,256,512...)}

type Complexe=record
              re,im:double;
	      end;
    

     Sequence=array[0..n-1] of complexe;
      
    
var        
	Image:array[0..n-1] of ^sequence;
         
{l'image est vue comme un tableau de ligne}
         
           inv:array[0..n-1] of word;

{Ce tableau est utilisé par l'algorithme de TFR,
il doit être initialisé de sorte que inv[k] soit
le miroir de k, c'est ce que fait la procédure qui suit}

        
Procedure LesInverses(a,x,y,b:word);
begin
if a>d   then inv[x]:=y                
	 else begin
              LesInverses(a shl 1,x or a, y or b, b shr 1);
              LesInverses(a shl 1, x, y, b shr 1);
	      end;
end;

Procedure Produit(var r:complexe;a,b:complexe);
{Produit de deux nombres complexes}
begin
r.re:=a.re*b.re-a.im*b.im;
r.im:=a.re*b.im+a.im*b.re;
end;
 
Procedure Fourier(var a:Sequence);
{Calcule la TFR de a sur place, la transformée inverse
s'obient en changeant -2 en +2}
var s,i,m:word;
        k:word;
      g,z:Complexe;
    t1,t2:Complexe;
        r:double;
begin
for s:=0 to n-1 do
   if s<=inv[s] then
		begin
		z:=a[s];
                a[s]:=a[inv[s]];
                a[inv[s]]:=z;
                end;

for s:=1 to maxlog do
    begin
    m:=1 shl s;
    z.re:=cos(-2*pi/m);
    z.im:=sin(-2*pi/m);
    g.re:=1;
    g.im:=0;
    for i:=0 to (m shr 1)-1 do
    	begin
        k:=i;
        while (k<n) do
        	begin
                (* Papillon *)
                produit(t1,a[k+(m shr 1)],g);
                t2:=a[k];
                a[k].re:=t2.re+t1.re;
                a[k].im:=t2.im+t1.im;

                a[k+(m shr 1)].re:=t2.re-t1.re;
                a[k+(m shr 1)].im:=t2.im-t1.im;
                k:=k+m;
		end;
        Produit(g,g,z);
        end;
    end;

r:=sqrt(n);
for i:=0 to n-1 do
	begin
	a[i].re:=a[i].re/r;
        a[i].im:=a[i].im/r;
        end;
end;



Procedure Transposition;
{transposition du tableau image}
var i,j:word;
      z:Complexe;
begin
for i:=0 to n-1 do
	for j:=i+1 to n-1 do
        	begin
                z:=Image[i]^[j];
                Image[i]^[j]:=Image[j]^[i];
		Image[j]^[i]:=z;
                end;

end;

 

Procedure FourierLigne;
var i:word;
begin
writeln('Transformée de Fourier...');
for i:=0 to n-1 do Fourier(Image[i]^);
        
end;

 
 

begin

{ chargement de l'image,
n'oubliez pas d'allouer vos pointeurs :-)
}

LesInverses(1,0,0,n shr 1);

{la transformée de l'image nxn s'obtient en quatre étape }

FourierLigne;Transposition;FourierLigne;Transposition;

{traitement de l'image,
n'oubliez pas de liberer vos pointeurs!}
end.