#include "snd.h"

#define NUM_CACHED_FFT_WINDOWS 8
#define NUM_FFT_WINDOWS 16

enum {rectangular,hanning,welch,parzen,bartlett,hamming,blackman2,blackman3,blackman4,
      exponential,riemann,kaiser,cauchy,poisson,gaussian,tukey};

static char *FFT_WINDOWS[NUM_FFT_WINDOWS] = 
     {"rectangular","hanning","welch","parzen","bartlett","hamming","blackman2","blackman3","blackman4",
      "exponential","riemann","kaiser","cauchy","poisson","gaussian","tukey"};

int default_fft_window(snd_state *ss) {if (ss) return(ss->global_fft_window); else return(blackman2);}

typedef struct {
  int s_type;
  int type;
  int inuse;
  int size;
  int zero_pad;
  float beta;
  float *window;
} fft_window_state;

typedef struct {
  int s_type;
  int n,nn,mmax,istep,m,i,size,wintype,t_type;
  double wr,wpr,wi,wpi,theta;
  int slice,inner,outer;
  void *chan;
  fft_window_state *wp;
  float *data;
  float *hwin;
  int fw_slot,hwin_size,wavelet_choice;
  int beg;
  int losamp,edit_ctr,dBing;
  int zero_pad;
  float cutoff;
  snd_state *ss;
} fft_state;

static fft_window_state *fft_windows[NUM_CACHED_FFT_WINDOWS];


static float bessi0(float x)
{ /* from Lau "A Numerical Library in C for Scientists and Engineers" CRC Press 1995 */
  float z,denominator,numerator;
  if (x == 0.0) return(1.0);
  if (fabs(x) <= 15.0) 
    {
      z=x*x;
      numerator=(z*(z*(z*(z*(z*(z*(z*(z*(z*(z*(z*(z*(z*(z*
							0.210580722890567e-22+0.380715242345326e-19)+
						     0.479440257548300e-16)+0.435125971262668e-13)+
					       0.300931127112960e-10)+0.160224679395361e-7)+
					 0.654858370096785e-5)+0.202591084143397e-2)+
				   0.463076284721000e0)+0.754337328948189e2)+
			     0.830792541809429e4)+0.571661130563785e6)+
		       0.216415572361227e8)+0.356644482244025e9)+
		 0.144048298227235e10);
      denominator=(z*(z*(z-0.307646912682801e4)+
		      0.347626332405882e7)-0.144048298227235e10);
      return(-numerator/denominator);
    } 
  else fprintf(stderr,snd_string_bessie_dumb);
  return(1.0);
}


/* -------------------------------- HANKEL TRANSFORM -------------------------------- */
/*
 * Abel transform followed by fft
 *
 * taken (with modifications) from cwplib abel.c and hankel.c by
 *   Dave Hale and Lydia Deng, Colorado School of Mines, 06/01/90 
 *   that code: Copyright (c) Colorado School of Mines, 1995. All rights reserved.
 * 
 * Original reference:
 *   Hansen, E. W., 1985, Fast Hankel transform algorithm:  IEEE Trans. on
 *   Acoustics, Speech and Signal Processing, v. ASSP-33, n. 3, p. 666-671.
 */

#define NSE 9
static float h[NSE] = {1.000000000000000000,0.610926299405048390,0.895089852938535935,1.34082948787002865,2.02532848558443890,
		       3.18110895533701843,5.90898360396353794,77.6000213494180286,528.221800846070892};    
static float lambda[NSE] = {0.000000000000000000,-2.08424632126539366,-5.78928630565552371,-14.6268676854951032,
			    -35.0617158334443104,-83.3258406398958158,-210.358805421311445,-6673.64911325382036,-34897.7050244132261};

typedef struct abeltStruct {int n; float **a,**b0,**b1;} abelt;
static abelt *at = NULL;

static void make_abel_transformer(int n)
{
  int i,j,nse=NSE;
  float **a,**b0,**b1,fi,hj,lambdaj,scale,temp;
  if ((!at) || (at->n != n))
    {
      if (at) {for (i=0;i<at->n;i++) {free(at->a[i]); free(at->b0[i]); free(at->b1[i]);} free(at->a); free(at->b0); free(at->b1);}
      else at = (abelt *)calloc(1,sizeof(abelt));
      a = (float **)calloc(n,sizeof(float *));
      b0 = (float **)calloc(n,sizeof(float *));
      b1 = (float **)calloc(n,sizeof(float *));
      for (i=0;i<n;i++) 
	{
	  a[i] = (float *)calloc(nse,sizeof(float));
	  b0[i] = (float *)calloc(nse,sizeof(float));
	  b1[i] = (float *)calloc(nse,sizeof(float));
	}
      for (i=1; i<n; ++i) 
	{
	  fi = (float)i+1.0;
	  for (j=0; j<nse; ++j) 
	    {
	      hj = h[j];
	      lambdaj = lambda[j];
	      a[i][j] = temp = pow(fi/(fi-1.0),lambdaj);
	      temp *= fi/(fi-1.0);
	      scale = 2.0*hj*(fi-1.0) / ((lambdaj+1.0)*(lambdaj+2.0));				
	      b0[i][j] = scale * (fi-1.0+(lambdaj+2.0-fi)*temp);
	      b1[i][j] = -scale * (lambdaj+1.0+fi-fi*temp);
	    }
	}
      at->n = n;
      at->a = a;
      at->b0 = b0;
      at->b1 = b1;
    }
}

static void abel (float *f, float *g)
{
  int i,j,n,nse=NSE;
  float **a,**b0,**b1,xi[NSE],sum,fi,fip1;
  n = at->n;
  a = at->a;
  b0 = at->b0;
  b1 = at->b1;
  fi = f[n-1];
  g[0] = 0.5*f[0]+fi;
  for (j=0,sum=0.0; j<nse; ++j)
    {
      xi[j] = b1[n-1][j]*fi;
      sum += xi[j];
    }
  g[n-1] = sum;
  for (i=n-2; i>0; --i) 
    {
      fip1 = fi;
      fi = f[i];
      g[0] += fi;
      for (j=0,sum=0.0; j<nse; ++j) 
	{
	  xi[j] = a[i][j]*xi[j] + b0[i][j]*fip1 + b1[i][j]*fi;
	  sum += xi[j];
	}
      g[i] = sum;
    }
  g[0] *= 2.0;
}

/*
;;; test cases:

(with-sound (:output "j") ;sum of bessel funcs
  (loop for i from 0 below 1024 and r from 0.0 by (/ 1.0 1024) do 
    (outa i (* .5 (+ (bes-jn 0 (* 500 pi r)) 
                     (bes-jn 0 (* 150 pi r)))))))

(with-sound (:output "j") ;jinc
  (loop for i from 0 below 1024 and r from 0.0 by (/ 1.0 1024) do
    (if (= r 0.0)
	(outa i .999)
      (outa i (/ (bes-j1 (* 40 pi r)) (* r 20 pi))))))
*/


/* -------------------------------- WAVELET TRANSFORM -------------------------------- */
/* 
 * taken from wavelets.cl in clm which is taken from
 * "Numerical Recipes in C" second edition, M. J. Shensa Naval Ocean Systems Center,
 *  Wickerhauser "Adapted Wavelet Analysis", and the UBC Imager Wavelet Package by Bob Lewis
 */

static int data1_size = 0;
static float *data1 = NULL;

