#ifndef NEXT

#include "snd-ufun.h"
#include <stdlib.h>
#include <stdio.h>
#include <unistd.h>
#include <math.h>
#include <string.h>

#define DEBUGGING 0


/* -------------------------------- EXTERNAL FUNCTION EXAMPLES -------------------------------- 
 * these are compiled and loaded into a shared library, accessible in Snd at run time via C-X C-L (dlopen/dlsym) 
 * see Snd's makefile for an example of how to make a shared object file.
 *
 * the examples here are:
 *
 *    rms -- a simple analysis function (returns rms of segment)
 *    fft_filter -- a simple editing function (smooths using fft-based filtering)
 *    autocorrelation -- a simple display function (performs autocorrelation on the visible data)
 * if DEBUGGING
 *    times2 -- grab 1024 points from cursor and multiply by 2
 *    extlist -- grab 128, delete, insert 256, make region of second 128 (demo of list of actions)
 *
 *
 * Each external function is of the form void *func(void *ptr).
 *
 * When called, the user function is passed a pointer to a snd_to_ufun structure,
 *    defined in snd-ufun.h.  The fields are:
 *
 *    snd_next_sample: a function that returns the next sample of the current channel's data
 *    snd_previous_sample: same but the previous sample is returned
 *    snd_eof: returns true when the end of the channel's data is reached.
 *       Each of these functions takes as its only argument the channel-associated pointer
 *       in the array of pointers:
 *    snd_ptrs: snd_data pointers (Snd internal format for its data acess)
 *    chans: the number of channels currently participating in the function call
 *    count: the numeric argument typed before the C-X C-H command or equivalent (i.e. C-U arg)
 *    nargs: the number of arguments typed in the function invocation (i.e. a_fun(4,3.14))
 *    args: the arguments themselves (floats)
 *    wbegs: the window starting sample (also the position of the snd_ptrs readers upon invocation)
 *    wends: window end sample numbers
 *    wcursors: cursor sample numbers
 *
 * The user function then calls snd_next_sample repeatedly to get as much data as it wants,
 * does some operation thereupon, and returns a pointer to a ufun_to_snd structure (defined
 * in snd-ufun.h).  It fields are:
 *    op -- the operation desired from Snd -- see below
 *    chans -- number of channels of returned data (normally the same as the input chans)
 *    lens -- returned data length per channel
 *    begs -- where this data is supposed to begin in the input channels (if relevant)
 *    results -- an array of arrays of floats containing the channel-related results (if not in a file)
 *    filename -- the file (a normal sound file) in which the results can be found (if not in an array)
 *    result_type -- either UFUN_ARRAY, UFUN_FILE, or UFUN_LIST
 *      the file option is intended for cases involving a lot of data (too much for in-core buffers)
 *
 * The 'op' field is currently one of:
 *
 *    SND_NO_OP -- no action taken
 *    SND_REPORT_RESULTS -- results are displayed in the associated minibuffers (only array results)
 *    SND_CHANGE_RESULTS -- results change (replace) current channel data
 *    SND_INSERT_RESULTS -- results are inserted into current channel data
 *    SND_ADD_RESULTS -- results are mixed into current channel data
 *    SND_DELETE -- some portion of current channel data is deleted (results and result_type ignored)
 *    SND_DISPLAY_RESULTS -- results are displayed alongside the time domain and fft windows (only array results)
 *    SND_REGIONIFY_RESULTS -- results are saved as a region
 *    SND_CHANGE_AND_REGIONIFY_RESULTS -- results replace current data and are saved separately as a region
 *    SND_INSERT_AND_REGIONIFY_RESULTS -- results are inserted into current data and saved as a region
 *    SND_ADD_AND_REGIONIFY_RESULTS -- results are mixed into current data and saved as a region
 *
 * The external function can also return a list of these structures, each of
 * which will be processes in its turn by Snd.  The list is of the ufun_list_to_snd
 * type, where:
 *
 *    result_type -- SND_UFUN_LIST
 *    len -- number of ufun_to_snd structs involved
 *    ops -- an array of ufun_to_snd structs, each being one of the things described above.
 *
 * the somewhat peculiar data accessors in the snd_to_ufun structure are made necessary by shared object limitations.
 */

