deeming 5
  # Written by Paul Bartholdi, adapted by Laurent Eyer March 1, 2004
  # $1 : Time
  # $2 : Signal
  # $3 : minimum searched frequency
  # $4 : maximum searched frequency 
  # $5 : frequency step

  centre $1 $2 
  set f=$3,$4,$5
  deem $1 $2 f
  expand 1.001
  erase
   window 1 3 1 3 
   limits $1 $2 
   limits $1 $fy2 $fy1
   box points $1 $2
   xlabel $1
   ylabel $2
     
   window 1 -3 1 2 
   limits f ff box 0 2 connect f ff
   window 1 -3 1 1
   limits f fg box connect f fg
   max ff indff maxff 
   xlabel Frequency (max at $(f[$indff]) (= 1/$(1/f[$indff])))
   window 1 1 1 1


centre 2   # Subtraction of the mean for the two parameters
  stats $1 $1_mean $1_sig $1_kur
  stats $2 $2_mean $2_sig $2_kur
  set $1 = $1 - $$1_mean
  set $2 = $2 - $$2_mean

deem 3
  # calcul du module normalise ff et fg  de la transformee de fourier
  # discrete (dft) pour les frequences $3, selon l'agorithme de Deeming.
  # $1:  time
  # $2:  signal
  # $3:  frequencies
  # $3f: normalised modulus ''\
  #                            > results
  # $3g: ........... alias   ./              $3ge same extended to 2*fmax
  #
  define  pi2 (2*pi)
  define  _dtemps (dimen($1))
  define  _dfreq  (dimen($3))
  define  _dt2    ($_dtemps * $_dtemps)
  #
  set fze = $3 + $3[$_dfreq-1]  IF( $3>0 )
  set fe  = $3 concat fze     # extended frequencies
  set ffr = 0 * $3            # real
  set ffi = 0 * $3            # imaginary
  set fgr = 0 * $3            # aliasing
  set fgi = 0 * $3
  set fger= 0 * fze
  set fgei= 0 * fze
  set $1p  = $1 * $pi2         # reutilise ci-dessous
  do it = 0, $_dtemps - 1 {
    define tit ($1p[$it])
    define sit ($2[$it])
    set  _c = cos( $3 * $tit )
    set  _s = sin( $3 * $tit )
    set  _ce = cos( fze * $tit )
    set  _se = sin( fze * $tit )
    set  ffr = ffr + _c * $sit
    set  ffi = ffi + _s * $sit
    set  fgr = fgr + _c
    set  fgi = fgi + _s
    set  fger = fger + _ce
    set  fgei = fgei + _se
  }
  set ff = 4 * ( ffr * ffr + ffi * ffi ) / $_dt2
  set fg = ( fgr * fgr + fgi * fgi ) / $_dt2
  set fge = fg concat ( ( fger * fger + fgei * fgei ) / $_dt2 )
  #  delete _c   delete _s   delete fz   delete $1p
  #  delete _ce  delete _se  delete fze
  if( $?POWER == 0 ) {  set ff  = sqrt(ff)  set fg  = sqrt(fg)  set fge = sqrt(fge) }

set_power    #  calcule la puissance du signal dans  ff  ( et  fg  )
  define POWER 1

clear_power  #  calcule l'amplitude du signal dans  ff  ( et  fg )
  define POWER delete

max 3     #  cherche le max d'un vecteur :  $3 = $1[$2] = max{ $1 }
  define __dmax (dimen($1)-1)
  set _im=0,$__dmax
  set __s=$1
  sort { __s _im }
  define $2 (_im[$__dmax])   define $3 (__s[$__dmax])
  delete __s  delete _im