FFT2D, IMAGE, SPECTRE, FILTRE PASSE-BAS PASSE-HAUT

//-------------------------------------------------
// FFT 1D - 2D
//-------------------------------------------------


//-------------------------------------------------
// calculs d'une tranformee de Fourier rapide
// non-normalisee sur l'ensemble des nombres
// dans le tableau <tab> indexes par <tab[0]>
// puis <tab[step]> puis <tab[2*step]> ...
// le nombre d'elements doit etre une puissance de deux
static void computing_FFT1D(P_COMPLEX res,P_COMPLEX tab,int step,int n,P_COMPLEX w)
{
if(1 == n)
  {
  // cas trivial
  res->re = tab->re;
  res->im = tab->im;
  }
else
  {
  int     k;
  COMPLEX w2,wPowk;
  MulComplex(&w2,w,w);
  // diviser pour regner
  computing_FFT1D(res,tab,2*step,n/2,&w2); // coefficients pairs
  computing_FFT1D(res+n/2,tab+step,2*step,n/2,&w2); // coefficients impairs
  wPowk.re = 1.;
  wPowk.im = 0.;
  for(k=0;k<n/2;k++)
    {
    // P(X) = A(X²) + X.B(X²)
    //  avec A(pair) et B(impair)
    COMPLEX A,B;
    A = res[k];
    B = res[k + n/2];
    MulComplex(&B,&B,&wPowk);
    AddComplex(res+k,&A,&B);
    SubComplex(res+k+n/2,&A,&B);
    MulComplex(&wPowk,&wPowk,w);
    }
  }
} // computing_FFT1D()


//-------------------------------------------------
// effectue une transformation de Fourier rapide
// normalisee sur le tableau <tab> de taille
// <size>, une puissance de deux
// le resultat sera dans <res>
static void FFT1D(P_COMPLEX res,P_COMPLEX tab,int size,int sgn)
{
COMPLEX   W;
int       k;
double    sqrt_size;

// obtention de la racine <size>-ieme de l'unite = exp(2.i.PI/size)
GetNthRootUnitComplex(&W,sgn,size);
// calcul de la FFT non-normalisee
computing_FFT1D(res,tab,1,size,&W);
// normalisation
sqrt_size = sqrt((double)size);
for(k=0;k<size;k++)
  {
  res[k].re /= sqrt_size;
  res[k].im /= sqrt_size;
  }
} // FFT1D()


//-------------------------------------------------
// effectue une transformation de Fourier en
// deux dimensions sur la matrice <mat>
// qui doit etre de taille n.n avec n une
// puissance de deux
// le resultat sera dans <res>
// on a besoin pour cette fonction d'une matrice
// temporaire et d'un tableau
// cette fonction retourne le module maximum du spectre
static double FFT2D(PP_COMPLEX res,PP_COMPLEX mat,int size,int sgn,PP_COMPLEX tmp,P_COMPLEX buffer)
{
double  A_max2;
int     x,y;

// FFT a une dimension selon les Y
for(x=0;x<size;x++)
  {
  FFT1D(buffer,mat[x],size,sgn);
  for(y=0;y<size;y++)
    {
    tmp[y][x] = buffer[y];
    }
  }
// FFT a une dimension selon les X
A_max2 = 0.;
for(x=0;x<size;x++)
  {
  P_COMPLEX Z;

  FFT1D(res[x],tmp[x],size,sgn);

  Z = res[x];
  for(y=0;y<size;y++)
    {
    double r2;
    r2 = Z->re*Z->re + Z->im*Z->im;
    if(r2 > A_max2)
      {
      A_max2 = r2;
      }
    Z ++;
    }
  }
return sqrt(A_max2);
} // FFT2D()