/* -------------------------------- RMS -------------------------------- */
/* rms: an example of an external analysis function (does not edit data or add to displays) */
/* 
 * after begin loaded via C-X C-L library.so rms, can be invoked via 
 *  C-X y or C-X C-Y: autocall every time X axis changes (- as arg = cancel auto call)
 *  C-X h or C-X C-H: call once -- this prompts for the desired external function
 *  C-h or C-H: call again without prompt (sort of like C-S and C-R)
 *  C-X h/H: call on selection
 *  the versions that do not have the control key on the second letter apply to the current selection.
 */

void *rms(void *ptr)
{
  ufun_to_snd *sf;                                        /* results passed back to Snd */
  snd_to_ufun *sp = (snd_to_ufun *)ptr;                   /* void pointer to simplify linkages */
  /* sp->nargs is the number of arguments in-coming, each in the float array sp->args */
  int i,j,len;
  float val,sum;
  isnd_fun *at_eof;
  fsnd_fun *next_sample;
  at_eof = sp->snd_eof;
  next_sample = sp->snd_next_sample;
  sf = (ufun_to_snd *)calloc(1,sizeof(ufun_to_snd));      /* create the structure to hold our results */
  sf->results = (float **)calloc(sp->chans,sizeof(float *));
  sf->lens = (int *)calloc(sp->chans,sizeof(int));
  for (i=0;i<sp->chans;i++)                               /* sp->chans = how many sync'd sounds are participating */
    {
      len=0;                                              /* how many samples are include in the sum */
      j=0;                                                /* if count argument used, j stops us when it is reached */
      sum = 0.0;                                          /* our running sum */
      sf->results[i] = (float *)calloc(1,sizeof(float));  /* will hold rms result for this channel */
      sf->lens[i] = 1;                                    /* length of returned results for this channel */
      while (!((*at_eof)(sp->snd_ptrs[i])))               /* eof at end of file or region etc */
	{
	  val = ((*next_sample)(sp->snd_ptrs[i]));        /* next_sample returns the next (current edit-wise) sample */
    	                                                  /* previous_sample returns the previous -- see snd-edits.c for others */
	  sum += (val*val);                               /* sum of the squares as per rms definition */
	  len++;
	  j++;
	  if ((sp->count > 1) && (j >= sp->count)) break; /* numeric arg (if any) = segment length */
	}
      if (len == 0) len = 1;                              /* perhaps we're already at EOF */
      sf->results[i][0] = sqrt(sum/(float)len);           /* rms calc */
    }
  sf->chans = sp->chans;                                  /* this many separate results in results arrays */ 
  sf->op = SND_REPORT_RESULTS;                            /* report results, but no other action taken */
  sf->result_type = UFUN_ARRAY;
  return((void *)sf);
}



/* -------------------------------- FFT_FILTER -------------------------------- */
/* fft_filter: an example of an editing function (does not add to display) */
/* 
 * the expectation here is that the portion of data being smoothed is not so large as to require a work proc (in X jargon)
 */

/* the standard fft... */

#define one_pi 3.141592653589793
static void shuffle (float* rl, float* im, int n)
{
  int i,m,j;
  float tempr,tempi;
  j=0;
  for (i=0;i<n;i++)
    {
      if (j>i) {tempr = rl[j]; tempi = im[j]; rl[j] = rl[i]; im[j] = im[i]; rl[i] = tempr; im[i] = tempi;}
      m = n>>1; 
      while ((m>=2) && (j>=m)) {j -= m; m = m>>1;}
      j += m;}}

static void c_fft (float* rl, float* im, int n, int isign, int ipow)
{
  int mmax,j,pow,prev,lg,i,ii,jj;
  float wrs,wis,tempr,tempi;
  double wr,wi,theta,wtemp,wpr,wpi;
  shuffle(rl,im,n);
  mmax = 2; prev = 1; pow = n*0.5; theta = (one_pi*isign);
  for (lg=0;lg<ipow;lg++)
    {
      wpr = cos(theta); wpi = sin(theta); wr = 1.0; wi = 0.0;
      for (ii=0;ii<prev;ii++)
	{
	  wrs = (float) wr;  wis = (float) wi; i = ii; j = ii + prev;
	  for (jj=0;jj<pow;jj++)
	    {
	      tempr = wrs*rl[j] - wis*im[j]; tempi = wrs*im[j] + wis*rl[j];
	      rl[j] = rl[i] - tempr; im[j] = im[i] - tempi; rl[i] += tempr; im[i] += tempi;
	      i += mmax; j += mmax;
	    }
	  wtemp = wr; wr = (wr*wpr) - (wi*wpi); wi = (wi*wpr) + (wtemp*wpi);
	}
      pow = pow*0.5; prev = mmax; theta = theta*0.5; mmax = mmax*2;}}