static void pwt(float *data, int n, int isign, float *cc, float *cr, int ncof)
{
  int n1,nh,joff;
  int ii,i,jr,nj,k,iinh;
  float val;
  if (n<4) return;
  if (n > data1_size)
    {
      if (data1) free(data1);
      data1 = (float *)calloc(n,sizeof(float));
      data1_size = n;
    }
  else 
    {
      for (i=0;i<n;i++) data1[i]=0.0;
    }
  n1 = n-1;
  nh = n>>1;
  joff = ncof*(n-0.5);
  if (isign>=0)
    {
      for (ii=0,i=1,iinh=nh;i<=n;i+=2,ii++,iinh++)
	{
	  nj = i+joff;
	  for (k=1;k<=ncof;k++)
	    {
	      jr = n1 & (nj+k);
	      val = cc[k-1] * data[jr];
	      data1[ii] += val;
	      data1[iinh] += val;
	    }
	}
    }
  else
    {
      for (ii=0,iinh=nh,i=1;i<=n;i+=2,ii++,iinh++)
	{
	  nj = i+joff;
	  for (k=1;k<=ncof;k++)
	    {
	      jr = n1 & (nj+k);
	      data1[jr] += data[ii] * cc[k-1];
	      data1[jr] += data[iinh] * cc[k-1];
	    }
	}
    }
  for (i=0;i<n;i++) data[i] = data1[i];
}

static void wavelet_transform(float *data, int n, int isign, float *cc, int cc_size)
{
  float *cr;
  int i,j,sig,nn;
  cr = (float *)calloc(cc_size,sizeof(float));
  sig = -1;
  for (i=0,j=cc_size-1;i<cc_size;i++,j--) 
    {
      cr[j] = sig*cc[i];
      sig = -sig;
    }
  if (n>4)
    {
      if (isign>=0)
	{
	  for (nn=n;nn>=4;nn>>=1) pwt(data,nn,isign,cc,cr,cc_size);
	}
      else
	{
	  for (nn=4;nn<=n;nn<<=1) pwt(data,nn,isign,cc,cr,cc_size);
	}
    }
  free(cr);
}

static float daub4[4] = {0.4829629131445341, 0.8365163037378079, 0.2241438680420134, -0.1294095225512604};
static float daub6[6] = {0.332670552950, 0.806891509311, 0.459877502118, -0.135011020010, -0.085441273882, 0.035226291886};
static float daub8[8] = {0.230377813309, 0.714846570553, 0.630880767930, -0.027983769417, -0.187034811719, 0.030841381836, 
			 0.032883011667, -0.010597401785};
static float daub10[10] = {0.160102397974, 0.603829269797, 0.724308528438, 0.138428145901, -0.242294887066, -0.032244869585,
			   0.077571493840, -0.006241490213, -0.012580751999, 0.003335725285};
static float daub12[12] = {0.111540743350, 0.494623890398, 0.751133908021, 0.315250351709, -0.226264693965, -0.129766867567,
			   0.097501605587, 0.027522865530, -0.031582039317, 0.000553842201, 0.004777257511, -0.001077301085};
static float daub14[14] = {0.077852054085, 0.396539319482, 0.729132090846, 0.469782287405, -0.143906003929, -0.224036184994,
			   0.071309219267, 0.080612609151, -0.038029936935, -0.016574541631, 0.012550998556, 0.000429577973,
			   -0.001801640704, 0.000353713800};
static float daub16[16] = {0.054415842243, 0.312871590914, 0.675630736297, 0.585354683654, -0.015829105256, -0.284015542962,
			   0.000472484574, 0.128747426620, -0.017369301002, -0.044088253931, 0.013981027917, 0.008746094047,
			   -0.004870352993, -0.000391740373, 0.000675449406, -0.000117476784};
static float daub18[18] = {0.038077947364, 0.243834674613, 0.604823123690, 0.657288078051, 0.133197385825, -0.293273783279,
			   -0.096840783223, 0.148540749338, 0.030725681479, -0.067632829061, 0.000250947115, 0.022361662124,
			   -0.004723204758, -0.004281503682, 0.001847646883, 0.000230385764, -0.000251963189, 0.000039347320};
static float daub20[20] = {0.026670057901, 0.188176800077, 0.527201188931, 0.688459039453, 0.281172343661, -0.249846424327,
			   -0.195946274377, 0.127369340336, 0.093057364604, -0.071394147166, -0.029457536822, 0.033212674059,
			   0.003606553567, -0.010733175483, 0.001395351747, 0.001992405295, -0.000685856695, -0.000116466855,
			   0.000093588670, -0.000013264203};
static float Battle_Lemarie[24] = {-0.0028284274, -0.004242641, 0.008485282, 0.008485282, -0.018384777, -0.016970564, 0.042426407, 
				   0.032526914, -0.11030867, -0.049497478, 0.4341636, 0.7665038, 0.4341636, -0.049497478, -0.11030867, 
				   0.032526914, 0.042426407, -0.016970564, -0.018384777, 0.008485282, 0.008485282, -0.004242641, -0.0028284274, 0.0};
static float Burt_Adelson[6] = {-0.07071068, 0.3535534, 0.8485282, 0.3535534, -0.07071068, 0.0};
static float Beylkin[18] = {0.099305765374353, 0.424215360812961, 0.699825214056600, 0.449718251149468,
			    -.110927598348234, -.264497231446384, 0.026900308803690, 0.155538731877093,
			    -.017520746266529, -.088543630622924, 0.019679866044322, 0.042916387274192,
			    -.017460408696028, -.014365807968852, 0.010040411844631, .0014842347824723,
			    -.002736031626258, .0006404853285212};
static float coif2[6] = {0.038580775, -0.12696913, -0.07716155, 0.6074917, 0.74568766, 0.2265843};
static float coif4[12] = {0.0011945726958388, -0.01284557955324, 0.024804330519353, 0.050023519962135, -0.15535722285996,
			  -0.071638282295294, 0.57046500145033, 0.75033630585287, 0.28061165190244, -0.0074103835186718,
			  -0.014611552521451, -0.0013587990591632};
static float coif6[18] = {-0.0016918510194918, -0.00348787621998426, 0.019191160680044, 0.021671094636352, -0.098507213321468,
			  -0.056997424478478, 0.45678712217269, 0.78931940900416, 0.38055713085151, -0.070438748794943,
			  -0.056514193868065, 0.036409962612716, 0.0087601307091635, -0.011194759273835, -0.0019213354141368,
			  0.0020413809772660, 0.00044583039753204, -0.00021625727664696};
static float sym2[5] = {-0.1767767, 0.3535534, 1.0606602, 0.3535534, -0.1767767};
static float sym3[4]= {0.1767767, 0.5303301, 0.5303301, 0.1767767};
static float sym4[10] = {0.033145633, -0.066291265, -0.1767767, 0.4198447, 0.994369, 0.4198447, -0.1767767, -0.066291265, 0.033145633, 0.0};
static float sym5[8] = {0.066291265, -0.19887379, -0.15467963, 0.994369, 0.994369, -0.15467963, -0.19887379, 0.066291265};
static float sym6[16] = {-0.0030210863, -0.009063259, -0.016831767, 0.07466399, 0.03133298, -0.30115914, -0.026499243, 
			 0.9516422, 0.9516422, -0.026499243, -0.30115914, 0.03133298, 0.07466399, -0.016831767, -0.009063259, -0.0030210863};

static char *wavelet_names[]={"daub4","daub6","daub8","daub10","daub12","daub14","daub16","daub18","daub20",
			      "battle_lemarie","burt_adelson","beylkin","coif2","coif4","coif6",
			      "sym2","sym3","sym4","sym5","sym6"};
static int wavelet_sizes[]={4,6,8,10,12,14,16,18,20,
			    24,6,18,6,12,18,
			    5,4,10,8,16};
