# Fourier transform for irregular sampling following
# T.J. Deeming
#
# Example: deeming Time Imag 0 5 0.0001 
# Geneva Observatory, March 2 2002, Paul Bartholdi
# Adapted by Laurent Eyer March 1, 2004

deeming 5
  #  
  # $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
   ylabel DFT
   window 1 -3 1 1
   limits f fg box connect f fg
   ylabel Spectral window
   max ff indff maxff 
   xlabel Frequency (max at $(f[$indff]) (= 1/$(1/f[$indff])))
   window 1 1 1 1


deem 3
  # Calculation of the normalized modulus $3f and $3g  of the discrete Fourier transform
  # discrete Fourier Transform (dtf) for the frequencies $3, from Deeming's algorithm.
  # $1:  time
  # $2:  signal
  # $3:  frequencies
  # $3f: Normalized modulus ''\
  #                            > resultats
  # $3g: ........... alias   ./            $3ge same extended to 2*fmax
  # If the variable POWER is defined, the power (for ff and fg) is computed
  # The suffixes _r and _i give the real and imaginary parts
  # of the imaginary vectors

  define  _dtemps (dimen($1))
  define  _dfreq  (dimen($3))
  define  _dt2    ($_dtemps * $_dtemps)
  #
  set $3e  = $3 concat ( $3 + $3[1] + $3[$_dfreq-1] )  # awaited frequencies
  set $3f_r = 0 * $3            # real
  set $3f_i = 0 * $3            # imaginary
  set $3g_r = 0 * $3            # aliasing
  set $3g_i = 0 * $3
  set $3ge_r= 0 * $3e
  set $3ge_i= 0 * $3e
  set $3e2pi = 2 * pi * $3e
  set p = 1 / ( f>0 ? f : f[1] )
  set _n_ = 0,$_dfreq-1
  do it = 0, $_dtemps - 1 {
    set  sit = $2[$it]
    set  _Z_ = $3e2pi * $1[$it]
    set  _ce = cos( _Z_ )
    set  _se = sin( _Z_ )
    set  _c  = _ce[_n_]
    set  _s  = _se[_n_]
    set  $3f_r = $3f_r + _c * sit
    set  $3f_i = $3f_i + _s * sit
    set  $3ge_r = $3ge_r + _ce
    set  $3ge_i = $3ge_i + _se
  }
  set $3f_r  = 2 * $3f_r  / $_dtemps
  set $3f_i  = 2 * $3f_i  / $_dtemps
  set $3ge_r = $3ge_r / $_dtemps
  set $3ge_i = $3ge_i / $_dtemps
  set $3f    = ( $3f_r * $3f_r + $3f_i * $3f_i )
  set $3ge   = ( $3ge_r * $3ge_r + $3ge_i * $3ge_i )
  set $3g    = $3ge[_n_]
  set $3g_r  = $3ge_r[_n_]
  set $3g_i  = $3ge_i[_n_]
  if( $?POWER == 0 ) {
    set $3f  = sqrt($3f)
    set $3g  = sqrt($3g)
    set $3ge = sqrt($3ge)
  }

set_power    #  compute the power in the vector ff  ( and  fg  )
  define POWER 1

clear_power  #  compute the amplitude in the vector  ff  ( and  fg )
  define POWER delete

centre 2     # subtract the mean of the two parameters
  set $1 = $1 - sum($1) / dimen($1)
  set $2 = $2 - sum($2) / dimen($2)

max 3        #  search for the maximum of a vector :  $3 = $1[$2] = max{ $1 }
  set _im LOCAL
  set __s LOCAL
  set _im=0,dimen($1)-1
  set __s=$1
  sort< __s _im >
  define $2 (_im[dimen($1)-1])
  define $3 ($1[$$2])