/* Standalone analysis program, reading fixed-point monaural sound files.
 * Currently set for maximum of 34 poles, & max anal segment of 500 samples,
 * meaning that the frame slices cannot exceed 250 samples.
 * Program size expands linearly with slice size, and as square of npoles.
 */

#include "../../sysdep.h"
#include <stdio.h>
#include <math.h>    /* removed include <sys/file.h> */
#include <pwd.h>

#ifdef SFIRCAM
#include <local/sfheader.h>
#endif SFIRCAM
#include "crack.h"
#include "lpsf.h"
#include "../../lpc.h"

#define INCR 4
#define NAMESIZE 1024		/* commandline string size limit */

static int NPOLE;
static int FRAME;
static int NPP1;

char	*buildsfname();
extern	int	remotein;
int debug=0, verbose=0;

#define POLEMAX	34
#define FRAMAX  500 
#define NDATA 	4	/* number of data values stored with frame */

main(argc, argv)
     int argc;
     char **argv;
{
        int	jj, slice, counter;
        float	coef[POLEMAX+NDATA], inskip, dur;
        char	*sfname;
	char    iname[NAMESIZE];
	char    oname[NAMESIZE];
        double	errn, rms1, rms2, cc[40];
        short	sigi[FRAMAX];			/* changed to short */
        long	i, nsamps, skipbytes, outskip;
        struct	SFDESC sfdesc;
        int	anal, sfd, nbytes, nblpc, firstframe, nread;

	int	n, hsize, esize;		/* new declarations */
#ifdef SFIRCAM
	char	inbuf[SIZEOF_HEADER];
	SFHEADER *sfh;
	int Hflag = 0;
#endif SFIRCAM
	LPHEADER *lph;
	char	*lpbuf, *tp;
	char    c;
	int rflag = 0;

	lpbuf = (char *) calloc(LPBUFSIZ, sizeof(char));
	lph = (LPHEADER *) lpbuf;
	tp = lph->text;

	/* DEFAULTS... */
	*iname = '\0';		/* '\0' ==> stdin */
	*oname = '\0';		/* '\0' ==> stdout */
	NPOLE = 34;	
	slice = 200;
	inskip = 0.0;
	dur = 1.0;		/* should implement an EOF dur */
	*tp = '\0';
	sfdesc.sf_chans = 1;
	sfdesc.sf_class = INT;
	sfdesc.sf_srate = 20000;

	if (debug) verbose++;
	if ( (c = crack(argc, argv, "q", 1)) == 'q' ) {	/* q flag present */
	/* query for everything, ignore any other flags */
	  verbose++;
	  fprintf(stderr, "Enter name of soundfile to analyze : ");
	  gets(iname);
#ifdef SFIRCAM
	  fprintf(stderr,"Is soundfile headerless? (y or n) ");
	  gets(oname); 
	  if ( *oname == 'y' || *oname == 'Y' ) {
	    Hflag++; fprintf(stderr,"Enter soundfile sampling rate : ");
	    scanf("%d", &(sfdesc.sf_srate)); rflag++;
	  }
#endif SFIRCAM
	  fprintf(stderr, "Enter name of output data file : ");
	  gets(oname);
	  fprintf(stderr, "Enter number of poles : ");
	  scanf("%d", &NPOLE);
	  fprintf(stderr, "Enter inter-frame offset in samples : ");
	  scanf("%d", &slice);
	  fprintf(stderr, "Specify input skip (secs) : ");
	  scanf("%f", &inskip);
	  fprintf(stderr, "Specify duration to analyze (secs) : ");
	  scanf("%f", &dur);
	  fflush(stdin);
	  fprintf(stderr, "Enter any useful facts about the file.\n");
	  fprintf(stderr, "Type 2 carriage returns to quit.\n");
	  fprintf(stderr, "> ");
#ifndef SYS5
	  getchar();			/* VAX, etc */
#endif
	start:				/* collect text until 2 cr's */
	  while ((*tp++ = getchar()) != '\n')	
	    if ((tp - lph->text) == (LPBUFSIZ - sizeof(LPHEADER) + 4))
	      goto foo;
	  fprintf(stderr, "> ");
	  if ((*tp++ = getchar()) != '\n') goto start;
	} 
	else if ( c == EOF ) usage("-q flag doesn't take an option");
	else {			/* no q flag present */
	  /* set options from commandline */
	  arg_index = 0;	/* re-crack */
	  while ((c = crack(argc, argv, 
#ifdef SFIRCAM
			    "r|p|i|s|d|o|C|H", 
#else SFIRCAM
			    "r|p|i|s|d|o|C|", 
#endif SFIRCAM
			    0)) != NULL) {
	    if (c == EOF) usage("error cracking commandline");
	    switch (c) {	/* cases w/o options */
#ifdef SFIRCAM
	    case 'H': Hflag++; break; /* no soundfile header */
#endif SFIRCAM	
	    default:		/* cases w/ options */
	      if (arg_option == NULL) usage("flag missing an option");
	      switch (c) {
		/* sound flags... */
	      case 'r':		/* sampling rate */
		sscanf(arg_option,"%d", &(sfdesc.sf_srate)); rflag++; break;
		/* analysis and i/o flags... */
	      case 'p':		/* npoles */
		sscanf(arg_option,"%d", &NPOLE); break;
	      case 'i':		/* frame size */
		sscanf(arg_option,"%d", &slice); break; 
	      case 's':		/* skiptime */
		sscanf(arg_option,"%f", &inskip); break;	
	      case 'd':		/* duration */
		sscanf(arg_option,"%f", &dur); break; 
	      case 'o':		/* output filename */
		if ( strlen(arg_option) >= NAMESIZE )
		  die("outputfile name too long!");
		strcpy(oname,arg_option); break; 
	      case 'C':		/* comment string */
		strncat(tp,arg_option,(LPBUFSIZ - sizeof(LPHEADER) + 4));
		tp += strlen(tp);
		break;
	      default: usage("unrecognized flag");
	      }
	    }
	  }
	  if ( arg_index != argc ) { /* soundfile name here */
	    if ( strlen(argv[arg_index]) >= NAMESIZE )
	      die("soundfile name too long!");
	    strcpy(iname,argv[arg_index]);
	  }
	}
	if (verbose) {
	  fprintf(stderr,"Writing to lpfile ");
	  if (*oname == '\0') fprintf(stderr,"<stdout>\n");
	  else fprintf(stderr,"%s\n",oname);
	  fprintf(stderr,"Reading from soundfile ");
	  if (*iname == '\0') fprintf(stderr,"<stdin>\n");
	  else fprintf(stderr,"%s\n",iname);
	  fprintf(stderr, "poles=%d interframeoffset=%d skip=%f duration=%f\n",
		  NPOLE, slice, inskip, dur);
	  fprintf(stderr, "lpheader comment=\n%s\n", lph->text);
	}
	/* finished with comandline parsing/querying */
      foo:
	if ( *oname == '\0' ) anal = fileno(stdout); /* write stdout */
	else {
	  if((anal = creat(oname,0644)) < 0)
	    die("cannot open output data file");
	}
	if ( *iname == '\0' ) sfd = fileno(stdin); /* read stdin */
	else {
	  sfname = buildsfname(iname);
#ifdef SFREMOTE
	  if (index(sfname,':') != NULL)
	    sfd = rsfopen(sfname,'<');
	  else 
#endif SFREMOTE
	    if ((sfd = open(sfname,O_RDONLY)) < 0)
	      dies("cannot open sound file %s", sfname);
	}

#ifdef SFIRCAM
	if ( Hflag == 0 ) {	/* soundfile header present */
	  if ((n = readin(sfd, inbuf, SIZEOF_HEADER)) < SIZEOF_HEADER)
	    die("soundfile header read error");
	
	  sfh = (SFHEADER *) inbuf;
	  if (sfmagic(sfh) != SF_MAGIC)
	    dies("%s is not a soundfile", iname);
	
	  sfdesc.sf_chans = sfchans(sfh);
	  sfdesc.sf_class = sfclass(sfh);
	  if ( rflag == 0 ) sfdesc.sf_srate = sfsrate(sfh);
	}
#endif SFIRCAM
	if (debug) {
	  fprintf(stderr, "sf_chans=%d sf_srate=%d sf_class=%d\n", 
		  sfdesc.sf_chans, sfdesc.sf_srate, sfdesc.sf_class);
	}
/* 	End new code		 */

	if(sfdesc.sf_chans != 1)
		die("takes monaural files only");
	if(sfdesc.sf_class != INT)
		die("takes integer (shortsam) files only");

	if (NPOLE > POLEMAX)
		dies("poles exceeds maximum allowed");
	FRAME = slice * 2;
	if (FRAME > FRAMAX)
	  die("framesize (inter-frame-offset*2) exceeds maximum allowed");

	lph->lpmagic = LP_MAGIC;
	lph->srate = sfdesc.sf_srate;
	lph->npoles = NPOLE;
	lph->nvals = NPOLE + 4;
	lph->framrate = lph->srate / (float) slice;
	lph->duration = dur;

	/* compute header size */
	hsize = tp - (char *) lph;	/* sizeof text */
	esize = INCR - (hsize % INCR);	/* round up to 4 byte boundary */
	if (esize == 4) esize = 0;	/* if we need it */
	lph->headersize = hsize + esize;
	if (lph->headersize >= (LPBUFSIZ - sizeof(LPHEADER) + 4))
	  lph->headersize = LPBUFSIZ;

/* 	More new code for lpheader stuff		 */
	if (write(anal, (char *) lph, lph->headersize) < lph->headersize)
		die("anallpc: can't write header");
/* 	End of lp header code	 */
	
	skipbytes = (long)(inskip * sfdesc.sf_srate) * (long)sfdesc.sf_class;
	if (remotein || *iname == '\0') { /* not seekable */
	  while (skipbytes > FRAMAX) { /* remote: spinread for lseek */
	    if ((nread = readin(sfd,(char *)sigi,FRAMAX)) < FRAMAX)
	      die("soundfile skip error");
	    skipbytes -= nread;
	  }
	  if ((nread = readin(sfd, (char *)sigi, skipbytes)) < skipbytes)
	    die("soundfile skip error");
	}
	else if (lseek(sfd,skipbytes,L_INCR) < 0)
		die("bad sflseek on skip");

/*	fprintf(stderr,"Specify starting output frame number\t");
/*	scanf("%d",&firstframe);			
/*  for now we'll just say to start at 0 ('C' mind set)  	*/

	firstframe = 0;

	NPP1 = NPOLE+1;
	nblpc = (NPOLE + NDATA)*FLOAT;
	outskip = firstframe * nblpc + lph->headersize;	/* headersize is new */

	if (*oname != '\0')
	  if ((lseek(anal,outskip,0)) < 0)
	    die("Bad lseek on analysis file");
	else if (firstframe != 0)
	  fprintf(stderr,"warning: can't lseek on stdout");


	nsamps = dur * sfdesc.sf_srate;
	nbytes = FRAME * sfdesc.sf_class;	/* bytes for whole frame */
	if ((nread = readin(sfd, (char *)sigi, nbytes)) != nbytes)
	  die("soundfile read error: couldn't fill first frame");
/*	for(i=0;i<FRAME;i++) printf("%d ",sigi[i]); 	*/

	if (debug) fprintf(stderr, "nbytes : %d\n", nbytes);
	nbytes = slice * sfdesc.sf_class; /* bytes for half frame */
				/* FRAME is slice * 2 */
	if (debug) fprintf(stderr, "nbytes shifted: %d\n", nbytes);
	i = 0;
	counter = 1;

	while (i < nsamps) {
		alpol(sigi,&errn,&rms1,&rms2,cc);
		coef[0] = (float)rms2; 
		coef[1] = (float)rms1; 
		coef[2] = (float)errn; 
		coef[3] = 0.;			  /* for pitch of this frame */

/* 		fprintf(stderr,"%d %f %f %f %f\n",
			counter++,coef[0],coef[1],coef[2],coef[3]);	 */
			
					  /* reverse the coefs & change sign */
		for(jj=NDATA; jj<NPOLE+NDATA; jj++)
			coef[jj] = (float)-cc[NPOLE-jj+(NDATA - 1)];  
						/* then write the anal frame */
		if ((n = write(anal, (char *)coef, nblpc)) != nblpc)
			die("write error");

		for(jj=0; jj<slice; jj++,i++)	   /* now move slice forward */
			sigi[jj] = sigi[jj+slice];	/*   & refill:	     */
		if ((n = readin(sfd, (char *)(sigi+slice), nbytes)) != nbytes)
		  break;	/* ran up against eof */
	}
}