static float *wavelet_data[]={daub4,daub6,daub8,daub10,daub12,daub14,daub16,daub18,daub20,
			      Battle_Lemarie,Burt_Adelson,Beylkin,coif2,coif4,coif6,
			      sym2,sym3,sym4,sym5,sym6};


/* tested with pulse train and oboe.snd */



/* -------------------------------- CHEBYSHEV TRANSFORM -------------------------------- */
/*
 * translated from clm's bessel.lisp which has all the references
 */

static void chebyshev_polynomials(double x, double *f, int n)
{
  int k;
  double r,s,h;
  r = x;
  s = 1.0;
  h = 0.0;
  f[0] = 1.0;
  f[1] = r;
  for (k=1;k<n;k++)
    {
      h = r;
      r = 2*r*x - s;
      f[k+1] = r;
      s = h;
    }
}

static void chebyshev_transform(float *data, int n)
{
  double *An,*Pn;
  int i,k;
  double x,rate,wx,wfx,ln2;
  rate = 2.0/(float)n;
  ln2 = log(n)/log(2);
  An = (double *)calloc(n,sizeof(double));
  Pn = (double *)calloc(n,sizeof(double));
  for (k=0,x=-1.0;k<n;k++,x+=rate)
    {
      if ((x == 1.0) || (x == -1.0))
	wx = ln2;
      else wx = 1.0/sqrt(1.0 - x*x);
      wfx = wx * data[k];
      chebyshev_polynomials(x,Pn,n-1);
      for (i=0;i<n;i++)	An[i] += (wfx * Pn[i]);
    }
  for (i=0;i<n;i++) data[i] = An[i];
  free(An);
  free(Pn);
}

/* 
;;; test case (bessel.lisp):

   (with-sound (:output "wave") 
     (let ((a (make-array 8 :initial-contents '(0 0 0 .25 0 0 0 0.5))))
       (loop for i from 0 below 256 do (outa i (* .5 (chebyshev-polynomial a (- (/ i 128) 1.0)))))))
*/


/* -------------------------------- LEGENDRE TRANSFORM -------------------------------- */

static void legendre_polynomials(double x, double *f, int n)
{
  int k;
  double r,s,h;
  r = x;
  s = 1.0;
  h = 0.0;
  f[0] = 1.0;
  f[1] = x;
  for (k=1;k<n;k++)
    {
      h = r;
      r = (r*x*(1+2*k) - s*k) / (float)(k+1);
      f[k+1] = r;
      s = h;
    }
}

static void legendre_transform(float *data, int n)
{
  int i,k;
  double x,rate;
  double *An,*Pn;
  rate = 2.0/(float)n;
  An = (double *)calloc(n,sizeof(double));
  Pn = (double *)calloc(n,sizeof(double));
  for (k=0,x=-1.0;k<n;k++,x+=rate)
    {
      legendre_polynomials(x,Pn,n-1);
      for (i=0;i<n;i++) An[i] += data[k] * Pn[i];
    }
  for (i=0;i<n;i++) data[i] = An[i];
  free(An);
  free(Pn);
}

/*
(with-sound (:output "wave") 
     (let ((a (make-array 8 :initial-contents '(0 0 0 .25 0 .5 0 0))))
       (loop for i from 0 below 500 and j from -1.0 by .004 do (outa i (* .5 (legendre-polynomial a j))))))
*/


/* -------------------------------- FFT DATA WINDOW -------------------------------- */

static int make_fft_window_1(float *window, int size, int type, float beta)
{
  /* mostly taken from
   *    Fredric J. Harris, "On the Use of Windows for Harmonic Analysis with the
   *    Discrete Fourier Transform," Proceedings of the IEEE, Vol. 66, No. 1,
   *    January 1978.
   *    Albert H. Nuttall, "Some Windows with Very Good Sidelobe Behaviour", 
   *    IEEE Transactions of Acoustics, Speech, and Signal Processing, Vol. ASSP-29,
   *    No. 1, February 1981, pp 84-91
   *
   * JOS had slightly different numbers for the Blackman-Harris windows.
   */
  int i,j,midn,midp1,midm1;
  float freq,rate,sr,angle,expn,expsum,I0beta,cx;
  I0beta = 0.0;
  midn = size >> 1;
  midp1 = (size+1)/2;
  midm1 = (size-1)/2;
  freq = two_pi/(float)size;
  rate = 1.0/(float)midn;
  angle = 0.0;
  expn = log(2)/(float)midn+1.0;
  expsum = 1.0;
  if (type == kaiser) I0beta = bessi0(beta); /* Harris multiplies beta by pi */
  switch (type)
    {
    case rectangular:
      for (i=0;i<size;i++) {window[i]=1.0;}
      break; 
    case hanning: /* Hann would be more accurate */
      for (i=0,j=size-1,angle=0.0;i<=midn;i++,j--,angle+=freq) {window[j]=(window[i]=0.5-0.5*cos(angle));}
      break; 
    case welch:
      for (i=0,j=size-1;i<=midn;i++,j--) {window[j]=(window[i]=1.0-sqr((float)(i-midm1)/(float)midp1));}
      break; 
    case parzen:
      for (i=0,j=size-1;i<=midn;i++,j--) {window[j]=(window[i]=1.0-fabs((float)(i-midm1)/(float)midp1));}
      break; 
    case bartlett:
      for (i=0,j=size-1,angle=0.0;i<=midn;i++,j--,angle+=rate) {window[j]=(window[i]=angle);}
      break; 
    case hamming:
      for (i=0,j=size-1,angle=0.0;i<=midn;i++,j--,angle+=freq) {window[j]=(window[i]=0.54-0.46*cos(angle));} 
      break; 
    case blackman2: /* using Chebyshev polynomial equivalents here */
      for (i=0,j=size-1,angle=0.0;i<=midn;i++,j--,angle+=freq) 
	{              /* (+ 0.42323 (* -0.49755 (cos a)) (* 0.07922 (cos (* a 2)))) */
	  cx = cos(angle);
	  window[j]=(window[i]=(.34401+(cx*(-.49755+(cx*.15844)))));
	}
      break; 
    case blackman3:
      for (i=0,j=size-1,angle=0.0;i<=midn;i++,j--,angle+=freq) 
	{              /* (+ 0.35875 (* -0.48829 (cos a)) (* 0.14128 (cos (* a 2))) (* -0.01168 (cos (* a 3)))) */
	  cx = cos(angle);
	  window[j]=(window[i]=(.21747+(cx*(-.45325+(cx*(.28256-(cx*.04672)))))));
	}
      break; 
    case blackman4:
      for (i=0,j=size-1,angle=0.0;i<=midn;i++,j--,angle+=freq) 
	{             /* (+ 0.287333 (* -0.44716 (cos a)) (* 0.20844 (cos (* a 2))) (* -0.05190 (cos (* a 3))) (* 0.005149 (cos (* a 4)))) */
	  cx = cos(angle);
	  window[j]=(window[i]=(.084037+(cx*(-.29145+(cx*(.375696+(cx*(-.20762+(cx*.041194)))))))));
	}
      break; 
    case exponential:
      for (i=0,j=size-1;i<=midn;i++,j--) {window[j]=(window[i]=expsum-1.0); expsum *= expn;}
      break;
    case kaiser:
      for (i=0,j=size-1,angle=1.0;i<=midn;i++,j--,angle-=rate) {window[j]=(window[i]=bessi0(beta*sqrt(1.0-sqr(angle)))/I0beta);}
      break;
    case cauchy:
      for (i=0,j=size-1,angle=1.0;i<=midn;i++,j--,angle-=rate) {window[j]=(window[i]=1.0/(1.0+sqr(beta*angle)));}
      break;
    case poisson:
      for (i=0,j=size-1,angle=1.0;i<=midn;i++,j--,angle-=rate) {window[j]=(window[i]=exp((-beta)*angle));}
      break;
    case riemann:
      sr = two_pi/(float)size;
      for (i=0,j=size-1;i<=midn;i++,j--) 
	{
	  if (i == midn) window[j]=(window[i]=1.0);
	  else window[j]=(window[i]=sin(sr*(midn-i))/(sr*(midn-i)));
	}
      break;
    case gaussian:
      for (i=0,j=size-1,angle=1.0;i<=midn;i++,j--,angle-=rate) {window[j]=(window[i]=exp(-.5*sqr(beta*angle)));}
      break;
    case tukey:
      cx = midn*(1.0-beta);
      for (i=0,j=size-1;i<=midn;i++,j--) 
	{
	  if (i >= cx) window[j]=(window[i]=1.0);
	  else window[j]=(window[i]=.5*(1.0-cos(one_pi*i/cx)));
	}
      break;
    default: snd_error("unknown window type"); break;
    }
  for (i=0;i<size;i++) {window[i] *= clm_sndflt;}
  return(1);
}

