: $Id: tsa.mod,v 1.18 2011/10/25 21:45:39 billl Exp $ 


NEURON {
 SUFFIX tsa
 GLOBAL INSTALLED
 GLOBAL verbose
}

PARAMETER {
  INSTALLED=0
  verbose=0
}

VERBATIM
#include "misc.h"
static double *x1x, *y1y, *z1z;
//void spctrm(double data[], double p[], int m, int k);
extern double dfftpow(double* x,int n,double* ppow,int powlen,int* fftlen);

#define WINDOW(j,a,b) (1.0-fabs((((j)-1)-(a))*(b)))       /* Bartlett */

double *nrvector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
	double *v;

	v=(double *)malloc((size_t) ((nh-nl+1+1)*sizeof(double)));
	if (!v) { nrerror("allocation failure in vector()"); return 0x0;}
	return v-nl+1;
}

void nrfree_vector(double *v, long nl, long nh)
/* free a double vector allocated with vector() */
{
  free((char*) (v+nl-1));
}

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

static double nrsqrarg;
#define SQR(a) ((nrsqrarg=(a)) == 0.0 ? 0.0 : nrsqrarg*nrsqrarg)


void nrfour1(double mdata[], unsigned long nn, int isign)
{
	unsigned long n,mmax,m,j,istep,i;
	double wtemp,wr,wpr,wpi,wi,theta;
	double tempr,tempi;

	n=nn << 1;
	j=1;
	for (i=1;i<n;i+=2) {
		if (j > i) {
			SWAP(mdata[j],mdata[i]);
			SWAP(mdata[j+1],mdata[i+1]);
		}
		m=n >> 1;
		while (m >= 2 && j > m) {
			j -= m;
			m >>= 1;
		}
		j += m;
	}
	mmax=2;
	while (n > mmax) {
		istep=mmax << 1;
		theta=isign*(6.28318530717959/mmax);
		wtemp=sin(0.5*theta);
		wpr = -2.0*wtemp*wtemp;
		wpi=sin(theta);
		wr=1.0;
		wi=0.0;
		for (m=1;m<mmax;m+=2) {
			for (i=m;i<=n;i+=istep) {
				j=i+mmax;
				tempr=wr*mdata[j]-wi*mdata[j+1];
				tempi=wr*mdata[j+1]+wi*mdata[j];
				mdata[j]=mdata[i]-tempr;
				mdata[j+1]=mdata[i+1]-tempi;
				mdata[i] += tempr;
				mdata[i+1] += tempi;
			}
			wr=(wtemp=wr)*wpr-wi*wpi+wr;
			wi=wi*wpr+wtemp*wpi+wi;
		}
		mmax=istep;
	}
}

void nrspctrm(double mdata[], double p[], int m, int k)
{
  // assume overlap = 1
	void nrfour1(double mdata[], unsigned long nn, int isign);
	int mm,m44,m43,m4,kk,joffn,joff,j2,j,c=0;
	double w,facp,facm,*w1,*w2,sumw=0.0,den=0.0;

	mm=m+m;
	m43=(m4=mm+mm)+3;
	m44=m43+1;
	w1=nrvector(1,m4);
	w2=nrvector(1,m);
	facm=m;
	facp=1.0/m;
	for (j=1;j<=mm;j++) sumw += SQR(WINDOW(j,facm,facp));
	for (j=1;j<=m;j++) p[j]=0.0;
	for (j=1;j<=m;j++) w2[j] = mdata[c++];
	for (kk=1;kk<=k;kk++) {
		for (joff = -1;joff<=0;joff++) {
			for (j=1;j<=m;j++) w1[joff+j+j]=w2[j];
			for (j=1;j<=m;j++) w2[j] = mdata[c++];
			joffn=joff+mm;
			for (j=1;j<=m;j++) w1[joffn+j+j]=w2[j];
		}
		for (j=1;j<=mm;j++) {
			j2=j+j;
			w=WINDOW(j,facm,facp);
			w1[j2] *= w;
			w1[j2-1] *= w;
		}
		nrfour1(w1,mm,1);
		p[1] += (SQR(w1[1])+SQR(w1[2]));
		for (j=2;j<=m;j++) {
			j2=j+j;
			p[j] += (SQR(w1[j2])+SQR(w1[j2-1])
				+SQR(w1[m44-j2])+SQR(w1[m43-j2]));
		}
		den += sumw;
	}
	den *= m4;
	for (j=1;j<=m;j++) p[j] /= den;
	nrfree_vector(w2,1,m);
	nrfree_vector(w1,1,m4);
}

double mymax(double x, double y){
  return x > y ? x : y;
}

/* ********************************************************************* */