alpol(sig, errn, rms1, rms2, c)
 double *errn, *rms1, *rms2, *c;
 short *sig;					/* sig now short */
{
	double a[POLEMAX][POLEMAX], v[POLEMAX], b[POLEMAX];
	double x[FRAMAX], y[FRAMAX];
	double *vp=v, *bp=b, *xp=x, *yp=y;
	double sum, sumx, sumy;
	int k1, i, l, k, limit, j;

	for (xp=x; xp-x < FRAME ;++xp,++sig) 
		*xp = (double) *sig;
	k1 = NPOLE + 1;
	for (i=0; i < NPOLE ;++i)  {
		sum = (double) 0.0;
		for (k=NPOLE; k < FRAME ;++k)
			sum += x[k-(i+1)] * x[k];
		v[i] = -sum;
		if (i != NPOLE - 1)  {
			limit = NPOLE - (i+1);
			for (l=0; l < limit ;++l)  {
				sum += x[NPOLE-(i+1)-(l+1)]* x[NPOLE-(l+1)] - x[FRAME-(i+1)-(l+1)]* x[FRAME-(l+1)];
				a[(i+1)+l][l] = a[l][(i+1)+l] = sum;
			}
		}
	}
	sum = (double) 0.0;
	for (k=NPOLE; k < FRAME ;++k)
		sum += pow(x[k], (double) 2.0);
	sumy = sumx = sum;
	for (l=0; l < NPOLE ;++l)  {
		sum += pow(x[NPOLE-(l+1)], (double) 2.0) - pow(x[FRAME-(l+1)], (double) 2.0);
		a[l][l] = sum;
	}
	gauss(a, v, b);
/*	filtn(x, y, b);   */
	for (i=0; i < NPOLE ;++i)
		sumy = sumy - b[i]*v[i];
	*rms1 = sqrt(sumx/(double) (FRAME - k1 + 1) );
	*rms2 = sqrt(sumy/(double) (FRAME - k1 + 1) );
	*errn = pow(((*rms2)/(*rms1)), (double) 2.0);
	for (bp=b; bp-b < NPOLE ;++bp,++c)
		*c = *bp;
	return(0);
}