static int compare_peaks(const void *pk1, const void *pk2)
{
  if (((fft_peak *)pk1)->freq > ((fft_peak *)pk2)->freq) return(1);
  else if (((fft_peak *)pk1)->freq == ((fft_peak *)pk2)->freq) return(0);
  return(-1);
}

int find_and_sort_peaks(float *buf, fft_peak *found, int num_peaks, int size)
{ /* in the fft peak finder below we assume data between 0 and 1 */
  int i,j,pks,minpk;
  float minval,la,ra,ca;
  float peaks[MAX_NUM_PEAKS];
  int inds[MAX_NUM_PEAKS];
  pks = 0;
  la = 0.0;
  ca = 0.0;
  ra = 0.0;
  minval = 0.00001;
  for (i=0;i<size;i++)
    {
      la = ca;
      ca = ra;
      ra = buf[i];
      if ((ca > minval) && (ca > ra) && (ca > la))
	{
	  if (pks < num_peaks)
	    {
	      inds[pks] = i-1;
	      peaks[pks++] = ca;
	    }
	  else
	    {
	      minval = peaks[0];
	      minpk = 0;
	      for (j=1;j<num_peaks;j++)
		{
		  if (peaks[j] < minval) 
		    {
		      minval = peaks[j];
		      minpk = j;
		    }
		}
	      if (ca > minval)
		{
		  inds[minpk] = i-1;
		  peaks[minpk] = ca;
		}
	    }
	}
    }
  for (i=0;i<pks;i++)
    {
      j = inds[i];
      ca = buf[j];
      found[i].amp = buf[j];
      found[i].freq = j;
    }
  qsort((void *)found,pks,sizeof(fft_peak),compare_peaks);
  return(pks);
}

int find_and_sort_fft_peaks(float *buf, fft_peak *found, int num_peaks, int fftsize2, int srate, float samps_per_pixel)
{
  /* we want to reflect the graph as displayed, so each "bin" is samps_per_pixel wide */
  int i,j,k,pks,minpk,hop,pkj,oldpkj;
  float minval,la,ra,ca,logca,logra,logla,offset,fscl,ascl;
  float peaks[MAX_NUM_PEAKS];
  int inds[MAX_NUM_PEAKS];
  fscl = (float)srate/(float)fftsize2;
  hop = (int)(samps_per_pixel+0.5);
  if (hop < 1) hop=1;
  if (num_peaks > MAX_NUM_PEAKS) num_peaks = MAX_NUM_PEAKS;
  pks = 0;
  la = 0.0;
  ca = 0.0;
  ra = 0.0;
  minval = (float)fftsize2/100000.0;
  ascl = 0.0;
  pkj = 0;
  for (i=0;i<fftsize2;i+=hop)
    {
      la = ca;
      ca = ra;
      oldpkj = pkj;
      ra = 0.0;
      for (k=0;k<hop;k++) 
	{
	  if (buf[i+k]>ra) {pkj = i+k; ra=buf[pkj];} /* reflect user's view of the graph */
	}
      if ((ca > minval) && (ca > ra) && (ca > la))
	{
          if (ascl < ca) ascl = ca;
	  if (pks < num_peaks)
	    {
	      inds[pks] = oldpkj;
	      peaks[pks] = ca;
	      pks++;
	    }
	  else
	    {
	      minval = peaks[0];
	      minpk = 0;
	      for (j=1;j<num_peaks;j++)
		{
		  if (peaks[j] < minval) 
		    {
		      minval = peaks[j];
		      minpk = j;
		    }
		}
	      if (ca > minval)
		{
		  inds[minpk] = oldpkj;
		  peaks[minpk] = ca;
		}
	    }
	}
    }
  /* now we have the peaks; turn these into interpolated peaks/amps, and sort */
  for (i=0;i<pks;i++)
    {
      j = inds[i];
      ca = buf[j]/ascl;
      la = buf[j-1]/ascl; 
      ra = buf[j+1]/ascl; 
      if (la<.00001) la=.00001;
      if (ra<.00001) ra=.00001;
      logla = log10(la);
      logca = log10(ca);
      logra = log10(ra);
      offset = (0.5*(logla-logra))/(logla+logra-2*logca); /* this assumes amps<1.0 (from XJS sms code) */
      found[i].amp = pow(10.0,logca-0.25*offset*(logla-logra));
      found[i].freq = fscl*(j+offset);
      if (found[i].freq < 0.0) found[i].freq = 0.0;
    }
  qsort((void *)found,pks,sizeof(fft_peak),compare_peaks);
  return(pks);
}


static int shuffle_fft_state(fft_state *fs)
{
  int i,j,m;
  float tempr;
  float *data;
  data = fs->data;
  if (!data) snd_error("no fft data!");
  fs->n=(fs->nn*2);
  j=1;
  for (i=1;i<fs->n;i+=2) 
    {
      if (j > i) 
	{
	  tempr=data[j];
	  data[j]=data[i];
	  data[i]=tempr;
	  tempr=data[j+1];
	  data[j+1]=data[i+1];
	  data[i+1]=tempr;
	}
      m=fs->n >> 1;
      while (m >= 2 && j > m) 
	{
	  j -= m;
	  m >>= 1;
	}
      j += m;
    }
  fs->mmax=2;
  fs->outer = 1;
  return(1);
}

static int snd_fft(fft_state *fs)
{
  double wtemp;
  float tempr,tempi;
  float *data;
  int k,j;
  data = fs->data;
  if (fs->n == fs->mmax) return(1);
  if (fs->outer)
    {
      fs->outer = 0;
      fs->istep=2*fs->mmax;
      fs->theta=6.28318530717959/fs->mmax;
      wtemp=sin(0.5*fs->theta);
      fs->wpr = -2.0*wtemp*wtemp;
      fs->wpi=sin(fs->theta);
      fs->wr=1.0;
      fs->wi=0.0;
      fs->m = 1;
      fs->inner = 1;
    }
  if (fs->inner)
    {
      fs->i = fs->m;
      fs->inner = 0;
    }
  k=0;
LOOP:
  j=fs->i+fs->mmax;
  tempr=fs->wr*data[j]-fs->wi*data[j+1];
  tempi=fs->wr*data[j+1]+fs->wi*data[j];
  data[j]=data[fs->i]-tempr;
  data[j+1]=data[fs->i+1]-tempi;
  data[fs->i] += tempr;
  data[fs->i+1] += tempi;
  k++;
  fs->i += fs->istep;
  if (fs->i > fs->n)
    {
      fs->inner = 1;
      fs->wr=(wtemp=fs->wr)*fs->wpr-fs->wi*fs->wpi+fs->wr;
      fs->wi=fs->wi*fs->wpr+wtemp*fs->wpi+fs->wi;
      fs->m += 2;
      if (fs->m >= fs->mmax)
	{
	  fs->outer = 1;
	  fs->mmax = fs->istep;
	}
    }
  else 
    if (k < 100) goto LOOP;
  return(0);
}

