//genesis - spectra_funcs2.g

/* Functions for calculating/displaying spectra

   Assumes global definitions for:
   float tmax   // simulation time in sec
   float out_dt // time step for output from the data source
   int graphics, debug //flags for using graphics or printing extra info

   A typical usage is:

   str timedata ="/time_data" // name of table object for data to analyze
   str data_source = "/EPSCsummer"  // name of object that generates time data
   include spectra_funcs2.g
   setclock  1  {out_dt}
   make_EPSCsummer
   // assume that make_Vmgraph created a graph to plot sum of EPSCs
   addmsg /EPSCsummer /data/Iksum PLOT output *EPSC_sum *red
   make_time_data_table
   addmsg {data_source} {timedata} INPUT output
   make_FTgraph {x} {y}
*/

/* If use_fft_processor =  0, the Fourier Transform will be performed with a
   GENESIS script function to perform the slow Discrete Fourier Transform.

   Otherwise, the function plot_FT invokes an external compiled program
   fft_processor (http://www.arachnoid.com/signal_processing/fft.html).
   It communicates with the program via files {fft_infile} and {fft_outfile}.
*/

int use_fft_processor = 0

int debug = 1

setclock  1  {out_dt}

function make_FTgraph(xpos, ypos)
    str form = "/spectra_form"
    float xpos, ypos
    float xmin = 0; float xmax = 200;
    float ymin = 0; float ymax = 1;
    create xform {form} [{xpos},{ypos},400,380]
    pushe {form}
    create xgraph power_spectrum  -hgeom 70% -bg white
    setfield ^ xmin {xmin} xmax {xmax} ymin 0 ymax {ymax}
    create xbutton do_FT -script calc_spectra \
	-label "Calculate frequency spectra of network activity"
    create xdialog min_time -wgeom 50% -label "Minimum time (sec)" \
	-value 0.0
    create xdialog max_time -wgeom 50% -ygeom 0:do_FT -xgeom 0:min_time \
        -label "Maximum time (sec)" -value {tmax}
    create xdialog max_freq -label "Maximum frequency (Hz)" \
	-value 200
    create xbutton DISMISS -script "xhide "{form}
    pope
    makegraphscale  {form}/power_spectrum
    xshow /spectra_form
end

// make table to hold time series data to Fourier analyse
function make_time_data_table
    int xdivs = {round {2*tmax/{getclock 1}}} + 2 // twice as large as needed
    float xmin = 0; float xmax = {tmax} // this is arbitrary
    if ({exists {timedata}}) // Get rid of any existing one with old xdivs
        delete {timedata}    // table doesn't seem to have a TABDELETE action
    end
    create table {timedata}
    call {timedata} TABCREATE {xdivs} {xmin} {xmax}
    setfield {timedata} step_mode 3
    useclock {timedata}  1
    // One could add a RC element here between {data_source} and {timedata}
    // as low-pass filter, but aliasing of high frequency spectral components
    // doesn't seem to be problem
end