gauss(aold, bold, b)
double aold[POLEMAX][POLEMAX], *bold, *b;
{
	double amax, dum, pivot;
	double c[POLEMAX], a[POLEMAX][POLEMAX];
	int i, j, k, l, istar, ii, lp;

	/* aold and bold untouched by this subroutine */
	for (i=0; i < NPOLE ;++i)  {
		c[i] = bold[i];
		for (j=0; j < NPOLE ;++j)
			a[i][j] = aold[i][j];
	}
	/* eliminate i-th unknown */
	for (i=0; i < NPOLE - 1 ;++i)  {        /* find largest pivot */
		amax = 0.0;
		for (ii=i; ii < NPOLE ;++ii)  {
			if (fabs(a[ii][i]) >= amax)  {
				istar = ii;
				amax = fabs(a[ii][i]);
			}
		}
		if (amax < 1e-20)
			die("gauss: ill-conditioned");

		for (j=0; j < NPOLE ;++j)  {    /* switch rows */
			dum = a[istar][j];
			a[istar][j] = a[i][j];
			a[i][j] = dum;
		}
		dum = c[istar];
		c[istar] = c[i];
		c[i] = dum;
		/* pivot */
		for (j=i+1; j < NPOLE ;++j)  {
			pivot = a[j][i] / a[i][i];
			c[j] = c[j] - pivot * c[i];
			for (k=0; k < NPOLE ;++k)
				a[j][k] = a[j][k] - pivot * a[i][k];
		}
	}       /* return if last pivot is too small */
	if (fabs(a[NPOLE-1][NPOLE-1]) < 1e-20)
		die("gauss: ill-conditioned");

	*(b+NPOLE-1) = c[NPOLE-1] / a[NPOLE-1][NPOLE-1];
	/* back substitute */
	for (k=0; k<NPOLE-1; ++k)  {
		l = NPOLE-1 -(k+1);
		*(b+l) = c[l];
		lp = l + 1;
		for (j = lp; j<NPOLE; ++j)
			*(b+l) += -a[l][j] * *(b+j);
		*(b+l) /= a[l][l];
	}
	return(1);
}