static int snd_fft_cleanup(fft_state *fs)
{
  double wtemp;
  int n2p3,i,i1,i2,i3,i4;
  float h1r,h1i,h2r,h2i;
  float *data;
  data = fs->data;
  fs->theta=3.141592653589793/((double)fs->nn);
  wtemp=sin(0.5*fs->theta);
  fs->wpr = -2.0*wtemp*wtemp;
  fs->wpi=sin(fs->theta);
  fs->wr=1.0+fs->wpr;
  fs->wi=fs->wpi;
  n2p3=2*fs->nn+3;
  for (i=2;i<=(fs->nn/2);i++) 
    {
      i1=i+i-1;
      i2=i1+1;
      i3=n2p3-i2;
      i4=1+i3;
      h1r=0.5*(data[i1]+data[i3]);
      h1i=0.5*(data[i2]-data[i4]);
      h2r=0.5*(data[i2]+data[i4]);
      h2i=(-0.5)*(data[i1]-data[i3]);
      data[i1]=h1r+fs->wr*h2r-fs->wi*h2i;
      data[i2]=h1i+fs->wr*h2i+fs->wi*h2r;
      data[i3]=h1r-fs->wr*h2r+fs->wi*h2i;
      data[i4] = -h1i+fs->wr*h2i+fs->wi*h2r;
      fs->wr=(wtemp=fs->wr)*fs->wpr-fs->wi*fs->wpi+fs->wr;
      fs->wi=fs->wi*fs->wpr+wtemp*fs->wpi+fs->wi;
    }
  h1r=data[1];
  data[1] = h1r+data[2];
  data[2] = h1r-data[2];
  return(1);
}

static int snd_fft_to_spectrum (fft_state *fs)
{
  int i,j;
  float *fft_data;
  snd_state *ss;
  float nan_check; /* sqrt is screwing up -- floating underflow */
  fft_data = fs->data;
  ss = fs->ss;
  /* realft returns the srate/2 value as the imaginary part of the DC component "as a convenience" */
#ifdef LINUX 
  /* this does not need to be called -- it just forces some portion of the math library to be loaded */
  /* without it, we occasionally get NaNs where we should not */
  if (isnan(fft_data[1])) fft_data[1]=0.0;
#endif
  if (ss->transform_type == hankel)
    {
      nan_check = 0.0;
      for (j=0;j<fs->size;j+=2)
	{
	  if (fft_data[j]>nan_check) 
	    nan_check = fft_data[j];
	  else
	    if ((-fft_data[j] ) > nan_check)
	      nan_check = (-fft_data[j]);
	}
      for (i=0,j=0;i<fs->size;j++,i+=2)
	{
	  fft_data[j] = fft_data[i] * nan_check;
	}
    }
  else
    {
      if (fft_data[1] < 0.0001) fft_data[0] = 0.0; else fft_data[0] = fft_data[1];
      if (fft_data[2] < 0.0001) fft_data[fs->size-1] = 0.0; else fft_data[fs->size-1] = fft_data[2];
#ifdef LINUX
      for (i=3,j=1;i<fs->size-3;i+=2,j++) /* changed 13-Jan-98 to fix srate/2 problem */
	{
	  nan_check = fft_data[i]*fft_data[i]+fft_data[i+1]*fft_data[i+1];
	  if (nan_check < 0.0001) fft_data[j] = 0.0; else fft_data[j] = sqrt(nan_check);
	}
#else
      for (i=3,j=1;i<fs->size-3;i+=2,j++)
	{
	  fft_data[j] = sqrt(fft_data[i]*fft_data[i]+fft_data[i+1]*fft_data[i+1]);
	}
#endif
    }
  return(1);
}



int fftwincpy(char *s1, int type)
{
  strcpy(s1,FFT_WINDOWS[type]);
  return(0);
}

int fft_window_needs_beta (int window)
{
  return(window>=kaiser);
}

static float current_betas[5] = {7.0,4.0,3.0,3.0,0.5};
static float beta_maxes[5] = {20.0,10.0,10.0,10.0,1.0};

float fft_window_beta(int window)
{
  if (fft_window_needs_beta(window)) 
    return(current_betas[window-kaiser]);
  else return(0.0);
}

float fft_window_beta_max(int window)
{
  if (fft_window_needs_beta(window)) 
    return(beta_maxes[window-kaiser]);
  else return(0.0);
}

void fft_set_current_window_beta(float beta, int window)
{
  if (fft_window_needs_beta(window)) 
    current_betas[window-kaiser]=beta;
}

static int make_fft_window(fft_state *fs)
{
  /* build fft window, taking int->float transition into account */
  float *window;
  fft_window_state *wp;
  int toploc;
  snd_state *ss;
  ss = fs->ss;
  switch (ss->transform_type)
    {
    case fourier:
      wp = (fft_window_state *)(fs->wp);
      window = wp->window;
      toploc = fs->size / (float)(1 + fs->zero_pad);
      return(make_fft_window_1(window,toploc,wp->type,wp->beta));
      break;
    case hankel: 
      if (fs->size != fs->hwin_size)
	{
	  fs->hwin_size = fs->size;
	  if (fs->hwin) free(fs->hwin);
	  fs->hwin = (float *)calloc(fs->size,sizeof(float));
	}
      break;
    case wavelet: case chebyshev: case legendre:
      break;
    }
  return(1);
}

static void free_fft_window(int i)
{
  if (fft_windows[i])
    {
      if (((fft_window_state *)fft_windows[i])->window) free(((fft_window_state *)fft_windows[i])->window);
      free(fft_windows[i]);
      fft_windows[i]=NULL;
    }
}

static fft_state *free_fft_state(fft_state *fs)
{
  if (fs) 
    {
      if (fs->fw_slot == -1) free(fs->wp); /* free window only if it's not in the cache */
      if (fs->hwin) free(fs->hwin);
      free(fs); 
    }
  return(NULL);
}

static void decrement_fft_window_use(fft_state *fs)
{
  if (fs->fw_slot != -1)
    fft_windows[fs->fw_slot]->inuse--;
}