void* fft_filter(void *ptr)
{
  /* args[0] = lowpass cutoff point on scale from 0.0 to 1.0 */
  /* args[1] = segment len, if given, otherwise go to window end */
  static float *data = NULL;
  static int data_size = 0;
  static float *idata;
  int i,j,k,m,n,len;
  ufun_to_snd *sf;                                        /* results passed back to Snd */
  float val,scaler;
  fsnd_fun *next_sample;
  snd_to_ufun *sp = (snd_to_ufun *)ptr;                   /* void pointer to simplify linkages */
  float cutoff,ymin,ymax,dmin,dmax;
  float fixup;
  if (sp->nargs > 0) cutoff = sp->args[0]; else cutoff = .5;
  if (sp->nargs > 1) len = pow(2.0,ceil(log(sp->args[1])/log(2.0))); else len = 0;
  if (!data)
    {
      if (len) data_size = len; else data_size = 512;     /* we'll make more room if it's needed */
      data = (float *)calloc(data_size,sizeof(float));
      idata = (float *)calloc(data_size,sizeof(float));
    }
  next_sample = sp->snd_next_sample;
  sf = (ufun_to_snd *)calloc(1,sizeof(ufun_to_snd));
  sf->results = (float **)calloc(sp->chans,sizeof(float *));
  sf->lens = (int *)calloc(sp->chans,sizeof(int));
  sf->begs = (int *)calloc(sp->chans,sizeof(int));        /* where the new data starts in each channel */
  for (i=0;i<sp->chans;i++)
    {
      j=0;
      ymin = 1.0; ymax = -1.0;
      m = sp->wbegs[i];
      while (m <= sp->wends[i])
	{
	  val = ((*next_sample)(sp->snd_ptrs[i]));
	  if (val < ymin) ymin = val;
	  if (val > ymax) ymax = val;
	  if (j == data_size)                             /* gotta make more room */
	    {
	      data_size *= 2;
	      data = (float *)realloc(data,data_size*sizeof(float));
	      idata = (float *)realloc(idata,data_size*sizeof(float));
	    }
	  data[j] = val;
	  j++;
	  if ((len) && (j >= len)) break;
	  m++;
	}
      /* now zero pad to next power of 2 */
      m = ceil(log((float)j)/log(2.0));
      k = pow(2.0,m);
      if (k > data_size)
	{
	  data_size = k;
	  data = (float *)realloc(data,data_size*sizeof(float));
	  idata = (float *)realloc(idata,data_size*sizeof(float));
	}
      c_fft(data,idata,k,1,m);
      for (n=(int)(cutoff*j);n<k;n++)
	{
	  data[n] = 0.0;
	  idata[n] = 0.0;
	}
      c_fft(data,idata,k,-1,m);
      sf->results[i] = (float *)calloc(j,sizeof(float));
      dmin = data[0]; dmax = data[0];
      for (n=1;n<j;n++) {if (data[n] < dmin) dmin = data[n]; if (data[n] > dmax) dmax = data[n];}
      scaler = (ymax - ymin)/(dmax - dmin);
      fixup = (scaler*dmin) - ymin;
      for (n=0;n<j;n++) 
	{
	  sf->results[i][n] = data[n]*scaler - fixup;
	}
      sf->lens[i] = j;
      sf->begs[i] = sp->wbegs[i];                         /* i.e. no change from whatever we were passed */
    }
  sf->chans = sp->chans;
  sf->op = SND_CHANGE_RESULTS;                            /* add our filtered data to the various edit trees */
  sf->result_type = UFUN_ARRAY;
  return((void *)sf);
}



