var
  f: array[0..NMax-1] of real; //valors de la funció

Procedure FFT;
Const
  NMax = 8; //Nombre màxim de dades
  MaxStages = 4; //Això ha de ser més gran que el lg2(NMax);

var
  N:integer; //Nombre de dades, nombre dels coeficients de Fourier N més petit que NMAX i cal que sigui una potència de 2, 2, 4, 8, 16, 32, 64, 128, ...
  nStages:integer;
  cFr,cFi: array[0..NMAx-1,0..MaxStages-1] of real; //Coeficients Fourier, part real i part imaginària

  //això serveix per calcular l'ordre en què s'ha de posar les dades de la funció
  ncoe:integer;
  coe,coev: array[1..NMax] of integer; //coeficients dels valors de la x on cal buscar l'f(x)

  s,s1,s2:string; //Per mostrar informació debug


  x,y:integer; //posició dins la matriu de coeficients
  i,l,lp:integer; //Comptadors
  a:real; //això ho fem servir per calcular amb menys decimals, caldria treure-ho. Es per mostrar valors de debug mentre desenvolupem que no sigui massa llargs
  b:integer;
  Wr,Wi: array [0..NMax-1,1..MaxStages-1] of real;
  un: array [0..NMax-1,1..MaxStages-1] of integer;

begin
  n:=NMax;
  //Primera col·lecció de dades per comprovar resultats f(0)=0, f(1)=1, ... https://youtu.be/FaWSGmkboOs
  for y:=0 to N-1 do f[y]:=y;

  //Segona col·lecció de dades per comprovar resultats https://youtu.be/RioJKiSBlyg
  {
  f[0]:=1;
  f[1]:=2;
  f[2]:=3;
  f[3]:=4;
  f[4]:=4;
  f[5]:=3;
  f[6]:=2;
  f[7]:=1;
  }

  //Posem  tots els valors 0
  for y:=0 to N-1 do for x:=0 to MaxStages-1 do cFr[y,x]:=0;
  for y:=0 to N-1 do for x:=0 to MaxStages-1 do cFi[y,x]:=0;

  //log2(dades) Això es pot fer d'una altra manera, però ens assegurem que el nombre de dades és una potència de 2. Màxim 4096 dades. Això es pot canviar per si cal que la finestra de dades, nombre de dades sigui més gran.
  if N=1 then begin
    nStages:=1;
    ncoe:=1;
  end else if N=2 then begin
    nStages:=2;
    ncoe:=2;
  end else if N=4 then begin
    nStages:=3;
  end else if N=8 then begin
    nStages:=4;
  end else if N=16 then begin
    nStages:=5;
  end else if N=32 then begin
    nStages:=6;
  end else if N=64 then begin
    nStages:=7;
  end else if N=128 then begin
    nStages:=8;
  end else if N=256 then begin
    nStages:=9;
  end else if N=512 then begin
    nStages:=10;
  end else if N=1024 then begin
    nStages:=11;
  end else if N=2048 then begin
    nStages:=12;
  end else if N=4096 then begin
    nStages:=13;
  end else begin
    ShowMessage('El nombre de dades no és una potència de 2 <= que 4096');
    exit;
  end;

  //Càlcul de l'ordre en que s'han d'avaluar les f. Això es fa posant el nombres en binari, inveertint i reordenant, però ho faig així perquè és entendor
  coe[1]:=0;
  coe[2]:=1;
  coev[1]:=0;
  coev[2]:=1;
  if nStages>2 then begin
    ncoe:=2;
    for i:=3 to nStages do begin
      for l:=1 to ncoe do begin
        coe[l]:=2*coev[l];
        coe[ncoe+l]:=2*coev[l]+1;
      end;
      for l:=1 to ncoe do begin
        coev[l]:=coe[l];
        coev[ncoe+l]:=coe[ncoe+l];
      end;
      ncoe:=2*ncoe;
    end;
  end;
  if N<>ncoe then begin
    ShowMessage('No coincideix el nombre de coeficients amb el nombre de dades, error no cotrolat');
    exit;
  end;

  //Això és informació per debug
  s:=IntToStr(N)+' - '+IntToStr(ncoe)+' - ';
  for i:=1 to N do s:=s+IntToStr(coe[i])+', ';
  ShowMessage(s); //ordre dels valors de les n que es posen a la fórmula. Per exemple si n=8, 0, 4, 2, 6, 1, 5, 3, 7, és l'ordre que cal posar a la fórmula

  //Posem els valors del coeficients del primer stage amb els valors f ordenats
  for i:=1 to N do cFr[i-1,0]:=f[coe[i]];
  //Això és informació per debug
  s:='';
  for i:=0 to N-1 do s:=s+RealToString(cFr[i,0])+', ';
  ShowMessage(s); //valors de les f. Per exemple si n=8, f(0), f(4), f(2), f(6), f(1), f(5), f(3), f(7).

  //calcula coeficients que multipliquen
  l:=2; //Stage
  lp:=1; //Stage petit, meitat de l'anterior, pas
  for x:=1 to nStages-1 do begin
    for y:=0 to N-1 do begin
      if y mod l >= lp then b:=-1 else b:=1;
      un[y,x]:=b;
      if b=1 then begin
        Wr[y,x]:=1;
        Wi[y,x]:=0;
      end else begin
        a:=cos(2*pi*(y mod lp)/l);
        if abs(a)<0.00001 then a:=0;
        Wr[y,x]:=a;
        a:=sin(2*pi*(y mod lp)/l);
        if abs(a)<0.00001 then a:=0;
        Wi[y,x]:=a;
      end;
    end;
    lp:=l;
    l:=l*2;
  end;

  //Per mostrar resultats.
  FormPrincipal.MemoResultats.Text:='';
  //FFT
  l:=2; //Stage
  lp:=1; //Stage petit, meitat de l'anterior, pas
  x:=0;
  while x0 then s2:='+'+s2+'i' else s2:=s2+'i';
    s:=s+s1+' '+s2+#9;
    FormPrincipal.MemoResultats.Lines.Append(s);
  end;
end;