/* Calculate and plot the spectra using a Fast Fourier Transform program
   fft_processor (http://www.arachnoid.com/signal_processing/fft.html).
   If this isn't available, set use_fft_processor =  0 to use a GENESIS
   script function to perform the slow Discrete Fourier Transform.
*/
function plot_FT(timedata, tmin, tmax, fmax)
    str timedata
    float tmin, tmax, fmax, delta_t, delta_f
    float fmin = 0.0
    float h_k, a_n, b_n, p_n, twoPIbyN, scalefactor, pmax
    int i, n, k, kmin, kmax, Npts, Nfreqs
    int maxindex = {getfield {timedata} output}  // index of last entry
    // Use the largest even number for Npts, in order to calc Npts/2 freqs
    int maxpts = 2*{trunc {(maxindex + 1)/2.0}}
    delta_t = {getclock 1}
    kmin = {trunc {tmin/delta_t}}
    kmax = {round {tmax/delta_t}}
    if ({kmax} > {maxpts})
            kmax = maxpts
    end
    delta_f = 1.0/((kmax - kmin)*delta_t)
    Npts = kmax - kmin + 1
    // Limit the range of frequencies calculated
    // fmax must be <= maxpts*delta_f -- this is not checked for!
    Nfreqs = {round {fmax/delta_f}} + 1
    if (debug)
      echo "Index of last time step: " {kmax}
      echo {Npts} " time points from " {tmin} " to " {tmax}
      echo {Nfreqs} " frequency points from 0 to " {(Nfreqs - 1)*delta_f}
      echo "time interval: " {delta_t} " sec; freq interval: " {delta_f} " Hz."
     end

    /* If graphics are used, set up the plots */
    if (graphics)
        str spectraplot = "/spectra_form/power_spectrum/FTpower"
        float ymin = 0; float ymax = 1.0;
        setfield /spectra_form/power_spectrum xmin {fmin} xmax {fmax} \
		ymin 0 ymax {ymax}
        if ({exists {spectraplot}})
	    delete {spectraplot}
        end
        create xplot {spectraplot}
	setfield {spectraplot} fg blue ysquish 0
        call /spectra_form/power_spectrum RESET
    end  // if (graphics)

    if (use_fft_processor) // Use the external program "fft_processor"
      str fft_infile = "fft_infile.txt"
      str fft_outfile = "fft_outfile.txt"
      int array_size // size of array given to fft_processor
      int sample_rate // Sample rate in Hz (integer for fft_processor)
      float new_delta_f // Frequency interval in {fft_outfile}
      int new_Nfreqs // Number of frequencies to plot from {fft_outfile}
      str ablist // Will hold a line read from {fft_outfile}

      // Make a new table of the size to hold the data from tmin thru tmax
      create table /temp_data

      /* Allocate space with xdivs xmin xmax.  Here xmin/xmax are relative to
         the start and end of the selected time range. i.e. xmin = 0,
         xmax = tmax - tmin.
      */
      call /temp_data TABCREATE {Npts - 1} 0 {(Npts - 1)*out_dt}

      i = 0 
      for (k = {kmin}; k <= {kmax}; k = k +1) // Loop over desired time interval
            h_k = {getfield {timedata} table->table[{k}]}
            setfield /temp_data table->table[{i}] {h_k}
            i = i + 1
      end

      // Number of points in the output file have to be a power of 2
      // Horrible curly brackets have to be just right!
      array_size = {pow 2 { {trunc { {log {Npts} }/{log 2} } } + 1 } }
      sample_rate = {round  {(array_size/out_dt)/Npts} }
      if (debug)
        echo "array_size = "{array_size}
        echo "sample_rate = "{sample_rate}
        echo "Before expansion /temp_data last entry " {Npts-1}
        echo {getfield /temp_data table->table[{Npts-1}]}
      end
      // Expand the table to array_size
      call /temp_data TABFILL {array_size - 1} 0

      openfile {fft_infile} w
      writefile  {fft_infile} {array_size}
      writefile  {fft_infile} {sample_rate}
      for (i = 0; i < {array_size}; i = i +1)
          writefile {fft_infile} {getfield /temp_data  table->table[{i}] } 0.0
      end
      closefile {fft_infile}
      if (debug)
        echo "After expansion /temp_data last entry " {array_size - 1}
        echo {getfield /temp_data table->table[{array_size-1}]}
      end
  
      delete /temp_data // delete table, allowing repeated calls to plot_FT

      // Unix shell command using fft_processor.  If this line produces an
      // error, set use_fft_processor = 0, and see the comments earlier
      cat fft_infile.txt | fft_processor > {fft_outfile}

      // now read the resulting fft_outfile
      openfile {fft_outfile} r
      array_size = {readfile  {fft_outfile}}
      sample_rate = {readfile  {fft_outfile}}
      new_delta_f = 1.0*sample_rate/array_size  // cast to float
      new_Nfreqs = {trunc {Nfreqs*delta_f/new_delta_f}}

      if (debug)
        echo "array_size = "{array_size}
        echo "sample_rate = "{sample_rate}
        echo "new_delta_f = "{new_delta_f}
        echo "new_Nfreqs = "{new_Nfreqs}
      end

      // Will store power spectrum in a table
      create table /temp_spectra
      call /temp_spectra TABCREATE {new_Nfreqs - 1} 0 \
          {(new_Nfreqs - 1)*new_delta_f}
      pmax = 0.0
      ablist = {readfile {fft_outfile} -l}  // ignore the zero freq data
      setfield /temp_spectra table->table[0] 0 // replace it with 0

      for (n = 1; n < {new_Nfreqs}; n = n + 1)
        // This is a trick to separately get the two values on one line
        // without doing another read
        ablist = {readfile {fft_outfile} -l}
        a_n = {getarg {arglist {ablist}} -arg 1}
        b_n = {getarg {arglist {ablist}} -arg 2}
        p_n = a_n*a_n + b_n*b_n
        setfield /temp_spectra table->table[{n}] {p_n}
        pmax = {max {pmax} {p_n}} // Find spectrum max for normalization
      end // loop over frequencies

      scalefactor = 1.0/pmax
      if (debug)
           echo "max p_n = " {pmax}
      end
      closefile {fft_outfile}

      if (graphics)
        for (n = 0; n < {new_Nfreqs}; n = n + 1)
          p_n = scalefactor * {getfield /temp_spectra table->table[{n}]}
          call {spectraplot} ADDPTS {n*new_delta_f} {p_n}
        end
      end
      delete /temp_spectra

    else  // Use script function for DFT
      /* Calculate the spectra with a "Slow Fourier Transform".  This is
	 much slower than the Fast Fourier Transform algorithm, which
	 cleverly combines terms in the double sum in the Discrete Fourier
	 Transform.  But it is simple to implemement in GENESIS, and doesn't
         require Npts to be a power of 2.
      */
      // Will store power spectrum in a table
      create table /temp_spectra
      call /temp_spectra TABCREATE {Nfreqs - 1} 0 {(Nfreqs - 1)*delta_f}
      pmax = 0.0
      twoPIbyN = 2*3.14159265/Npts
      if (debug)
          echo {getdate}
      end
      // Now do the slow double summation
      // First, sum over frequencies, skipping the zero frequency component
      setfield /temp_spectra table->table[0] 0 // replace it with 0
      for (n = 1; n < {Nfreqs}; n = n + 1)
        a_n = 0.0
        b_n = 0.0
        for (k = {kmin}; k <= {kmax}; k = k +1) // Sum over time data
            h_k = {getfield {timedata} table->table[{k}]}
            a_n = a_n + h_k*{cos {twoPIbyN*n*k}}
            b_n = b_n + h_k*{sin {twoPIbyN*n*k}}
        end // time loop
        p_n = a_n*a_n + b_n*b_n
        setfield /temp_spectra table->table[{n}] {p_n}
        pmax = {max {pmax} {p_n}} // Find spectrum max for normalization
      end // loop over frequencies
      if (debug)
          echo {getdate}
      end
      scalefactor = 1.0/pmax
      if (graphics)
        for (n = 0; n < {Nfreqs}; n = n + 1)
          p_n = scalefactor * {getfield /temp_spectra table->table[{n}]}
          call {spectraplot} ADDPTS {n*delta_f} {p_n}
        end
      end
    end  // if(use_fft_processor)
end

function calc_spectra
    str dialog = "/spectra_form"
    float tmin = {getfield {dialog}/min_time value}
    float tmax = {getfield {dialog}/max_time value}
    float fmax = {getfield {dialog}/max_freq value}
    plot_FT {timedata} {tmin} {tmax} {fmax}
end