filtn(x, y, b)
double x[], b[], y[];
{
	double sum;
	int i, j;
	double *yp;

	for (yp=y; yp-y < NPOLE ;++yp)
		*yp = (double) 0.0;
	yp = y;
	for (i=NPOLE; i < FRAME; ++i)  {
		sum = x[i];
		for (j=0; j < NPOLE ;++j)  {
			sum = sum + b[j] * x[i-(j+1)];
		}
		*(yp+i) = sum;
	}
	return(0);
}

iarrayprint(array, size, name)
short array[];
int size;
char *name;
{
	 int i;
	 printf("\n%s : \n", name);
	 for (i = 0; i < size; i++)	{
		  if ((i % 10) == 0)
		    printf("\n");
		  printf("%d ", array[i]);
	 }
}

farrayprint(array, size, name)
 double array[];
 int size;
 char *name;
{
	 int i;
	 printf("\n%s : \n", name);
	 for (i = 0; i < size; i++)	{
		  if ((i % 10) == 0)
		    printf("\n");
		  printf("%.1f ", array[i]);
	 }
}
	
usage(msg)
     char *msg;
{
#ifdef SFIRCAM
  fprintf (stderr, "%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s",
#else SFIRCAM
  fprintf (stderr, "%s%s%s%s%s%s%s%s%s%s%s%s%s%s",
#endif SFIRCAM
   "   Usage:  anallpc [flag][option] ... [soundfile]\n",
   "[flag][option] from among:\n",
/*   "-cNCHANS          number of channels of sound input (default 1)\n",*/
#ifdef SFIRCAM
   "-rSRATE           sampling rate of sound input (default read from\n",
   "                     soundfile header; else 20000 if -H specified)\n",
#else
   "-rSRATE           sampling rate of sound input (default 20000)\n",
#endif
   "-pNPOLES          number of poles for analysis (default 34)\n",
   "-iINTERFRAMEOFSET offset between frames in samples (default 200)\n",
   "                     (framesize will be twice INTERFRAMEOFSET)\n",
   "-sSKIPTIME        initial seconds of sound to skip over (default 0.0)\n",
   "-dDURATION        duration in seconds to analyze (default 1.0)\n",
   "-oOUTPUTFILE      analysis output file (stdout if absent, default)\n",
   "-CCOMMENT_STRING  comment field of lp header (default empty)\n",
#ifdef SFIRCAM
   "-H                no soundfile header, SRATE taken from -r flag\n",
   "                       (default reads SRATE from soundfile header)\n",
#endif
   "soundfile         input (reads floats on stdin if absent, default)\n",
   "\n",
   "      or:  anallpc -q   (queries the user for commandline options)\n",
   "see also:  man anallpc\n\n"
   );
	   die(msg);
}