clear;

////////////////////////////
// Partie 1: TFD rapide  ///
////////////////////////////


//Fonction permettant de realiser le decoupage
// d'un vecteur X de taille pair en deux vecteurs
// Le premier aura les valeurs des indices impair
// Le second aura les valeurs des indices pair
function [X0, X1] = decouper(X);
 dimensionX = length(X);
 for k=0:dimensionX/2-1;
   X0(k+1) = X(2*k+1);
   X1(k+1) = X(2*k+2);
 end;
endfunction;

//Tests de la fonction
[dec1, dec2] = decouper([1,2,4,5]');


//Transformee Fourier Discrete Rapide
// Precondition: X de dimension 2^q
function T= TFDrapide(X)
 dimensionX = length(X);
 if (dimensionX==1) then
    //cas d'arret
    T(1) = X(1);
  else
    //cas general
    [X0, X1] = decouper(X);
    TFDX0 = TFDrapide(X0);
    TFDX1 = TFDrapide(X1);
    for k=0:dimensionX/2-1;
      T(k+1) = TFDX0(k+1) + exp(-2* %i * %pi * k/ dimensionX) * TFDX1(k+1);
      T(k+1+dimensionX/2) = TFDX0(k+1) - exp(-2* %i * %pi * k/ dimensionX) * TFDX1(k+1);
    end;
  end;
endfunction;



///////////////////////////////////
// Partie 2: TFD rapide inverse  //
///////////////////////////////////


//Transformee Fourier Discrete Inverse Rapide
function T = TFDrapideInverse(X)
 dimensionX = length(X);
 T = TFDrapide(X);
 for k=0:dimensionX-1;
   T(k+1) = conj(T(k+1)) / dimensionX;
 end;
endfunction;



//////////////////////////////////
// Partie 3: Module de spectre  //
//////////////////////////////////


//Calcul le module du spectre d'un vecteur X
//  qui sera allonge a une taille 256
//Precondition: taille de X inferieur a 256
function S = modspectre(X)
  dimensionX = length(X);
  for k=1:256;
    S(k) = X(max(int((k /256) * dimensionX),1));
  end;
  S = TFD(S);
  for k=1:256;
    S(k) = abs(S(k));
  end;
endfunction;


/////////////////////////////
// FONCTION POUR LES TESTS //
/////////////////////////////

//version moins optimise du calcul de la TFD
// par simple application de la definition de la TFD
function T = TFD(X);
  dimensionX = length(X);
 for k=0:dimensionX-1;
   somme=0;
   for n=0:dimensionX-1;
     somme = somme + X(n+1)* exp(-2*%i*%pi*k*n/dimensionX);
   end;
   T(k+1) = somme;
 end;
endfunction;


//version moins optimise du calcul de la TFD inverse
// par simple application de la definition de la TFD inverse
function T = TFDinverse(X);
 dimensionX = length(X);
 for k=0:dimensionX-1;
   somme=0;
   for n=0:dimensionX-1;
     somme = somme + X(n+1)* exp(2*%i*%pi*k*n/dimensionX);
   end;
   T(k+1) = somme / dimensionX;
 end;
endfunction;


//fonction utilise pour tester la fonction modspectre
function X = creerSin()

 for k=1:2048;
   X(k) = sin(k/5);
 end;

endfunction;



////////////
// TESTS ///
////////////

//version rapide
timer();
 s3 =  TFDrapide(creerSin());
//version peu optimise
disp(timer());

timer();
S3 = TFD(creerSin());
disp(timer());

//version de scilab

 sc = fft( [45,14,14,15,14,15,16,17]',-1);

//version rapide
 //s4 = TFDrapideInverse([45,14,14,15,14,15,16,17]');
//version moins optimise
 //s4 = TFDinverse([45,14,14,15,14,15,16,17]');

//test de modspectre
 spectre = modspectre(creerSin());