static int set_up_fft_window(fft_state *fs)
{
  int i,empty,ok,unused;
  fft_window_state *wp;
  snd_state *ss;
  ss = fs->ss;
  if (ss->transform_type != fourier) return(1);
  /* first look to see if it already exists */
  empty = -1;
  ok = -1;
  unused = -1;
  for (i=0;i<NUM_CACHED_FFT_WINDOWS;i++)
    {
      wp = fft_windows[i];
      if (!wp) 
	{
	  if (empty == -1) empty = i;
	}
      else
	{
	  if (!wp->inuse) unused = i;
	  if (wp->size == fs->size) 
	    {

	      if ((wp->type == fs->wintype) && (wp->beta == fft_window_beta(wp->type)) && (wp->zero_pad == fs->zero_pad))
		{
		  fs->wp = wp;
		  fs->fw_slot = i;
		  wp->inuse++;
		  return(2);  /* skip making window */
		}
	      if (ok == -1) ok = i;
	    }
	}
    }
  if (empty == -1) empty = ok;
  if (empty == -1) empty = unused;
  if (empty == -1) 
    {
      wp = (fft_window_state *)calloc(1,sizeof(fft_window_state));
      fs->fw_slot = -1;
    }
  else
    {
      if (empty == unused) free_fft_window(empty);
      if (!fft_windows[empty])
	{
	  fft_windows[empty] = (fft_window_state *)calloc(1,sizeof(fft_window_state));
	}
      wp = fft_windows[empty];
      fs->fw_slot = empty;
    }
  fs->wp = wp;
  wp->s_type = make_snd_pointer_type(FFT_WINDOW_STATE);
  wp->size = fs->size;
  wp->zero_pad = fs->zero_pad;
  wp->type = fs->wintype;
  wp->beta = fft_window_beta(fs->wintype);
  wp->inuse++;
  if (!wp->window) wp->window = (float *)calloc(fs->size,sizeof(float));
  return(1);
}

 
/*-------------------------------- FFT_INFO -------------------------------- */

static fft_info *make_fft_info(int size, int window, float beta)
{
  fft_info *fp;
  fp = (fft_info *)calloc(1,sizeof(fft_info));
  fp->s_type = make_snd_pointer_type(FFT_INFO);
  fp->size = size;
  fp->window = window;
  fp->beta = beta;
  fp->ok = 1;
  fp->data = (float *)calloc(size+1,sizeof(float)); /* +1 for complex storage or starts at 1 or something */
  return(fp);
}

fft_info *free_fft_info(fft_info *fp)
{
  fp->chan = NULL;
  if (fp->data) free(fp->data);
  if (fp->axis) free_axis_info(fp->axis);
  free(fp);
  return(NULL);
}



/* -------------------------------- FFT_WINDOW_STATE, FFT_STATE -------------------------------- 
 *
 * FFT as work proc, using the "realft" version, and splitting out as many steps as possible.
 * the basic fft here is taken from Xavier Serra's SMS program (CLM's is based on two 0-based arrays,
 * which in this case is not so useful).  Number of splits depends on the FFT size.
 */

static int snd_fft_set_up(fft_state *fs)
{
  /* allocate arrays if needed */
  fft_info *fp;
  chan_info *cp;
  snd_state *ss;
  axis_info *fap;
  float max_freq,max_val,min_val;
  char *xlabel;
  cp = (chan_info *)(fs->chan);
  ss = cp->state;
  fp = cp->fft;
  if (!fp)                              /* associated channel hasn't done any ffts yet, so there's no struct */
    {
      cp->fft = make_fft_info(fs->size,fs->wintype,0.0);
      fp = cp->fft;
    }
  else
    {
      if ((!fp->ok) || (!fp->data) || (fs->size > fp->size))
	{
	  fp->size = fs->size;
	  if (fp->data) free(fp->data);
	  fp->data = (float *)calloc(fp->size+1,sizeof(float));
	  fp->ok = 1;
	}
    }
  fp->current_size = fs->size; /* protect against parallel size change via fft size menu */
  switch (ss->transform_type)
    {
    case fourier: 
      if (ss->logxing)
	{
	  xlabel = snd_string_log_freq;
	  max_freq = 1.0;
	}
      else
	{
	  xlabel = snd_string_frequency;
	  max_freq = ((float)(snd_SRATE(cp)) * 0.5 * ss->sonogram_cutoff);
	}
      break;
    case wavelet: 
      xlabel = wavelet_names[ss->wavelet_choice]; 
      max_freq = fs->size * ss->sonogram_cutoff;
      break;
    case hankel: 
      xlabel = "Hankel spectrum"; 
      max_freq = fs->size * ss->sonogram_cutoff;
      break;
    case chebyshev: 
      xlabel = "Chebyshev spectrum"; 
      max_freq = fs->size * ss->sonogram_cutoff;
      break;
    case legendre: 
      xlabel = "Legendre spectrum"; 
      max_freq = fs->size * ss->sonogram_cutoff;
      break;
    }
  if (ss->fft_style == normal_fft)
    {
      fap = fp->axis;
      if (ss->dBing) 
	{
	  max_val = 0.0; 
	  min_val = -60.0;
	}
      else 
	{
	  if ((!(ss->normalize_fft)) && (fap)) 
	    {
	      max_val = fap->ymax; 
	      min_val = fap->ymin;
	    }
	  else
	    {
	      max_val = 1.0;
	      if (ss->transform_type == fourier) min_val = 0.0; else min_val = -1.0;
	    }
	}
      fp->axis = make_axis_info(cp,
				0.0,max_freq,
				min_val,max_val,
				xlabel,
				0.0,max_freq,
				min_val,max_val,
				fap);
    }
  else make_sonogram_axes(cp);
  fs->data = fp->data;
  return(1);
}

void make_sonogram_axes(chan_info *cp)
{
  snd_state *ss;
  fft_info *fp;
  axis_info *ap;
  float max_freq,yang;
  char *xlabel;
  fp = cp->fft;
  if (fp)
    {
      ap = cp->axis;
      ss = cp->state;
      if (ss->transform_type == fourier)
	{
	  if ((ss->logxing) || (ss->fft_style == spectrogram))
	    max_freq = ss->sonogram_cutoff;
	  else max_freq = ss->sonogram_cutoff * (float)snd_SRATE(cp) * 0.5;
	}
      else max_freq = fp->current_size * ss->sonogram_cutoff;
      yang = fmod(ss->spectro_yangle,360.0);
      if (yang < 0.0) yang += 360.0;
      if (ss->fft_style == spectrogram)
	{
	  if (yang < 45.0) xlabel = snd_string_frequency;
	  else if (yang < 135.0) xlabel = snd_string_time;
	  else if (yang < 225.0) xlabel = snd_string_frequency_reversed;
	  else if (yang < 315.0) xlabel = snd_string_time_reversed;
	  else xlabel = snd_string_frequency;
	}
      else xlabel = snd_string_time;
      fp->axis = make_axis_info(cp,ap->x0,ap->x1,0.0,max_freq,xlabel,ap->x0,ap->x1,0.0,max_freq,fp->axis);
    }
}

static int apply_fft_window(fft_state *fs)
{
  /* apply the window, reading data if necessary, resetting IO blocks to former state */
  int i,ind0,toploc;
  float *window,*fft_data;
  int val;   /* this was a float??? 20-Jan-98 */
  snd_fd *sf;
  chan_info *cp;
  snd_state *ss;
  ss = fs->ss;
  cp = (chan_info *)(fs->chan);
  fft_data = fs->data;
  ind0 = graph_low_SAMPLE(cp) + fs->beg;
  sf = init_sample_read(ind0,cp,READ_FORWARD);
  switch (ss->transform_type)
    {
    case fourier:
      window = (float *)((fft_window_state *)(fs->wp))->window;
      toploc = fs->size / (float)(1 + fs->zero_pad);
      for (i=0;i<toploc;i++) 
	{
	  NEXT_SAMPLE(val,sf);
	  fft_data[i] = window[i] * val;
	}
      if (toploc < fs->size)
	{
	  for (i=toploc;i<fs->size;i++) fft_data[i] = 0.0;
	}
      decrement_fft_window_use(fs);
      break;
    case hankel:
      toploc = fs->size; 
      for (i=0;i<toploc;i++) fs->hwin[i] = next_sample(sf);
      make_abel_transformer(toploc);
      abel(fs->hwin,fft_data);
      break;
    case wavelet:
      for (i=0;i<fs->size;i++) 
	{
	  fft_data[i] = next_sample(sf);
	}
      wavelet_transform(fft_data,fs->size,0,wavelet_data[ss->wavelet_choice],wavelet_sizes[ss->wavelet_choice]);
      /* wavelet_transform(fft_data,fs->size,-1,wavelet_data[ss->wavelet_choice],wavelet_sizes[ss->wavelet_choice]); */
      /* tested by using inverse and viewing result in fft window compared to original wave */
      free_snd_fd(sf);
      return(5);  /* wavelet is done by this point (not a work-proc) */
      break;
    case chebyshev:
      for (i=0;i<fs->size;i++) 
	{
	  fft_data[i] = next_sample(sf);
	}
      chebyshev_transform(fft_data,fs->size);
      free_snd_fd(sf);
      return(5);  /* done by this point (not a work-proc) */
      break;
    case legendre:
      for (i=0;i<fs->size;i++) 
	{
	  fft_data[i] = next_sample(sf);
	}
      legendre_transform(fft_data,fs->size);
      free_snd_fd(sf);
      return(5);  /* done by this point (not a work-proc) */
      break;
    }
  free_snd_fd(sf);
  return(1);
}