void mysvd(int m, int n, double** u, double w[], double** v, int* ierr)
/*
 *   This subroutine is a translation of the Algol procedure svd,
 *   Num. Math. 14, 403-420(1970) by Golub and Reinsch.
 *   Handbook for Auto. Comp., Vol II-Linear Algebra, 134-151(1971).
 *
 *   This subroutine determines the singular value decomposition
 *        t
 *   A=usv  of a real m by n rectangular matrix, where m is greater
 *   then or equal to n.  Householder bidiagonalization and a variant
 *   of the QR algorithm are used.
 *  
 *
 *   On input.
 *
 *      m is the number of rows of A (and u).
 *
 *      n is the number of columns of A (and u) and the order of v.
 *
 *      u contains the rectangular input matrix A to be decomposed.
 *
 *   On output.
 *
 *      w contains the n (non-negative) singular values of a (the
 *        diagonal elements of s).  they are unordered.  if an
 *        error exit is made, the singular values should be correct
 *        for indices ierr+1,ierr+2,...,n.
 *
 *      u contains the matrix u (orthogonal column vectors) of the
 *        decomposition.
 *        if an error exit is made, the columns of u corresponding
 *        to indices of correct singular values should be correct.
 *
 *      v contains the matrix v (orthogonal) of the decomposition.
 *        if an error exit is made, the columns of v corresponding
 *        to indices of correct singular values should be correct.
 *
 *      ierr is set to
 *        zero       for normal return,
 *        k          if the k-th singular value has not been
 *                   determined after 30 iterations,
 *        -1         if memory allocation fails.
 *
 *   Questions and comments should be directed to B. S. Garbow,
 *   Applied Mathematics division, Argonne National Laboratory
 *
 *   Modified to eliminate machep
 *
 *   Translated to C by Michiel de Hoon, Human Genome Center,
 *   University of Tokyo, for inclusion in the C Clustering Library.
 *   This routine is less general than the original svd routine, as
 *   it focuses on the singular value decomposition as needed for
 *   clustering. In particular,
 *     - We require m >= n
 *     - We calculate both u and v in all cases
 *     - We pass the input array A via u; this array is subsequently
 *       overwritten.
 *     - We allocate for the array rv1, used as a working space,
 *       internally in this routine, instead of passing it as an
 *       argument. If the allocation fails, svd sets *ierr to -1
 *       and returns.
 *   2003.06.05
 */
{ int i, j, k, i1, k1, l1, its;
  double c,f,h,s,x,y,z;
  int l = 0;
  double g = 0.0;
  double scale = 0.0;
  double anorm = 0.0;
  double* rv1 = (double*)malloc(n*sizeof(double));
  if (!rv1)
  { *ierr = -1;
    return;
  }
  *ierr = 0;
  /* Householder reduction to bidiagonal form */
  for (i = 0; i < n; i++)
  { l = i + 1;
    rv1[i] = scale * g;
    g = 0.0;
    s = 0.0;
    scale = 0.0;
    for (k = i; k < m; k++) scale += fabs(u[k][i]);
    if (scale != 0.0)
    { for (k = i; k < m; k++)
      { u[k][i] /= scale;
        s += u[k][i]*u[k][i];
      }
      f = u[i][i];
      g = (f >= 0) ? -sqrt(s) : sqrt(s);
      h = f * g - s;
      u[i][i] = f - g;
      if (i < n-1)
      { for (j = l; j < n; j++)
        { s = 0.0;
          for (k = i; k < m; k++) s += u[k][i] * u[k][j];
          f = s / h;
          for (k = i; k < m; k++) u[k][j] += f * u[k][i];
        }
      }
      for (k = i; k < m; k++) u[k][i] *= scale;
    }
    w[i] = scale * g;
    g = 0.0;
    s = 0.0;
    scale = 0.0;
    if (i<n-1)
    { for (k = l; k < n; k++) scale += fabs(u[i][k]);
      if (scale != 0.0)
      { for (k = l; k < n; k++)
        { u[i][k] /= scale;
          s += u[i][k] * u[i][k];
        }
        f = u[i][l];
        g = (f >= 0) ? -sqrt(s) : sqrt(s);
        h = f * g - s;
        u[i][l] = f - g;
        for (k = l; k < n; k++) rv1[k] = u[i][k] / h;
        for (j = l; j < m; j++)
        { s = 0.0;
          for (k = l; k < n; k++) s += u[j][k] * u[i][k];
          for (k = l; k < n; k++) u[j][k] += s * rv1[k];
        }
        for (k = l; k < n; k++)  u[i][k] *= scale;
      }
    }
    anorm = mymax(anorm,fabs(w[i])+fabs(rv1[i]));
  }
  /* accumulation of right-hand transformations */
  for (i = n-1; i>=0; i--)
  { if (i < n-1)
    { if (g != 0.0)
      { for (j = l; j < n; j++) v[j][i] = (u[i][j] / u[i][l]) / g;
        /* double division avoids possible underflow */
        for (j = l; j < n; j++)
        { s = 0.0;
          for (k = l; k < n; k++) s += u[i][k] * v[k][j];
          for (k = l; k < n; k++) v[k][j] += s * v[k][i];
        }
      }
    }
    for (j = l; j < n; j++)
    { v[i][j] = 0.0;
      v[j][i] = 0.0;
    }
    v[i][i] = 1.0;
    g = rv1[i];
    l = i;
  }
  /* accumulation of left-hand transformations */
  for (i = n-1; i >= 0; i--)
  { l = i + 1;
    g = w[i];
    if (i!=n-1)
      for (j = l; j < n; j++) u[i][j] = 0.0;
    if (g!=0.0)
    { if (i!=n-1)
      { for (j = l; j < n; j++)
        { s = 0.0;
          for (k = l; k < m; k++) s += u[k][i] * u[k][j];
          /* double division avoids possible underflow */
          f = (s / u[i][i]) / g;
          for (k = i; k < m; k++) u[k][j] += f * u[k][i];
        }
      }
      for (j = i; j < m; j++) u[j][i] /= g;
    }
    else
      for (j = i; j < m; j++) u[j][i] = 0.0;
    u[i][i] += 1.0;
  }
  /* diagonalization of the bidiagonal form */
  for (k = n-1; k >= 0; k--)
  { k1 = k-1;
    its = 0;
    while(1)
    /* test for splitting */
    { for (l = k; l >= 0; l--)
      { l1 = l-1;
        if (fabs(rv1[l]) + anorm == anorm) break;
        /* rv1[0] is always zero, so there is no exit
         * through the bottom of the loop */
        if (fabs(w[l1]) + anorm == anorm)
        /* cancellation of rv1[l] if l greater than 0 */
        { c = 0.0;
          s = 1.0;
          for (i = l; i <= k; i++)
          { f = s * rv1[i];
            rv1[i] *= c;
            if (fabs(f) + anorm == anorm) break;
            g = w[i];
            h = sqrt(f*f+g*g);
            w[i] = h;
            c = g / h;
            s = -f / h;
            for (j = 0; j < m; j++)
            { y = u[j][l1];
              z = u[j][i];
              u[j][l1] = y * c + z * s;
              u[j][i] = -y * s + z * c;
            }
          }
          break;
        }
      }
      /* test for convergence */
      z = w[k];
      if (l==k) /* convergence */
      { if (z < 0.0)
        /*  w[k] is made non-negative */
        { w[k] = -z;
          for (j = 0; j < n; j++) v[j][k] = -v[j][k];
        }
        break;
      }
      else if (its==30)
      { *ierr = k;
        break;
      }
      else
      /* shift from bottom 2 by 2 minor */
      { its++;
        x = w[l];
        y = w[k1];
        g = rv1[k1];
        h = rv1[k];
        f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
        g = sqrt(f*f+1.0);
        f = ((x - z) * (x + z) + h * (y / (f + (f >= 0 ? g : -g)) - h)) / x;
        /* next qr transformation */
        c = 1.0;
        s = 1.0;
        for (i1 = l; i1 <= k1; i1++)
        { i = i1 + 1;
          g = rv1[i];
          y = w[i];
          h = s * g;
          g = c * g;
          z = sqrt(f*f+h*h);
          rv1[i1] = z;
          c = f / z;
          s = h / z;
          f = x * c + g * s;
          g = -x * s + g * c;
          h = y * s;
          y = y * c;
          for (j = 0; j < n; j++)
          { x = v[j][i1];
            z = v[j][i];
            v[j][i1] = x * c + z * s;
            v[j][i] = -x * s + z * c;
          }
          z = sqrt(f*f+h*h);
          w[i1] = z;
          /* rotation can be arbitrary if z is zero */
          if (z!=0.0)
          { c = f / z;
            s = h / z;
          }
          f = c * g + s * y;
          x = -s * g + c * y;
          for (j = 0; j < m; j++)
          { y = u[j][i1];
            z = u[j][i];
            u[j][i1] = y * c + z * s;
            u[j][i] = -y * s + z * c;
          }
        }
        rv1[l] = 0.0;
        rv1[k] = f;
        w[k] = x;
      }
    }
  }
  free(rv1);
  return;
}

int mycompare( const void* a, const void* b ) {
  double* arg1 = (double*) a;
  double* arg2 = (double*) b;
  if( *arg1 < *arg2 ) return -1;
  else if( *arg1 == *arg2 ) return 0;
  else return 1;
}  

double* myspectrum(double* v1,int dc,int* anslen){
  // n data pts will be divided into k partitions of size m
  // the spectrum will have m values from 0 to m/2 cycles/dt.
  //  n = (2*k+1)*m

  // data set
  //  Vect* v1 = vector_arg(1);

  //  int dc = v1->capacity();
  int mr;
  //  if (ifarg(2)) mr = (int)(*getarg(2)); else mr = dc/8;
  mr = dc / 8;

  // make sure the length of partitions is integer power of 2
  int m = 1;
  while(m<mr) m*=2;
  *anslen = m;

  int k = (int) (ceil(((double)(dc)/m-1.)/2.));
  int n = (2*k+1)*m;

  int i;

  double *mdata = (double *)calloc(n,(unsigned)sizeof(double));
  for (i=0;i<dc;++i) mdata[i] = v1[i]; //v1->elem(i);

  //if (ans->capacity() < m) ans->resize(m);
  double* ans = (double *)calloc(m,(unsigned)sizeof(double));
  //spctrm(&mdata[0], &ans->elem(0)-1, m, k);
  nrspctrm(&mdata[0], &ans[0], m, k);

  free((char *)mdata);

  return ans;
}

double myspct(void* v){
  double* x;
  int n = vector_instance_px(v,&x) , anslen = 0 , i = 0;
  double* ans = myspectrum(x,n,&anslen);
  for(i=0;i<anslen;i++) printf("%g ",ans[i]);
  printf("\n");
  double* p;
  int outlen = vector_arg_px(1,&p);
  for(i=0;i<outlen && i<anslen;i++) p[i]=ans[i];
  free(ans);
  return 1.0;
}

int iinrange(double* x,int sz,double dmin,double dmax){
  int i = 0, cnt = 0;
  for(i=0;i<sz;i++) if(x[i]>=dmin && x[i]<=dmax) cnt++;
  return cnt;
}