/* -------------------------------- AUTOCORRELATION -------------------------------- */

void *autocorrelation(void *ptr)
{
  /* here we'll ignore any arguments and just take the current windowful as our basic data set */
  /* we'll compute the autocorrelation thereof, then ask Snd to display it alongside the time domain waveform */
  /* we'll also ask that each time the window changes, this function will get called to keep its display up to date */
  float *rdata,*idata;
  ufun_to_snd *sf; 
  fsnd_fun *next_sample;
  snd_to_ufun *sp = (snd_to_ufun *)ptr;
  int data_len,fft_len;
  int i,j,ipow;
  float tempr;
  next_sample = sp->snd_next_sample;
  sf = (ufun_to_snd *)calloc(1,sizeof(ufun_to_snd));
  sf->results = (float **)calloc(sp->chans,sizeof(float *));
  sf->lens = (int *)calloc(sp->chans,sizeof(int));
  sf->begs = (int *)calloc(sp->chans,sizeof(int));
  for (i=0;i<sp->chans;i++)
    {
      data_len = sp->wends[i] - sp->wbegs[i];              /* one window's worth of data */
      ipow = ceil(log(data_len)/log(2.0));
      fft_len = pow(2.0,ipow);
      rdata = (float *)calloc(fft_len,sizeof(float));      /* this could be statically allocated or whatever */
      idata = (float *)calloc(fft_len,sizeof(float));
      for (j=0;j<data_len;j++)
	{
	  rdata[j]=((*next_sample)(sp->snd_ptrs[i]));      /* read in the window's data */
	}
      c_fft(rdata,idata,fft_len,1,ipow);                   /* fft it */
      for (j=0;j<fft_len;j++)                              /* complex multiply with complex conjugate of itself */ 
	{
	  tempr = rdata[j];
	  rdata[j] = tempr*tempr + idata[j]*idata[j];
	  idata[j] = 0.0;
	}
      c_fft(rdata,idata,fft_len,-1,ipow);                  /* and inverse fft */
      sf->results[i] = (float *)calloc(data_len,sizeof(float));      
      for (j=0;j<data_len;j++)
	{
	  sf->results[i][j] = rdata[j]/(float)fft_len;     /* scale results in some plausible manner */
	}
      sf->lens[i] = data_len;
      sf->begs[i] = sp->wbegs[i];
      free(rdata);
      free(idata);
    }
  sf->chans = sp->chans;
  sf->op = SND_DISPLAY_RESULTS;
  sf->result_type = UFUN_ARRAY;                            /* currently only this works (not UFUN_FILE here) */
  sf->xlabel = (char *)calloc(9,sizeof(char));
  strcpy(sf->xlabel,"lag time");
  sf->peak_scaler = 1.0;
  sf->peaks_ok = 1;                                        /* if show_peaks option is chosen (Options Menu), include us */
  return((void *)sf);
}


/* -------------------------------- test case -------------------------------- */
#if DEBUGGING

void *times2(void *ptr)
{
  /* test case for file as return value -- grabs 1024 samples from cursor in each channel and multiplies by 2 */
  ufun_to_snd *sf;
  snd_to_ufun *sp = (snd_to_ufun *)ptr;
  char hdr[32];
  int i,j,len,fd,jump;
  float *vals;
  isnd_fun *at_eof;
  fsnd_fun *next_sample;
  fsnd_fun *previous_sample;
  at_eof = sp->snd_eof;
  next_sample = sp->snd_next_sample;
  previous_sample = sp->snd_previous_sample;
  sf = (ufun_to_snd *)calloc(1,sizeof(ufun_to_snd));
  sf->lens = (int *)calloc(sp->chans,sizeof(int));
  sf->begs = (int *)calloc(sp->chans,sizeof(int));
  sf->filename = "hiho.data";
  vals = (float *)calloc(sp->chans * 1024,sizeof(float));
  /* the data file is assumed to be a sound file of some type, so call it a NeXT float sound file */
  fd = creat("hiho.data",0666);
  (*(int *)hdr) = 0x2e736e64;
  (*(int *)(hdr+4)) = 28;
  (*(int *)(hdr+8)) = len*4*sp->chans;
  (*(int *)(hdr+12)) = 6;           /* float data, interleaved channels (below) */
  (*(int *)(hdr+16)) = 22050;       /* ignored */
  (*(int *)(hdr+20)) = sp->chans;
  write(fd,hdr,28);
  for (i=0;i<sp->chans;i++)
    {
      /* first get to cursor location */
      jump = sp->wcursors[i] - sp->wbegs[i];
      if (jump > 0) for (j=0;j<jump;j++) (*next_sample)(sp->snd_ptrs[i]);
      else for (j=jump;j>0;j--) (*previous_sample)(sp->snd_ptrs[i]);
      /* now grab 1024 samples, multiply by 2, and return results as a file */
      for (len=i,j=0;j<1024;j++,len+=sp->chans)
	{
	  vals[len] = 2 * ((*next_sample)(sp->snd_ptrs[i])); 
	}
      sf->lens[i] = 1024;
      sf->begs[i] = sp->wcursors[i];
    }
  write(fd,(char *)vals,len*4*sp->chans);
  close(fd);
  sf->chans = sp->chans;                                 
  sf->op = SND_CHANGE_RESULTS;                           
  sf->result_type = UFUN_FILE;
  return((void *)sf);
}

