* fichier :  tfr.dgibi
************************************************************************
* Section : Mathematiques Traitement du signal
* Section : Dynamique
************************************************************************
*
*   Transformee de Fourier rapide (FFT via operateur TFR)
*   et transformee inverse  (FFT-1 via operateur TFRI)
*   
*   BP, 2018-10-04
*
************************************************************************

************************************************************************
* Test de plusieurs signaux
************************************************************************

* nombre de points consideres (TFR demande une puissance de 2)
*   p = 10;
  p = 18;
  n = 2**p;

* listes des temps (de dime n)
  t = prog 0. pas (4./n) npas (n-1);
* quelques listes utiles 
  c180 = cos (180*t);
  s180 = sin (180*t);
* creneau
  cre_c = masq c180 'EGSUPE' 0.;
  cre_s = masq s180 'EGSUPE' 0.;
* bruite
  xh = (0.7*(sin (100*t))) + (sin (360*t));
  xb = brui 'BLAN' 'UNIF' 0. 0.5 (n );
  bru1 = xh + xb;
* rampe 
  ram1  = 1.*t;
* pour le trace  
  tt = tabl; tt . 2 = mot 'TIRR';
  tt . 3 = mot 'POIN ';
  tt . 4 = mot 'MARQ SS ROND';
* erreur
  err_p = prog ;

* boucle sur les signaux a traiter  
REPE BB 6;

* quel cas traite t'on ?
  si (&BB ega 1) ; x = c180; finsi;
  si (&BB ega 2) ; x = s180; finsi;
  si (&BB ega 3) ; x = cre_c; finsi;
  si (&BB ega 4) ; x = ram1 ; finsi;
  si (&BB ega 5) ; x = bru1 ; finsi;
  si (&BB ega 6) ; x = cre_s; finsi;

* signal de depart  
  ev_x = evol azur manu 't' t 'x' x;

* Transformee de Fourier 
* Tx = TFR ev_x p 'MOPH';
  Tx = TFR ev_x p 'REIM';
*   si (&BB ega 3); list  Tx; sinon; list resum Tx; finsi;

* Transformee de Fourier inverse 
  ev_y = (TFRI Tx) COUL 'BRON';
*   si (&BB ega 3); list  ev_y; sinon; list resum ev_y; finsi;
*   si (&BB ega 3); dess (ev_x et ev_y) tt; finsi;
  
* test de bon fonctionnement  
  t_y = EXTR (EXTR ev_y 'COUR' 1) 'ABSC' ;
  x   = EXTR (EXTR ev_x 'COUR' 1) 'ORDO' ;
  y   = EXTR (EXTR ev_y 'COUR' 1) 'ORDO' ;
  err_x = MAXI 'ABS' (x - y)  ;
  err_p = err_p et err_x;

  MESS '>>> Cas' ' ' &BB ' : erreur =' err_x;
  
FIN BB ;

TEMP IMPR MAXI CPU;
xCPU1 = TEMP 'CPU' ;
mess '>>> temps CPU = ' xCPU1 ' ms ';

************************************************************************
* Test des options
************************************************************************
* Transformee de Fourier Module PHase
  Tx1 = TFR ev_x p 'MOPH' BLEU;
* Transformee de Fourier option FMAX
  Tx2 = TFR ev_x p 'FMAX' 256. VERT;
* Transformee de Fourier option FMIN (pas souvent utilise a priori)
  Tx3 = TFR ev_x p 'FMIN' 10.5 ORAN;
* Transformee de Fourier options  FMIN FMAX simultanee
  Tx4 = TFR ev_x p  'FMIN' 10.5 'FMAX' 99.5 BRUN;

* Transformee de Fourier inverse sur Tx1 ...
  ev_y1 = (TFRI Tx1) COUL 'BLEU';
  y1 = IPOL t ev_y1;
  err_x1 = MAXI 'ABS' (x - y1)  ;
  err_p = err_p et err_x1;
* ... et Tx2 seulement
  ev_y2 = (TFRI Tx2) COUL 'VERT';
  y2 = IPOL t ev_y2;
* pour ce cas, phenomene de Gibbs visible :
* on calcule une erreur moyenne plutot que max et on est + tolerant  
*   err_x2 = MAXI 'ABS' (x - y2)  ;
  err_x2 = (SOMM (ABS (x - y2))) / (dime x);  
*   si (p < 16); dess (ev_x et ev_y et ev_y1 et ev_y2) tt; finsi;
  

************************************************************************
* TEST DE NON REGRESSION
************************************************************************
* pour info : lors de la creation du cas-test,
*               +----------+-------------------------+--------------+
*               |       linux 64                     |     AIX      |
*               +----------+-------------------------+-------+------+
*               | castem18 | castem apres EVOL #9942 | avant | apres|
* +-------------+----------+-------------------------+-------+------+
* | precision   | ~1.E-10  | ~1.E-15                 |~1.E-10|~1E-15|
* | performance |   540ms  |  190ms                  | 1090ms| 930ms|
* +-------------+----------+-------------------------+-------+------+
* perf linux 32 apres EVOL : 306ms 

* PRECISION
XPREC = VALE 'PREC';  mess 'precision=' XPREC;
XTOL  = XPREC * 100;
ZPREC1 = (maxi 'ABS' err_p) > XTOL;
ZPREC2 = err_x2 > 1.E-2;
SI ZPREC1; MESS 'PB DE PRECISION !'; FINSI;
SI ZPREC2; MESS 'PB DE PRECISION avec option FMAX !'; FINSI;

* * PERFORMANCE
* ZPERF1 = xCPU1 > ??? difficile car tres dependant-machine ;
* SI ZPERF1; MESS 'PB DE PEROFRMANCE !'; FINSI;
* SI (ZPREC1 OU ZPERF1 ou ZPREC2);
SI (ZPREC1 OU ZPREC2);
  ERRE 5;
FINSI ;


FIN ;