double inrange(void* v){
  double* x , dmin, dmax;
  int sz = vector_instance_px(v,&x) , i , cnt = 0;
  return (double) iinrange(x,sz,*getarg(1),*getarg(2));
}

extern double dfftpow(double* x,int n,double* ppow,int powlen,int* fftlen);

ENDVERBATIM

: get averagecorrelations between channels, but first bandpass filter them
: tscoravgband(output vector,vec-list of eeg time series,window size,increment,vec:60/120Hz power,vec:saturation,
: vec:delta power,vec:theta,vec:beta,vec:gamma,vec:ripple,sampling rate,vec:low-cutoffs,vec:high-cutoffs, [optional
: output list of vectors for filtered time-series]
:  --
: low and high cutoffs are different bands that will be kept, so can do theta+gamma, i.e. lo=4,31 , hi=12,100
: or can just do theta : lo=4 , hi=12
FUNCTION tscoravgband () {
  VERBATIM

  int i;

  double dRet = 0.0;
  ListVec* pTS = 0;

  double* pcvin = 0;
  double** pfilter = 0x0;

  double** pTSN = 0; // normalized time series

  double* vTCor = 0; //correlation matrix abs sum avg for a time point
  int corsz = vector_arg_px(1,&vTCor);

  double* pTMP = 0x0;

  ListVec* pTSOut = 0;

  pTS = AllocListVec(*hoc_objgetarg(2)); //list of eeg time series
  
  if(!pTS){
    printf("tscoravg ERRA: problem initializing 1st 2 args!\n");
    goto CLEANUP;
  }
  int iChans = pTS->isz; // # of channels
  if(iChans < 2) {
    printf("tscoravg ERRB: must have at least 2 EEG channels!\n");
    goto CLEANUP;
  }
  int iWinSz = (int) *getarg(3); //window size for analysis
  int iSamples = pTS->plen[0]; //make sure all time series have same length
  int iINC = ifarg(4)?(int)*getarg(4):1; 
  double *phz = 0, *psat = 0, *pfft = 0; int sz = 0 , iNoiseTest = 0;
  if(ifarg(5) && ifarg(6)){ //do saturation and 60/120 Hz Frequency test
    if((sz=vector_arg_px(5,&phz))!=corsz){
      printf("tscoravg ERRG: invalid size for phz should be %d, is %d!\n",corsz,sz);
      goto CLEANUP;
    } else if((sz=vector_arg_px(6,&psat))!=corsz){
      printf("tscoravg ERRH: invalid size for psat should be %d, is %d!\n",corsz,sz);
      goto CLEANUP;
    } else {	      
      iNoiseTest = 1;
      pfft = (double*)calloc(iWinSz+1,sizeof(double));
    }
  }
  int iOutSz = (iSamples - iWinSz + 1) / iINC;//output List size

  double *pdelta=0,*ptheta=0,*pbeta=0,*pgamma=0,*pripple=0;
  if(ifarg(7) && vector_arg_px(7,&pdelta)<iOutSz){ printf("bad pdelta size\n"); goto CLEANUP; } //delta 0-3
  if(ifarg(8) && vector_arg_px(8,&ptheta)<iOutSz){ printf("bad ptheta size\n"); goto CLEANUP; } //theta 4-12
  if(ifarg(9)&& vector_arg_px(9,&pbeta)<iOutSz){ printf("bad pbeta size\n"); goto CLEANUP; } //beta 13-29
  if(ifarg(10)&& vector_arg_px(10,&pgamma)<iOutSz){ printf("bad pgamma size\n"); goto CLEANUP; } //gamma 30-100
  if(ifarg(11)&& vector_arg_px(11,&pripple)<iOutSz){ printf("bad pripple size\n"); goto CLEANUP; } //ripple 101-300

  double sampr = *getarg(12);

  double *plohz=0x0,*phihz=0x0;
  int iFilters = 1;
  if((iFilters=vector_arg_px(13,&plohz))!=(i=vector_arg_px(14,&phihz))) {
    printf("invalid arg 13,14 lohz/hihz filter sizes: %d %d\n",iFilters,i); goto CLEANUP;
  } else if(iFilters<1) { printf("no filters specified!\n"); goto CLEANUP; }
  int N = 1, M = 1025;
  while(N < iWinSz) N*=2;
  if(N>M) M=N+1;
  pfilter = getdouble2D(iFilters,N);
  //normalize low/high cutoff frequencies by sampling rate (should be between 0 and 0.5)
  for(i=0;i<iFilters;i++){ plohz[i]/=sampr; phihz[i]/=sampr; }
  for(i=0;i<iFilters;i++){ wsfirBP(pfilter[i], M, 1, plohz[i], phihz[i], N); wrap(pfilter[i],N,M); }
  pcvin = (double*) calloc(N,sizeof(double));

  if(verbose) printf("iWinSz=%d iINC=%d iSamples=%d\n",iWinSz,iINC,iSamples);
  for(i=1;i<iChans;i++) {
    if(pTS->plen[i] != iSamples){
      printf("tscoravg ERRC: time series of unequal size %d %d!\n",iSamples,pTS->plen[i]);
      goto CLEANUP;
    }
  }

  if(corsz < iOutSz){
    printf("tscoravg ERRD: need output vector of at least %d as arg 1, have only %d!\n",iOutSz,corsz);
    goto CLEANUP;
  }
  if(verbose) printf("iOutSz=%d\n",iOutSz);

  pTSN = getdouble2D(iChans,2*N); if(verbose) printf("got pTSN\n");
  if(!pTSN){
    printf("tscor ERRD: out of memory!\n");
    goto CLEANUP;
  }

  if(ifarg(15) && (pTSOut = AllocListVec(*hoc_objgetarg(15)))) { //optional - copy filtered output
    if(pTSOut->isz!=iChans) { printf("incorrect output vec-list size: %d %d\n",pTSOut->isz,iChans); goto CLEANUP; }
    for(i=0;i<iChans;i++) if(pTSOut->plen[i] < iSamples) { printf("incorrect output vec-size %d %d\n",pTSOut->plen[i],iSamples);
      goto CLEANUP; }}

  pTMP = (double*) calloc(2*N,sizeof(double));

  double* psums = (double*)calloc(iChans,sizeof(double)); //stores sum of values in channel's time window
  double* psums2 = (double*)calloc(iChans,sizeof(double));//stores sum of squared values in channel's time window

  double dVal = 0.0;
  int iStartIDX = 0 , ierr = 0 , iOutIDX = 0 , iFilt = 0;
  for(iOutIDX=0;iOutIDX<corsz;iOutIDX++) vTCor[iOutIDX] = 0.;
  if( phz ) for(iOutIDX=0;iOutIDX<corsz;iOutIDX++) phz[iOutIDX]=0.;
  if( psat ) for(iOutIDX=0;iOutIDX<corsz;iOutIDX++) psat[iOutIDX]=0.;
  iOutIDX=0;
  double CP = iChans*(iChans-1.0)/2.0;
  for(iStartIDX=0;iOutIDX<iOutSz && iStartIDX+iWinSz<iSamples;iStartIDX+=iINC,iOutIDX++){
    if(iOutIDX%1000==0) printf("%d\n",iOutIDX);
    int IDX = iStartIDX , iChan = 0 , iEndIDX = iStartIDX + iWinSz;
    int iCurSz = iEndIDX - iStartIDX;
    if(phz) phz[iOutIDX]=0.; if(psat)psat[iOutIDX]=0.;
    if(pdelta) pdelta[iOutIDX]=0.; if(ptheta) ptheta[iOutIDX]=0.;
    if(pbeta) pbeta[iOutIDX]=0.; if(pgamma) pgamma[iOutIDX]=0.; if(pripple) pripple[iOutIDX]=0.;

    double dsatavg = 0., dhzavg = 0.;
    for(iChan=0;iChan<iChans;iChan++){       
      i=0;      
      int iSat = 0; 
      double* pc = &pTS->pv[iChan][iStartIDX]; //saturation test
      for(IDX=iStartIDX;IDX<iEndIDX;IDX++,pc++)if((*pc>=900. && *pc<=1010.) || (*pc>=-1010. && *pc<=-900.)) iSat++;
      dsatavg += (double)iSat/ (double)iWinSz;
      int fftlen = 0; memset(pfft,0,sizeof(double)*(iWinSz+1));
      if(dfftpow(&pTS->pv[iChan][iStartIDX],iWinSz,pfft,iWinSz+1,&fftlen)){       //60,120 Hz Frequency test
	int f = 0; pfft[0]=0.0; double fsum = 0.0; double* ppf = pfft;
	for(f=0;f<fftlen && f<=sampr/2;f++) fsum += *ppf++; //normalize power btwn 0-sampr/2Hz to 1.0
	for(f=55;f<=65 && f<fftlen;f++) dhzavg += pfft[f] / fsum; //60Hz contribution
	for(f=115;f<=125 && f<fftlen;f++) dhzavg += pfft[f] / fsum; //120Hz contribution
	if(pdelta) {
	  fsum=0.0;
	  for(f=0;f<fftlen;f++) fsum+=pfft[f];
	  for(f=0;f<=3;f++) pdelta[iOutIDX] += pfft[f]/fsum;
	  if(ptheta) for(f=4;f<=12;f++) ptheta[iOutIDX] += pfft[f]/fsum;
	  if(pbeta) for(f=13;f<=29;f++) pbeta[iOutIDX] += pfft[f]/fsum;
	  if(pgamma) for(f=30;f<=100;f++) pgamma[iOutIDX] += pfft[f]/fsum;
	  if(pripple) for(f=101;f<=300 && f<fftlen;f++) pripple[iOutIDX] += pfft[f]/fsum;
	}
      }	
    }
    phz[iOutIDX] = dhzavg / (double)iChans;
    psat[iOutIDX] = dsatavg / (double)iChans;
    if(pdelta) pdelta[iOutIDX] /= (double)iChans;
    if(ptheta) ptheta[iOutIDX] /= (double)iChans;
    if(pbeta) pbeta[iOutIDX] /= (double)iChans;
    if(pgamma) pgamma[iOutIDX] /= (double)iChans;
    if(pripple) pripple[iOutIDX] /= (double)iChans;

    //bandpass filter and normalize
    for(iChan=0;iChan<iChans;iChan++){
      memset(pTMP,0,sizeof(double)*2*N);
      double* pC = &pTS->pv[iChan][iStartIDX], *pIN = pcvin;
      if(0) for(IDX=iStartIDX;IDX<iEndIDX;IDX++) *pIN++ = *pC++; //copy input
      for(iFilt=0;iFilt<iFilters;iFilt++) {        
        if(1) {
          memset(pTSN[iChan],0,sizeof(double)*N);//zero it out first...
          pIN = pcvin; pC = &pTS->pv[iChan][iStartIDX];
          for(IDX=iStartIDX;IDX<iEndIDX;IDX++) *pIN++ = *pC++; //copy input
        }
        convlv(pcvin-1,N,pfilter[iFilt]-1,M,1,pTSN[iChan]-1);//do filtering
        for(IDX=0;IDX<N;IDX++) pTMP[IDX] += pTSN[iChan][IDX]; 
      }
      for(IDX=0;IDX<N;IDX++) pTSN[iChan][IDX] = pTMP[IDX]; //copy filtered output to pTSN
      if(pTSOut)for(IDX=0;IDX<N && IDX+iStartIDX<pTSOut->plen[iChan];IDX++) pTSOut->pv[iChan][iStartIDX+IDX]=pTSN[iChan][IDX];

      psums[iChan]=psums2[iChan]=0.;
      pC = pTSN[iChan]; //pointer to start of channel
      for(i=0;i<N;i++,pC++){ 
	psums[iChan] += *pC;//for average
	psums2[iChan] += (*pC * *pC); //for std-dev
      }
      double davg = psums[iChan] / N;
      double dstdev = psums2[iChan]/N - davg*davg; 
      if(dstdev > 0. ) dstdev = sqrt(dstdev); else dstdev=1.; //make sure no problems with nan
      if(dstdev <= 0.) dstdev = 1.0;
      pC = pTSN[iChan]; //pointer to start of channel
      for(i=0;i<N;i++,pC++) *pC = (*pC - davg)/dstdev; //normalization
    }
    //get average of pairwise correlations across channels
    int iC1,iC2;
    double dpsum = 0.0;
    for(iC1=0;iC1<iChans;iC1++){
      for(iC2=0;iC2<iC1;iC2++){
        double r = 0.0, *p1 = pTSN[iC1], *p2 = pTSN[iC2];
	for(i=0;i<N;i++) r += (*p1++ * *p2++);
        r /= (double) N;
        dpsum += r;
      }
    }
    dpsum /= CP;
    if(dpsum > 1.0 || dpsum < -1.0) printf("out of bounds = %g!\n",dpsum);
    vTCor[iOutIDX] = dpsum;
  }
  dRet = 1.0;
  printf("\n");

CLEANUP:
  if(pTS) FreeListVec(&pTS);                       if(pTSN) freedouble2D(&pTSN,iChans);
  free(psums);                                     free(psums2); 
  if(pfft) free(pfft);
  if(pfilter) freedouble2D(&pfilter,iFilters);     if(pcvin) free(pcvin);
  if(pTMP) free(pTMP);                             if(pTSOut) FreeListVec(&pTSOut);
  return dRet;
  ENDVERBATIM
}


: normalize values in list of vectors
FUNCTION tsnorm () {
  VERBATIM
  double dRet = 0.0;
  ListVec* pEIG = AllocListVec(*hoc_objgetarg(1));
  if(!pEIG || pEIG->isz < 1 || pEIG->plen[0]<1){
    printf("tsnorm ERRA: problem initializing 1st arg!\n");
    goto TSNCLEANUP;
  }
  double* dmean = 0, *dstdev = 0;
  int cols = vector_arg_px(2,&dmean);
  if(cols!=vector_arg_px(3,&dstdev) || cols!=pEIG->plen[0]){
    printf("tsnorm ERRB: problem initializing 2nd,3rd arg!\n");
    goto TSNCLEANUP;
  }
  int i,j;
  double val;
  for(i=0;i<pEIG->isz;i++)for(j=0;j<cols;j++)pEIG->pv[i][j]=(pEIG->pv[i][j]-dmean[j])/dstdev[j];
  dRet = 1.0;
TSNCLEANUP:
  if(pEIG) FreeListVec(&pEIG);
  return dRet;
  ENDVERBATIM
}

VERBATIM
int ddiff(double* pin,double* pout,int sz){
  int i;
  for(i=0;i<sz-1;i++) pout[i]=pin[i+1]-pin[i];
  return 1;
}
ENDVERBATIM

: get averages from tscorfull into Vectors
: tsgetmeans( 1 cor vec list, 2 phzlist, 3 psatlist, 4 pderivlist, 5 avg cor vec A, 6 avg cor vec B, 7 avg cor vec C, 8 phzt, 9 psatt, 10 pderivt , 11 numchannels , 12 phza factor, 13 vfid, 14 fid, 15 vprctbadch)
: arg 5 will have average correlation without getting rid of any channels on any epochs
: arg 6 will have average correlation with getting rid of saturated and flat channels
: arg 7 will have avg. cor. with getting rid of saturated,flat channels and channels with 60/120Hz power > phzt
: pfctr not used currently
: 13 vector of file ids so can selectively clean up particular files
: 14 file id to clean up -- iff == -1, do all
: 15 % of 'bad' channels for each time period -- channels excluded for "avgC" / # of channels
FUNCTION tsgetmeans () {
VERBATIM
  int iC1,iC2,iChans,idx,jdx;
  double dRet = 0.0;
  ListVec *plvcor = 0, *plvhz = 0, *plvsat = 0, *plvderiv = 0;
  double* pavgA=0,*pavgB=0,*pavgC=0,*pprctbadch=0,*pfid=0,*pbadc=0;
  int vsz = 0, fid=-1; 
  ////////////////////////////////////////////////////////////////////////////////////////////////
  //
  //                                     set up input args.
  //
  if(!(plvcor = AllocListVec(*hoc_objgetarg(1)))){
    printf("tsgetmeans ERRA: couldn't get cor vec list!\n");
    goto MCLEANUP;
  }
  if(!(plvhz = AllocListVec(*hoc_objgetarg(2)))){
    printf("tsgetmeans ERRB: couldn't get phz vec list!\n");
    goto MCLEANUP;
  }
  if(plvhz->isz < plvcor->isz){
    printf("tsgetmeans ERRC: phz vec list size < cor vec list size : %d %d\n",plvhz->isz,plvcor->isz);
    goto MCLEANUP;
  }
  if(!(plvsat = AllocListVec(*hoc_objgetarg(3)))){
    printf("tsgetmeans ERRD: couldn't get psat vec list!\n");
    goto MCLEANUP;
  }
  if(plvsat->isz < plvcor->isz){
    printf("tsgetmeans ERRE: psat vec list size < cor vec list size : %d %d\n",plvsat->isz,plvcor->isz);
    goto MCLEANUP;
  }  
  if(!(plvderiv = AllocListVec(*hoc_objgetarg(4)))){
    printf("tsgetmeans ERRF: couldn't get pderiv vec list!\n");
    goto MCLEANUP;
  }
  if(plvderiv->isz < plvcor->isz){
    printf("tsgetmeans ERRG: pderiv vec list size < cor vec list size : %d %d\n",plvderiv->isz,plvcor->isz);
    goto MCLEANUP;
  }  
  if(vector_arg_px(5,&pavgA) < plvcor->isz || 
     vector_arg_px(6,&pavgB) < plvcor->isz || 
     vector_arg_px(7,&pavgC) < plvcor->isz ||
     vector_arg_px(13,&pfid) < plvcor->isz ||
     vector_arg_px(15,&pprctbadch)< plvcor->isz){
    printf("tsgetmeans ERRH: output vecs must have size >= %d!\n",plvcor->isz);
    goto MCLEANUP;
  }
  double phzt = *getarg(8), psatt = *getarg(9) , pderivt = *getarg(10);
  iChans = (int)*getarg(11);
  double pfctr = *getarg(12);
  fid = (int)*getarg(14);
  pbadc = (double*) malloc(iChans*sizeof(double));
  //
  //
  ////////////////////////////////////////////////////////////////////////////////////////////////

  /////////
  // do the calculations  
  double dSumA,dSumB,dSumC;
  int cntA,cntB,cntC,cntBadCh;
  for(idx=0;idx<plvcor->isz;idx++){
    if(fid>=0 && pfid[idx]!=fid) continue;
    cntA=cntB=cntC=jdx=0;
    dSumA=dSumB=dSumC=0.;
    for(iC1=0;iC1<iChans;iC1++){ //mark bad channels
      if(plvsat->pv[idx][iC1] >= psatt ||
         plvderiv->pv[idx][iC1] >= pderivt)
        pbadc[iC1]=1;
      else
        pbadc[iC1]=0;
    }
    for(iC1=0;iC1<iChans;iC1++){
      if(pbadc[iC1])continue;
      for(iC2=0;iC2<iC1;iC2++,jdx++){
        dSumA += plvcor->pv[idx][jdx];
        cntA++;
        if(pbadc[iC2]) continue;
        dSumB += plvcor->pv[idx][jdx];
        cntB++;
      }
    }
    pavgA[idx] = dSumA / (double) cntA;
    if(cntB) pavgB[idx] = dSumB / (double) cntB; else pavgB[idx] = -666.;
    jdx=0;
    for(iC1=0;iC1<iChans;iC1++) if(plvhz->pv[idx][iC1] >= phzt) pbadc[iC1]=1; //mark bad channels
    for(iC1=0;iC1<iChans;iC1++){
      if(pbadc[iC1]) continue;
      for(iC2=0;iC2<iC1;iC2++,jdx++){
        if(pbadc[iC2]) continue;
        dSumC += plvcor->pv[idx][jdx];
        cntC++;
      }
    }
    if(cntC) pavgC[idx] = dSumC / (double) cntC; else pavgC[idx] = -666.;
    cntBadCh=0;
    for(iC1=0;iC1<iChans;iC1++) if(pbadc[iC1]) cntBadCh++;
    pprctbadch[idx]= (double)cntBadCh/(double)iChans;
  }

  dRet = 1.0;
MCLEANUP:
 if(pbadc) free(pbadc);
 return dRet;
ENDVERBATIM
}

: get full time series correlation matrix for each time window into list of vectors
: tscorfull( 1 cor vec list, 2 time series list, 3 winsz, 4 incsz, 5 phzlist, 6 psatlist, 7 pderivlist)
FUNCTION tscorfull () {
  VERBATIM

  int i;

  double dRet = 0.0;
  ListVec* pTS = 0, *plvcor = 0, *plvhz = 0, *plvsat = 0, *plvderiv = 0;

  double** pTSN = 0; // normalized time series
  double* pdiff = 0; // for straight line test
  double* psums = 0, *psums2 = 0, *pfft = 0;

  int corsz=0; //# of vectors in output correlation vector list
  plvcor = AllocListVec(*hoc_objgetarg(1)); //each vec is correlation mat in vector form for a time point
  if(!plvcor || (corsz=plvcor->isz)<1){
    printf("tscorfull ERRA: problem initializing arg1 plvcor!\n");
    goto FCLEANUP;
  }

  pTS = AllocListVec(*hoc_objgetarg(2)); //list of eeg time series  
  if(!pTS){
    printf("tscorfull ERRA: problem initializing arg2 pTS!\n");
    goto FCLEANUP;
  }

  int iChans = pTS->isz; // # of channels
  if(iChans < 2) {
    printf("tscorfull ERRB: must have at least 2 timeseries!\n");
    goto FCLEANUP;
  }

  int iMatSz = iChans*(iChans-1)/2;

  for(i=0;i<corsz;i++){   //make sure all output vectors have right size
    if(plvcor->plen[i]<iMatSz){
      printf("tscorfull ERRC: each cor vec must have sz >= %d, [%d] has only %d\n",iMatSz,i,plvcor->plen[i]);
      goto FCLEANUP;
    }
  }

  int iWinSz = (int) *getarg(3); //window size for analysis
  int iINC = ifarg(4)?(int)*getarg(4):1;  //increment size for analysis

  plvhz = AllocListVec(*hoc_objgetarg(5));
  plvsat = AllocListVec(*hoc_objgetarg(6));
  plvderiv = AllocListVec(*hoc_objgetarg(7));

  pfft = (double*)calloc(iWinSz+1,sizeof(double));
  pdiff = (double*)calloc(iWinSz,sizeof(double));

  int iSamples = pTS->plen[0]; //make sure all time series have same length
  int iOutSz = (iSamples - iWinSz + 1) / iINC;//output List size

  if(verbose) printf("iWinSz=%d iINC=%d iSamples=%d\n",iWinSz,iINC,iSamples);
  for(i=1;i<iChans;i++) {
    if(pTS->plen[i] != iSamples){
      printf("tscorfull ERRD: time series of unequal size %d %d!\n",iSamples,pTS->plen[i]);
      goto FCLEANUP;
    }
  }

  if(corsz < iOutSz){
    printf("tscorfull ERRE: need output list vector of at least sz %d as arg 1, have only %d!\n",iOutSz,corsz);
    goto FCLEANUP;
  }   if(verbose) printf("iOutSz=%d\n",iOutSz);
  if(plvhz->isz<iOutSz || plvsat->isz<iOutSz || plvderiv->isz<iOutSz){
    printf("tscorfull ERRF: need output plvhz,plvsat,plvderiv size of at least %d!\n",iOutSz);
    goto FCLEANUP;
  }
  for(i=0;i<iOutSz;i++){
    if(plvhz->plen[i]<iChans || plvsat->plen[i]<iChans || plvderiv->plen[i]<iChans){
      printf("tscorfull ERRG: each vec in plvhz,plvsat,plvderiv must be sz >= %d!\n",iChans);
      goto FCLEANUP;
    }
  }

  pTSN = getdouble2D(iChans,iWinSz); if(verbose) printf("got pTSN\n"); //normalized time series
  if(!pTSN){
    printf("tscor ERRD: out of memory!\n");
    goto FCLEANUP;
  }

  psums = (double*)calloc(iChans,sizeof(double));   //for sum(X)
  psums2 = (double*)calloc(iChans,sizeof(double));  //for sum(X^2)
  double* pc = 0;
  double dVal = 0.0;
  int iStartIDX = 0 , ierr = 0 , iOutIDX = 0;
  FillListVec(plvcor,0.0); FillListVec(plvhz,0.); FillListVec(plvsat,0.); FillListVec(plvderiv,0.);//zero results
  for(iStartIDX=0;iOutIDX<iOutSz && iStartIDX+iWinSz<iSamples;iStartIDX+=iINC,iOutIDX++){
    if(iOutIDX%1000==0) printf("%d\n",iOutIDX);
    int IDX = iStartIDX , iChan = 0 , iEndIDX = iStartIDX + iWinSz;
    int iCurSz = iEndIDX - iStartIDX;
    for(iChan=0;iChan<iChans;iChan++){ //get sum(X) , sum(X^2) for ALL channels
      psums[iChan]=psums2[iChan]=0.;
      pc = &pTS->pv[iChan][iStartIDX];
      for(IDX=iStartIDX;IDX<iEndIDX;IDX++,pc++){
        psums[iChan] += *pc; 
        psums2[iChan] += (*pc * *pc);
      }
    }
    for(iChan=0;iChan<iChans;iChan++){       //do noise checks and normalize for ALL channels
      i=0;
      double davg = psums[iChan] / iCurSz;
      double dstdev = psums2[iChan]/iCurSz - davg*davg; 
      if(dstdev > 0. ) dstdev = sqrt(dstdev); else dstdev=1.; //make sure no problems with nan
      if(dstdev <= 0.) dstdev = 1.0;
      //saturation test
      int iSat = 0;
      pc = &pTS->pv[iChan][iStartIDX];
      for(IDX=iStartIDX;IDX<iEndIDX;IDX++,i++,pc++){	
        if((*pc>=900. && *pc<=1010.) || (*pc>=-1010. && *pc<=-900.)) iSat++;
        pTSN[iChan][i] = (*pc - davg) / dstdev;	//normalization
      }
      plvsat->pv[iOutIDX][iChan] = (double)iSat/ (double)iWinSz;
      //60,120 Hz Frequency test
      int fftlen = 0; memset(pfft,0,sizeof(double)*(iWinSz+1));
      if(dfftpow(&pTS->pv[iChan][iStartIDX],iWinSz,pfft,iWinSz+1,&fftlen)){
        int f = 0; pfft[0]=0.0; double fsum = 0.0;
        for(f=0;f<fftlen && f<=1000;f++) fsum += pfft[f]; //normalize power btwn 0-1000Hz to 1.0
        pc = &plvhz->pv[iOutIDX][iChan];
        for(f=55;f<=65 && f<fftlen;f++)   *pc += pfft[f]; //60Hz contribution
        for(f=115;f<=125 && f<fftlen;f++) *pc += pfft[f]; //120Hz contribution
        *pc /= fsum;
      }
      //discrete derivative test to measure how flat of a line signal looks like
      ddiff(&pTS->pv[iChan][iStartIDX],pdiff,iWinSz);
      plvderiv->pv[iOutIDX][iChan] = (double)iinrange(pdiff, iWinSz-1, -1.0, 1.0) / (double)(iWinSz-1.0);
    }

    //form correlation vector
    int iC1,iC2; pc = &plvcor->pv[iOutIDX][0];
    for(iC1=0;iC1<iChans;iC1++){
      for(iC2=iC1+1;iC2<iChans;iC2++){
	for(i=0;i<iWinSz;i++) 
          *pc += pTSN[iC1][i]*pTSN[iC2][i];
        *pc /= (double) iWinSz;
        pc++;
      }
    }
  }
  dRet = 1.0;
  printf("\n");

FCLEANUP:
  if(pTS) FreeListVec(&pTS);
  if(pTSN) freedouble2D(&pTSN,iChans);
  if(psums) free(psums);
  if(psums2) free(psums2);
  if(pfft) free(pfft);
  if(pdiff) free(pdiff);
  return dRet;
  ENDVERBATIM
}

: doesnt use eigenvalues 
FUNCTION tscoravg () {
  VERBATIM

  int i;

  double dRet = 0.0;
  ListVec* pTS = 0;

  double** pTSN = 0; // normalized time series
  double** pCorrel = 0; //correlation matrix

  double* vTCor = 0; //correlation matrix abs sum avg for a time point
  int corsz = vector_arg_px(1,&vTCor);

  pTS = AllocListVec(*hoc_objgetarg(2)); //list of eeg time series
  
  if(!pTS){
    printf("tscoravg ERRA: problem initializing 1st 2 args!\n");
    goto CLEANUP;
  }

  int iChans = pTS->isz; // # of channels
  if(iChans < 2) {
    printf("tscoravg ERRB: must have at least 2 EEG channels!\n");
    goto CLEANUP;
  }
  int iWinSz = (int) *getarg(3); //window size for analysis


  int iSamples = pTS->plen[0]; //make sure all time series have same length

  int iINC = ifarg(4)?(int)*getarg(4):1; 

  double *phz = 0, *psat = 0, *pfft = 0; int sz = 0 , iNoiseTest = 0;
  if(ifarg(5) && ifarg(6)){ //do saturation and 60/120 Hz Frequency test
    if((sz=vector_arg_px(5,&phz))!=corsz){
      printf("tscoravg ERRG: invalid size for phz should be %d, is %d!\n",corsz,sz);
      goto CLEANUP;
    } else if((sz=vector_arg_px(6,&psat))!=corsz){
      printf("tscoravg ERRH: invalid size for psat should be %d, is %d!\n",corsz,sz);
      goto CLEANUP;
    } else {	      
      iNoiseTest = 1;
      pfft = (double*)calloc(iWinSz+1,sizeof(double));
    }
  }

  int idoabs = ifarg(7)?(int)*getarg(7):0; //do abs or not

  int iOutSz = (iSamples - iWinSz + 1) / iINC;//output List size

  double *pdelta=0,*ptheta=0,*pbeta=0,*pgamma=0,*pripple=0;
  if(ifarg(8) && vector_arg_px(8,&pdelta)<iOutSz){ printf("bad pdelta size\n"); goto CLEANUP; } //delta 0-3
  if(ifarg(9) && vector_arg_px(9,&ptheta)<iOutSz){ printf("bad ptheta size\n"); goto CLEANUP; } //theta 4-12
  if(ifarg(10)&& vector_arg_px(10,&pbeta)<iOutSz){ printf("bad pbeta size\n"); goto CLEANUP; } //beta 13-29
  if(ifarg(11)&& vector_arg_px(11,&pgamma)<iOutSz){ printf("bad pgamma size\n"); goto CLEANUP; } //gamma 30-100
  if(ifarg(12)&& vector_arg_px(12,&pripple)<iOutSz){ printf("bad pripple size\n"); goto CLEANUP; } //ripple 101-300

  if(verbose) printf("iWinSz=%d iINC=%d iSamples=%d\n",iWinSz,iINC,iSamples);
  for(i=1;i<iChans;i++) {
    if(pTS->plen[i] != iSamples){
      printf("tscoravg ERRC: time series of unequal size %d %d!\n",iSamples,pTS->plen[i]);
      goto CLEANUP;
    }
  }

  if(corsz < iOutSz){
    printf("tscoravg ERRD: need output vector of at least %d as arg 1, have only %d!\n",iOutSz,corsz);
    goto CLEANUP;
  }
  if(verbose) printf("iOutSz=%d\n",iOutSz);

  pCorrel = getdouble2D(iChans,iChans); if(verbose) printf("got pCorrel\n");
  pTSN = getdouble2D(iChans,iWinSz); if(verbose) printf("got pTSN\n");
  if(!pCorrel || !pTSN){
    printf("tscor ERRD: out of memory!\n");
    goto CLEANUP;
  }

  double* psums = (double*)calloc(iChans,sizeof(double));
  double* psums2 = (double*)calloc(iChans,sizeof(double));

  double dVal = 0.0;
  int iStartIDX = 0 , ierr = 0 , iOutIDX = 0;
  for(iOutIDX=0;iOutIDX<corsz;iOutIDX++) vTCor[iOutIDX] = 0.;
  if( phz ) for(iOutIDX=0;iOutIDX<corsz;iOutIDX++) phz[iOutIDX]=0.;
  if( psat ) for(iOutIDX=0;iOutIDX<corsz;iOutIDX++) psat[iOutIDX]=0.;
  iOutIDX=0;
  for(iStartIDX=0;iOutIDX<iOutSz && iStartIDX+iWinSz<iSamples;iStartIDX+=iINC,iOutIDX++){
    if(iOutIDX%1000==0) printf("%d\n",iOutIDX);
    int IDX = iStartIDX , iChan = 0 , iEndIDX = iStartIDX + iWinSz;
    int iCurSz = iEndIDX - iStartIDX;
    if(phz) phz[iOutIDX]=0.; if(psat)psat[iOutIDX]=0.;
    if(pdelta) pdelta[iOutIDX]=0.; if(ptheta) ptheta[iOutIDX]=0.;
    if(pbeta) pbeta[iOutIDX]=0.; if(pgamma) pgamma[iOutIDX]=0.; if(pripple) pripple[iOutIDX]=0.;
    if(iStartIDX==0 || iINC > 1){
      for(iChan=0;iChan<iChans;iChan++){ //normalize data
	psums[iChan]=psums2[iChan]=0.;
	for(IDX=iStartIDX;IDX<iEndIDX;IDX++){
	  dVal = pTS->pv[iChan][IDX];
	  psums[iChan] += dVal;
	  psums2[iChan] += dVal*dVal;
	}
      }
    } else {
      for(iChan=0;iChan<iChans;iChan++){
	dVal = pTS->pv[iChan][iStartIDX-1];
	psums[iChan] -= dVal;
	psums2[iChan] -= dVal*dVal;
	dVal = pTS->pv[iChan][iEndIDX-1];
	psums[iChan] += dVal;
	psums2[iChan] += dVal*dVal;
      }
    }
    double dsatavg = 0., dhzavg = 0.;
    for(iChan=0;iChan<iChans;iChan++){       
      i=0;
      double davg = psums[iChan] / iCurSz;
      double dstdev = psums2[iChan]/iCurSz - davg*davg; 
      if(dstdev > 0. ) dstdev = sqrt(dstdev); else dstdev=1.; //make sure no problems with nan
      if(dstdev <= 0.) dstdev = 1.0;
      if(iNoiseTest){	
	//saturation test
	int iSat = 0;
	for(IDX=iStartIDX;IDX<iEndIDX;IDX++,i++){	
	  dVal = pTS->pv[iChan][IDX];
	  if((dVal>=900. && dVal<=1010.) || (dVal>=-1010. && dVal<=-900.)) iSat++;
	  pTSN[iChan][i] = (dVal - davg) / dstdev;	
	}
	dsatavg += (double)iSat/ (double)iWinSz;
	//60,120 Hz Frequency test
	int fftlen = 0; memset(pfft,0,sizeof(double)*(iWinSz+1));
	if(dfftpow(&pTS->pv[iChan][iStartIDX],iWinSz,pfft,iWinSz+1,&fftlen)){
	  int f = 0; pfft[0]=0.0; double fsum = 0.0;
	  for(f=0;f<fftlen && f<=1000;f++) fsum += pfft[f]; //normalize power btwn 0-1000Hz to 1.0
	  for(f=55;f<=65 && f<fftlen;f++) dhzavg += pfft[f] / fsum; //60Hz contribution
	  for(f=115;f<=125 && f<fftlen;f++) dhzavg += pfft[f] / fsum; //120Hz contribution
	  if(pdelta) {
	    for(f=0;f<=3;f++) pdelta[iOutIDX] += pfft[f]/fsum;
	    if(ptheta) for(f=4;f<=12;f++) ptheta[iOutIDX] += pfft[f]/fsum;
	    if(pbeta) for(f=13;f<=29;f++) pbeta[iOutIDX] += pfft[f]/fsum;
	    if(pgamma) for(f=30;f<=100;f++) pgamma[iOutIDX] += pfft[f]/fsum;
	    if(pripple) for(f=101;f<=300 && f<fftlen;f++) pripple[iOutIDX] += pfft[f]/fsum;
	  }
	}	
      } else {
	for(IDX=iStartIDX;IDX<iEndIDX;IDX++,i++){	
	  dVal = pTS->pv[iChan][IDX];
	  pTSN[iChan][i] = (dVal - davg) / dstdev;	
	}
      }
    }
    if(iNoiseTest){
      phz[iOutIDX] = dhzavg / (double)iChans;
      psat[iOutIDX] = dsatavg / (double)iChans;
      if(pdelta) pdelta[iOutIDX] /= (double)iChans;
      if(ptheta) ptheta[iOutIDX] /= (double)iChans;
      if(pbeta) pbeta[iOutIDX] /= (double)iChans;
      if(pgamma) pgamma[iOutIDX] /= (double)iChans;
      if(pripple) pripple[iOutIDX] /= (double)iChans;
    }

    //form correlation matrix
    int iC1,iC2;
    for(iC1=0;iC1<iChans;iC1++){
      for(iC2=0;iC2<=iC1;iC2++){
	if(iC1==iC2){
	  pCorrel[iC1][iC2]=1.0; //diagonals == 1.0
	  continue;
	} else pCorrel[iC1][iC2]=0.; //init to 0.
	for(i=0;i<iWinSz;i++){
	  pCorrel[iC1][iC2] += pTSN[iC1][i]*pTSN[iC2][i];
	}
	pCorrel[iC1][iC2] /= (double)iWinSz;
	pCorrel[iC2][iC1] = pCorrel[iC1][iC2];
      }
    }
    double dpsum = 0.0; int icc = 0; //get avg of abs val of correlations
    if(idoabs){
      for(iC1=0;iC1<iChans;iC1++){
	for(iC2=0;iC2<iC1;iC2++){
	  dpsum += fabs(pCorrel[iC1][iC2]);
	  icc++;
	}
      }
    } else {
      for(iC1=0;iC1<iChans;iC1++){
	for(iC2=0;iC2<iC1;iC2++){
	  dpsum += pCorrel[iC1][iC2];
	  icc++;
	}
      }
    }
    dpsum /= (double) icc;
    if(dpsum > 1. || dpsum < -1.){
      printf("out of bounds = %g!\n",dpsum);
    }
    vTCor[iOutIDX] = dpsum;
  }
  dRet = 1.0;
  printf("\n");

CLEANUP:
  if(pTS) FreeListVec(&pTS);
  if(pCorrel) freedouble2D(&pCorrel,iChans);
  if(pTSN) freedouble2D(&pTSN,iChans);
  free(psums);
  free(psums2);
  if(pfft) free(pfft);
  return dRet;
  ENDVERBATIM
}


FUNCTION tscor () {
  VERBATIM

  int i;

  double dRet = 0.0;

  ListVec* pEIG = 0, *pTS = 0;

  double** pTSN = 0; // normalized time series
  double** pCorrel = 0; //correlation matrix
  double** pV = 0;
  double* pW = 0; // eigs

  pEIG = AllocListVec(*hoc_objgetarg(1)); //list of eigenvalue vectors, one for each time
  pTS = AllocListVec(*hoc_objgetarg(2)); //list of eeg time series
  
  if(!pEIG || !pTS){
    printf("tscor ERRA: problem initializing 1st 2 args!\n");
    goto CLEANUP;
  }

  int iChans = pTS->isz; // # of channels
  if(iChans < 2) {
    printf("tscor ERRB: must have at least 2 EEG channels!\n");
    goto CLEANUP;
  }
  int iWinSz = (int) *getarg(3); //window size for analysis
  int iINC = (int) *getarg(4); //increment for analysis

  int iSamples = pTS->plen[0]; //make sure all time series have same length

  if(verbose) printf("iWinSz=%d iINC=%d iSamples=%d\n",iWinSz,iINC,iSamples);
  for(i=1;i<iChans;i++) {
    if(pTS->plen[i] != iSamples){
      printf("tscor ERRC: time series of unequal size %d %d!\n",iSamples,pTS->plen[i]);
      goto CLEANUP;
    }
  }
  int iOutSz = (iSamples - iWinSz + 1) / iINC;//output List size
  if(pEIG->isz < iOutSz){
    printf("tscor ERRD: need List of size at least %d as arg 1, have only %d!\n",iOutSz,pEIG->isz);
    goto CLEANUP;
  }
  if(verbose) printf("iOutSz=%d\n",iOutSz);

  pCorrel = getdouble2D(iChans,iChans); if(verbose) printf("got pCorrel\n");
  pTSN = getdouble2D(iChans,iWinSz); if(verbose) printf("got pTSN\n");
  pV = getdouble2D(iChans,iChans);   if(verbose) printf("got pV\n");
  pW = (double*) malloc(sizeof(double)*iChans); if(verbose) printf("got pW\n");
  if(!pCorrel || !pTSN || !pV || !pW){
    printf("tscor ERRD: out of memory!\n");
    goto CLEANUP;
  }

  double dVal = 0.0;
  int iStartIDX = 0 , ierr = 0 , iTIDX = 0;
  for(iStartIDX=0;iStartIDX+iWinSz<iSamples;iStartIDX+=iINC,iTIDX++){
    if(iStartIDX%100==0) printf("%d\n",iStartIDX);
    int IDX = iStartIDX , iChan = 0 , iEndIDX = iStartIDX + iWinSz;
    int iCurSz = iEndIDX - iStartIDX;
    for(iChan=0;iChan<iChans;iChan++){ //normalize data
      double dSum = 0.0, dSum2 = 0.0;
      for(IDX=iStartIDX;IDX<iEndIDX;IDX++){
        dVal = pTS->pv[iChan][IDX];
	dSum += dVal;
	dSum2 += dVal*dVal;
      }
      dSum /= iCurSz; //mean
      dSum2 /= iCurSz; 
      dSum2 -= dSum*dSum;
      dSum2 = sqrt(dSum2); //standard deviation
      i=0;
      for(IDX=iStartIDX;IDX<iEndIDX;IDX++,i++){
	pTSN[iChan][i] = (pTS->pv[iChan][IDX] - dSum) / dSum2;
      }
    }
    //form correlation matrix
    int iC1,iC2;
    for(iC1=0;iC1<iChans;iC1++){
      for(iC2=0;iC2<=iC1;iC2++){
	if(iC1==iC2){
	  pCorrel[iC1][iC2]=1.0; //diagonals == 1.0
	  continue;
	} else pCorrel[iC1][iC2]=0.;
	for(i=0;i<iWinSz;i++){
	  pCorrel[iC1][iC2] += pTSN[iC1][i]*pTSN[iC2][i];
	}
	pCorrel[iC1][iC2] /= (double)iWinSz;
	pCorrel[iC2][iC1] = pCorrel[iC1][iC2];
      }
    }
    //perform SVD to get eigenvalues into pW
    mysvd(iChans,iChans,pCorrel,pW,pV,&ierr); 
    qsort(pW,iChans,sizeof(double),mycompare); //sort eigenvalues in ascending order
//    for(i=0;i<iChans;i++) pEIG->pv[iTIDX][i] = pW[i]; //store sorted eigenvalues
    memcpy(pEIG->pv[iTIDX],pW,sizeof(double)*iChans); //store sorted eigenvalues
  }
  dRet = 1.0;
  printf("\n");

CLEANUP:
  if(pEIG) FreeListVec(&pEIG);
  if(pTS) FreeListVec(&pTS);
  if(pCorrel) freedouble2D(&pCorrel,iChans);
  if(pTSN) freedouble2D(&pTSN,iChans);
  if(pW) free(pW);
  if(pV) freedouble2D(&pV,iChans);
  return dRet;
  ENDVERBATIM
}

VERBATIM

double getmean(double* p,int n) {
  if(n<1) return 0.0;
  int i = 0;
  double sum = 0.0, *pp=p;
  for(i=0;i<n;i++,pp++) sum += *pp;
  return sum / (double) n;
}

double getstd(double* p,int n,double X) {
  if(n<1) return 0.;
  int i;
  double X2=0., *pp=p;
  for(i=0;i<n;i++,pp++) X2 += (*pp * *pp);
  X2 = X2/(double)n - X*X;
  if(X2>0.) return sqrt(X2);
  return 0.;
}

double dnormv(double* p,int n) {
  if(n<1) return 0.0;
  double X = getmean(p,n);
  double s = getstd(p,n,X);
  double* pp = p;
  int i;
  if(s>0.)
    for(i=0;i<n;i++,pp++) *pp = (*pp - X)/s; 
  else
    for(i=0;i<n;i++,pp++) *pp -= X;
  return 1.0;
}

static double normv(void* v){ 
  double* x;
  int sz = vector_instance_px(v,&x);
  if(sz<1){
    printf("normv ERRA: empty input size!\n");
    return 0.;
  }
  return dnormv(x,sz);
}

double dcopynzidx(double* pin,double* pidx,double* pout,int sz) {
  int szout,i;
  szout=0;
  for(i=0;i<sz;i++) if(pidx[i]) pout[szout++]=pin[i];
  return (double) szout;
}

static double copynzidx (void* v) {
  double *x,*y,*z,ret; int sz;
  sz = vector_instance_px(v,&x);
  if(sz!=vector_arg_px(1,&y)) y=vector_newsize(vector_arg(1),sz);
  if(sz!=vector_arg_px(2,&z)) z=vector_newsize(vector_arg(2),sz);
  ret=dcopynzidx(x,y,z,sz);
  vector_resize(vector_arg(2),(int)ret);
  return ret;
}
ENDVERBATIM

: assumes v1,v2 are normalized
FUNCTION tsmul () {
  VERBATIM
  double* x, *y, *px, *py, sum = 0.;
  int xsz,ysz,i;
  xsz = vector_arg_px(1,&x);
  ysz = vector_arg_px(2,&y);
  px=x; 
  py=y;
  if(xsz!=ysz || xsz<1 || ysz<1){
    printf("tsmul ERRB: input vecs are invalid sizes: %d %d!\n",xsz,ysz);
    return -2.;
  }
  for(i=0;i<xsz;i++,px++,py++) sum += *px * *py;
  return sum / (double) xsz;
  ENDVERBATIM
}

PROCEDURE install () {
  if (INSTALLED==1) {
    printf("already installed tsa.mod")
  } else {
    INSTALLED=1
    VERBATIM
      //    install_vector_method("myspct", myspct);
    install_vector_method("inrange",inrange);
    install_vector_method("normv",normv);
    install_vector_method("copynzidx",copynzidx);
    ENDVERBATIM
  }
}