static int display_snd_fft(fft_state *fs)
{
  chan_info *cp;
  snd_state *ss;
  snd_info *sp;
  cp = (chan_info *)(fs->chan);
  ss = cp->state;
  sp = cp->sound;
  if (ss->fft_style == normal_fft)
    {
      set_chan_fft_in_progress(cp,0);
      if (ss->global_fft_size >= 65536) report_in_minibuffer(sp,"");
      display_channel_data(cp,sp,cp->state);
    }
  return(-1);
}

void *make_fft_state(chan_info *cp, int simple)
{
  /* in simple fft case, (snd-xchn.c) calls this and passes it to fft_in_slices */
  /* we can cause the current fft to be re-used by setting slice to 8 */
  fft_state *fs = NULL;
  snd_state *ss;
  axis_info *ap;
  int reuse_old = 0;
  ss = cp->state;
  ap = cp->axis;
  if ((simple) && (cp->fft_data))
    {
      fs = (fft_state *)(cp->fft_data);
      if ((fs->losamp == ap->losamp) && 
	  (fs->size == ss->global_fft_size) &&
	  (fs->wintype == ss->global_fft_window) &&
	  (fs->dBing == ss->dBing) &&
	  (fs->zero_pad == ss->zero_pad) &&
	  (fs->cutoff == ss->sonogram_cutoff) &&
	  (fs->t_type == ss->transform_type) &&
	  (fs->wavelet_choice == ss->wavelet_choice) &&
	  (fs->edit_ctr == cp->edit_ctr))
	reuse_old = 1;
    }
  if (reuse_old)
    fs->slice = 8;
  else
    {
      if (cp->fft_data) cp->fft_data = free_fft_state((fft_state *)(cp->fft_data));
      fs = (fft_state *)calloc(1,sizeof(fft_state));
      fs->s_type = make_snd_pointer_type(FFT_STATE);
      fs->slice = 0;
      fs->chan = cp;
      fs->cutoff = ss->sonogram_cutoff;
      fs->size = ss->global_fft_size;
      fs->zero_pad = ss->zero_pad;
      fs->wintype = ss->global_fft_window;
      fs->dBing = ss->dBing;
      fs->wp = NULL;
      fs->losamp = ap->losamp;
      fs->edit_ctr = cp->edit_ctr;
      fs->ss = ss;
      fs->hwin_size = 0;
      fs->hwin = NULL;
      fs->wavelet_choice = ss->wavelet_choice;
      fs->t_type = ss->transform_type;
    }
  fs->nn = fs->size/2;
  fs->beg = 0;
  if (simple) cp->fft_data = fs; else cp->fft_data = NULL;
  return((void *)fs);
}

Boolean fft_in_slices(void *fftData)
{
  /* return true when done */
  /* slices are: 
   *    create arrays if needed
   *    window/load data
   *    shuffle 
   *    step n times through the fft (100 to 200 per iteration) -- wait for +1 here
   *    return true 
   *
   * since we can be running multiple FFTs at once, not to mention other work procedures,
   * all FFT state needs to be in clientData.
   * 
   * Each slice function returns 0 => call me again, 1 => go to next slice, -1 => quit work proc altogether
   */
  fft_state *fs;
  int res = 0;
  fs = (fft_state *)fftData;
  switch (fs->slice)
    {
    case 0: res = snd_fft_set_up(fs);            break;
    case 1: res = set_up_fft_window(fs);         break;
    case 2: res = make_fft_window(fs);           break;
    case 3: res = apply_fft_window(fs);          break;
    case 4: res = shuffle_fft_state(fs);         break;
    case 5: res = snd_fft(fs);                   break;
    case 6: res = snd_fft_cleanup(fs);           break;
    case 7: res = snd_fft_to_spectrum(fs);       break;
    case 8: res = display_snd_fft(fs);           break;
    default: snd_error("impossible fft slice!"); break;
    }
  if (res == -1) return(TRUE);
  fs->slice += res;
  return(FALSE);
}

/* -------------------------------- SONOGRAM -------------------------------- */
/*
 * calls calculate_fft for each slice, each bin being a gray-scaled rectangle in the display
 */

/* as we run the ffts, we need to save the fft data for subsequent redisplay/printing etc */
/* many of these can be running in parallel so the pointers can't be global */
/* display_channel_data above needs to be smart about updates here -- just new data */

/* this work proc calls a loop by pixels (hop depends on pixel/samps decision)
   each pixel(group) sets up the fft_state pointer with beg reflecting hop
   then loops, each time called, calling fft_in_slices until it returns true.
   then grab that data, update the channel display, look to see if we're
   behind the times, if so cleanup and exit, else jump back to outer loop.
   */

typedef struct {
  int slice;
  int outlim,outer;
  fft_state *fs;
  chan_info *cp;
  int spectrum_size;
  sono_info *scp;
  int beg,hop,losamp,hisamp;
  int done;
  int window;
  int msg_ctr;
  int edit_ctr;
  float old_scale;
  int old_style,old_logxing,t_type,w_choice;
} sonogram_state;

void *make_sonogram_state(chan_info *cp)
{
  sonogram_state *sg;
  fft_state *fs;
  snd_state *ss;
  ss = cp->state;
  sg = (sonogram_state *)calloc(1,sizeof(sonogram_state));
  sg->cp = cp;
  sg->done = 0;
  fs = (fft_state *)make_fft_state(cp,0); /* 0=>not a simple one-shot fft */
  sg->fs = fs;
  sg->msg_ctr = 8;
  sg->t_type = ss->transform_type;
  sg->w_choice = ss->wavelet_choice;
  return((void *)sg);
}

void free_sono_info (chan_info *cp)
{
  int i;
  sono_info *si;
  si = (sono_info *)(cp->sonogram_data);
  if (si)
    {
      if (si->data)
	{
	  for (i=0;i<si->total_slices;i++)
	    {
	      if (si->data[i]) free(si->data[i]);
	    }
	  free(si->data);
	}
      free(si);
      cp->sonogram_data = NULL;
    }
}

