source: Sophya/trunk/SophyaPI/DemoPIApp/fft2d.cc@ 2600

int nlig = 128; int ncol = 128;      
TMatrix<r_8> in(nlig,ncol);          
int l,c;                             
for(l=40; l<80; l++)                 
  for(c=40; c<80; c++) {             
    double r2 = (l-60.5)*(l-60.5)+(c-60.5)*(c-60.5);   
    in(l,c) = 30.*exp(-r2/9.);        
}                                
//  Le dernier argument de RandomSequence est le sigma du bruit   
//  Si on veut rajouter du bruit, decommenter les 3 lignes suivantes  
  TMatrix<r_8> noise(nlig,ncol);                                
  noise = RandomSequence(RandomSequence::Gaussian, 0., 0.5);     
  in += noise;                                                  
TMatrix< complex<r_8> > fftin;       
FFTWServer ffts;                     
ffts.setNormalize(true);             
cout << " FFTServer info string= " << ffts.getInfo() << endl;
cout << " Calling ffts.FFTForward(in, fftin) : " << endl;    
ffts.FFTForward(in, fftin);                                  
//  On garde les tableau in et fftin avec KeepObj()             
KeepObj(in);                                                 
KeepObj(fftin);                                              
//  On fait un filtre passe-haut     
//  Il faut savoir que les frequences negatives sont rangees    
//  dans la seconde moitie du tableau en numero de ligne        
int kx, ky;                                     
TMatrix< complex<r_8> > fftfb(fftin.NRows(), fftin.NCols());           
TMatrix< complex<r_8> > fftfd(fftin.NRows(), fftin.NCols());           
// fb : Filtrage brutal , fd :  Filtrage doux 
for(ky=0; ky<fftin.NRows(); ky++)             
  for(kx=0; kx<fftin.NCols(); kx++) {
    int kyy = ky;
    // deuxieme  partie du tableau - Frequences negatives en y 
    if (ky > fftin.NRows()/2) kyy = fftin.NRows()-ky;
    double kr = sqrt((double)kx*kx+kyy*kyy);              
    complex<double> att;
    // Filtrage brutal  fonction en escalier 
    if (kr < 10)  att = complex<double>(0.5, 0.);
    else att = complex<double>(1., 0.);
    fftfb(ky, kx) = fftin(ky, kx) * att ;   
    // Filtrage doux sous forme de 1-exp(-alpha k_r)
    att = complex<double>(1.-exp(-kr*0.1), 0.);   
    fftfd(ky, kx) = fftin(ky, kx) * att ;   
  }
TMatrix<r_8> infb(nlig,ncol);          
TMatrix<r_8> infd(nlig,ncol);          
// FFTBackward -- Attention, le tableau des coefficient de Fourier 
// est modifie - On en fait une copie 
TMatrix< complex<r_8> > fftf(fftin.NRows(), fftin.NCols());   
// Pour le filtre brutal
fftf = fftfb;
ffts.FFTBackward(fftf, infb);       
// Pour le filtre doux
fftf = fftfd;
ffts.FFTBackward(fftf, infd);       
//  On garde les tableaux fft filtre et signal filtre
KeepObj(infb);                                             
KeepObj(infd);                                             
KeepObj(fftfb);                                                  
KeepObj(fftfd);