void *extlist(void *ptr)
{
  /* test case for list of ufuns as return value -- grab 128, delete, insert 256, regionify */
  ufun_to_snd *sf;
  ufun_list_to_snd *sflst;
  snd_to_ufun *sp = (snd_to_ufun *)ptr;
  int i,j,k;
  float *vals,*regvals;
  isnd_fun *at_eof;
  fsnd_fun *next_sample;
  fsnd_fun *previous_sample;
  at_eof = sp->snd_eof;
  next_sample = sp->snd_next_sample;

  /* make the ufun_list_to_snd structure */
  sflst = (ufun_list_to_snd *)calloc(1,sizeof(ufun_list_to_snd));
  sflst->result_type = SND_UFUN_LIST;
  sflst->len = 3;
  sflst->ops = (ufun_to_snd **)calloc(3,sizeof(ufun_to_snd *));

  /* grab 128 samples, repeat them to make 256 */
  vals = (float *)calloc(256,sizeof(float));
  for (i=0;i<128;i++)
    {
      vals[i] = ((*next_sample)(sp->snd_ptrs[0]));
      vals[i+128] = 2*vals[i];
    }

  /* now delete the original 128 */
  sf = (ufun_to_snd *)calloc(1,sizeof(ufun_to_snd));
  sf->lens = (int *)calloc(1,sizeof(int));
  sf->begs = (int *)calloc(1,sizeof(int));
  sf->chans = 1;
  sf->lens[0] = 128;
  sf->begs[0] = sp->wbegs[0];
  sf->op = SND_DELETE;
  sflst->ops[0] = sf;

  /* now insert 256 */
  sf = (ufun_to_snd *)calloc(1,sizeof(ufun_to_snd));
  sf->lens = (int *)calloc(1,sizeof(int));
  sf->begs = (int *)calloc(1,sizeof(int));
  sf->results = (float **)calloc(1,sizeof(float *));
  sf->chans = 1;
  sf->lens[0] = 256;
  sf->begs[0] = sp->wbegs[0];
  sf->results[0] = vals;
  sf->op = SND_INSERT_RESULTS;
  sf->result_type = UFUN_ARRAY;
  sflst->ops[1] = sf;

  /* now make a region of the second half of the vals array */
  regvals = (float *)calloc(128,sizeof(float));
  for (i=0;i<128;i++) regvals[i] = vals[i+128];
  sf = (ufun_to_snd *)calloc(1,sizeof(ufun_to_snd));
  sf->lens = (int *)calloc(1,sizeof(int));
  sf->begs = (int *)calloc(1,sizeof(int));
  sf->results = (float **)calloc(1,sizeof(float *));
  sf->chans = 1;
  sf->lens[0] = 128;
  sf->begs[0] = 0;
  sf->results[0] = regvals;
  sf->op = SND_REGIONIFY_RESULTS;
  sf->result_type = UFUN_ARRAY;
  sflst->ops[2] = sf;
  return((void *)sflst);
}

#endif
#endif