static int set_up_sonogram(sonogram_state *sg)
{
  /* return 1 to go on, 2 to quit early */
  sono_info *si;
  axis_info *ap,*fap;
  fft_info *fp;
  chan_info *cp;
  snd_state *ss;
  sonogram_state *lsg = NULL;
  int i,old_state,tempsize;
  cp = sg->cp;
  ss = cp->state;
  if (!(cp->fft)) 
    {
      old_state = ss->fft_style;
      ss->fft_style = normal_fft;
      snd_fft_set_up(sg->fs); 
      display_channel_data(cp,cp->sound,cp->state);
      ss->fft_style = old_state; /* sigh... we just want to know the x limits of the damned thing */
    }
  ap = cp->axis;
  sg->slice = 0;
  sg->outer = 0;
  sg->beg = ap->losamp;
  sg->losamp = ap->losamp;
  sg->hisamp = ap->hisamp;
  sg->window = ss->global_fft_window;
  fp = cp->fft;
  fap = fp->axis;
  sg->outlim = (fap->x_axis_x1-fap->x_axis_x0);
  if (ss->fft_style == spectrogram) sg->outlim /= ss->spectro_hop;
  sg->hop = ceil((float)(ap->hisamp - ap->losamp+1)/(float)(sg->outlim));
  /* if fewer samps than pixels, draw rectangles */
  if (ss->transform_type == fourier)
    sg->spectrum_size = (ss->global_fft_size)/2;
  else sg->spectrum_size = ss->global_fft_size;
  sg->edit_ctr = cp->edit_ctr;
  si = (sono_info *)(cp->sonogram_data);
  if (!si)
    {
      si = (sono_info *)calloc(1,sizeof(sono_info));
      cp->sonogram_data = si;
      si->total_bins = sg->spectrum_size;
      si->total_slices = pow(2.0,ceil(log(sg->outlim)/log(2.0)));
      si->begs = (int *)calloc(si->total_slices,sizeof(int));
      si->data = (float **)calloc(si->total_slices,sizeof(float *));
      for (i=0;i<si->total_slices;i++) si->data[i]=(float *)calloc(si->total_bins,sizeof(float));
    }
  else
    if ((si->total_slices < sg->outlim) || (si->total_bins < sg->spectrum_size))
      {
	for (i=0;i<si->total_slices;i++) if (si->data[i]) {free(si->data[i]); si->data[i] = NULL;}
	tempsize = pow(2.0,ceil(log(sg->outlim)/log(2.0)));
	if (si->total_slices < tempsize) 
	  {
	    free(si->data);
	    si->total_slices = tempsize;
	    si->begs = (int *)realloc(si->begs,si->total_slices*sizeof(int));
	    si->data = (float **)calloc(si->total_slices,sizeof(float *));
	  }
	if (si->total_bins < sg->spectrum_size) si->total_bins = sg->spectrum_size;
	for (i=0;i<si->total_slices;i++) si->data[i]=(float *)calloc(si->total_bins,sizeof(float));
      }
  sg->scp = si;
  si->target_bins = sg->spectrum_size;
  si->active_slices = 0;
  si->target_slices = sg->outlim;
  si->scale = 0.0;
  allocate_sono_rects(ss,si->total_bins);              /* was total_slices by mistake (21-Mar-97) */
  if (cp->last_sonogram)                               /* there was a previous run */
    {
      lsg = (sonogram_state *)(cp->last_sonogram);
      if ((lsg->done) &&                               /* it completed all ffts */
	  (lsg->outlim == sg->outlim) &&               /* the number of ffts is the same */
	  (lsg->spectrum_size == sg->spectrum_size) && /* ditto fft sizes */
	  (lsg->losamp == sg->losamp) &&               /* begins are same */
	  (lsg->hisamp == sg->hisamp) &&               /* ends are same */
	  (lsg->window == sg->window) &&               /* data windows are same */
	  (lsg->t_type == sg->t_type) &&               /* transform types are the same */
	  (lsg->w_choice == sg->w_choice) &&           /* wavelets are the same */
	  (lsg->edit_ctr == sg->edit_ctr))             /* underlying data is the same */
	{
	  sg->outer = sg->outlim;                      /*fake up the run */
	  si->active_slices = si->target_slices;
	  sg->old_scale = lsg->old_scale;
	  si->scale = sg->old_scale;
	  if ((lsg->old_style != ss->fft_style) ||
	      (lsg->old_logxing != ss->logxing))
	    make_sonogram_axes(cp);                    /* may need to fixup frequency axis labels */
	  sg->old_style = ss->fft_style;
	  sg->old_logxing = ss->logxing;
	  return(2);                                   /* so skip the ffts! */
	}
    }
  return(1);
}

static char sono_str[32];

static int run_all_ffts(sonogram_state *sg)
{
  int res;
  fft_state *fs;
  sono_info *si;
  chan_info *cp;
  axis_info *ap;
  snd_state *ss;
  float val;
  int i;
  /* return 0 until done with all ffts, then 1 -- 1 causes cleanup whether done or not */
  /* check for losamp/hisamp change? */
  res = fft_in_slices(sg->fs);
  if (res)
    {
      /* slice is done -- store it and prepare to start the next slice */
      fs = sg->fs;
      cp = sg->cp;
      ss = cp->state;
      si = (sono_info *)(cp->sonogram_data);
      si->begs[si->active_slices] = sg->beg + fs->beg;
      sg->msg_ctr--;
      if (sg->msg_ctr == 0)
	{
	  sprintf(sono_str,"%s: %d%%",
		  (ss->fft_style == sonogram) ? "sonogram" : "spectrogram",
		  (int)(100.0*((float)(si->active_slices)/(float)(si->target_slices))));
	  report_in_minibuffer(cp->sound,sono_str);
	  sg->msg_ctr = 8;
	}
      if (ss->transform_type == fourier)
	{
	  for (i=0;i<sg->spectrum_size;i++) 
	    {
	      val = fs->data[i];
	      if (val > si->scale) si->scale = val;
	      si->data[si->active_slices][i] = val;
	    }
	}
      else
	{
	  for (i=0;i<sg->spectrum_size;i++) 
	    {
	      val = fs->data[i];
	      if (val < 0.0) val = -val;  /* kinda dubious but I can't think of a good alternative */
	      if (val > si->scale) si->scale = val;
	      si->data[si->active_slices][i] = val;
	    }
	}
      si->active_slices++;
      sg->outer++;
      if ((sg->outer == sg->outlim) || (!(cp->ffting)) || (ss->fft_style == normal_fft)) return(1);
      fs->beg += sg->hop;
      fs->slice = 0;
      ap = cp->axis;
      if ((sg->losamp != ap->losamp) || (sg->hisamp != ap->hisamp)) 
	{
	  fs->beg = 0;
	  return(-1);
	}
    }
  return(0);
}

static int cleanup_sonogram(sonogram_state *sg)
{
  chan_info *cp;
  /* data has already been placed on cp->sonogram_state, so we need only clear the fft_state struct */
  if (sg)
    {
      cp = sg->cp;
      if (sg->fs) sg->fs = free_fft_state(sg->fs);
      cp->fft_data = NULL;
      set_chan_fft_in_progress(cp,0);
      display_channel_data(cp,cp->sound,cp->state);
      if (cp->last_sonogram) free(cp->last_sonogram);
      if (sg->outer == sg->outlim) sg->done = 1;
      sg->old_scale = (sg->scp)->scale;
      cp->last_sonogram = sg;
      sono_str[0]='\0';
      report_in_minibuffer(cp->sound,sono_str);
    }
  return(1);
}

Boolean sonogram_in_slices(void *sono)
{
  sonogram_state *sg = (sonogram_state *)sono;
  int res = 0;
  switch (sg->slice)
    {
    case 0: res = set_up_sonogram(sg); break; /* return 1 to go on, 2 to quit early */
    case 1: res = run_all_ffts(sg);    break; /* return 0 until done with all ffts, then 1 -- 1 causes cleanup whether done or not */
    case 2: res = cleanup_sonogram(sg); return(TRUE); break;
    }
  sg->slice += res;
  return(FALSE);
}