#include <stdio.h>
#include <stdlib.h>
#include <helper.h>
#include <pv.h>
#include <sound.h>
#include "burrowing.h"

int main(argc, argv)
    int argc; char **argv;
{
    int 	i,
		R = 44100,
		N = 1024,
		N2,
		Nw = 0,
		Nw2, 
		D = 0, 
		I = 0,
		in,
		on,
		eof = 0,
		useRms = 0,
		smoothWidth = 0,
		*bitshuffle,
		usage();

    float 	*Hwin,
		*Wanal,
		*Wsyn,
		*input1,
                *input2,
		*buffer1,
                *buffer2,
		*channel1,
                *channel2,
		*output,
		**smoothWindow,
		*smoothRms,
		rms = 0.,
		thresh = .001,
		currentThresh,
    		mult = .001,
		*trigland;

    char	ch,
		file1[1024],
                file2[1024];

    FILE	*fp1,
                *fp2;

    MrSoundStruct     soundStruct1,
		      soundStruct2;


    if (argc < 3)
	usage(1);

    synt = 0.0;

    while( (ch= crack( argc, argv, "R|N|M|D|I|f|F|m|t|w|rh", 0  )) != 0 ) {
	switch(ch) {
	    case 'R':	R = atoi(arg_option);
			break;
	    case 'N':	N = atoi(arg_option);
			break;
	    case 'M':	Nw = atoi(arg_option);
			break;
	    case 'D':	D = atoi(arg_option);
			break;
	    case 'I':	I = atoi(arg_option);
			break;
	    case 'f':	strcpy(file1, arg_option);
			break;
	    case 'F':	strcpy(file2, arg_option);
			break;
	    case 'm':   mult = atof(arg_option);
	                break;
	    case 't':	thresh = atof(arg_option);
			break;
	    case 'w':	smoothWidth = atof(arg_option);
			break;
	    case 'r':	useRms = 1;
			break;
	    case 'h':	usage(1);
	}
    }

    if (Nw == 0)
	Nw = N;

    if (D == 0)
      D = N / 2;

    if (I == 0)
      I = D;

    if (smoothWidth < 1)
      smoothWidth = 1;

    myPI = 4.*atan(1.);
    myTWOPI = 8.*atan(1.);
    N2 = N>>1;
    Nw2 = Nw>>1;

    Wanal = (float *) space( Nw, sizeof(float) );
    Wsyn = (float *) space( Nw, sizeof(float) );
    Hwin = (float *) space( Nw, sizeof(float) );
    input1 = (float *) zerospace( Nw, sizeof(float) );
    input2 = (float *) zerospace( Nw, sizeof(float) );
    output = (float *) zerospace( Nw, sizeof(float) );
    buffer1 = (float *) space( N, sizeof(float) );
    buffer2 = (float *) space( N, sizeof(float) );
    channel1 = (float *) space( N+2, sizeof(float) );
    channel2 = (float *) space( N+2, sizeof(float) );

    /* make space for magnitude smoothing function */

    smoothWindow = (float **) space( smoothWidth+1, sizeof(float *) );

    for (i=0; i < smoothWidth+1; i++)
      *(smoothWindow+i) = (float *) zerospace( N2+1, sizeof(float) );

    smoothRms = (float *) zerospace( N2+1, sizeof(float) );


    if ( (Nw / D) == 2 )
      makehanning( Hwin, Wanal, Wsyn, Nw, N, I, 0, 1 );

    else
      makehanning( Hwin, Wanal, Wsyn, Nw, N, I, 0, 0 );
    
  /* FFT cosine funk */
    trigland = (float *) space( N2, sizeof(float) );

  /* FFT bit shuffle */
    bitshuffle = (int *) space( 3 + (int) sqrt( (float) N2 ), sizeof(int) );

    init_rdft( N, bitshuffle, trigland );


    in = -Nw;
    if ( D )
	on = (in*I)/D;
    else
	on = in;

    /* open input files */

    if ( (fp1 = GetMrSoundStream( file1, &soundStruct1 )) == NULL ) {
      fprintf(stderr,"Cannot open %s\n",file1);
      exit(-1);
    }

    if ( soundStruct1.dataFormat != SND_FORMAT_FLOAT ) {
      fprintf(stderr,"%s is not a 32-bit floating point sound file\n", file1);
      exit(-1);
    }


    if ( (fp2 = GetMrSoundStream( file2, &soundStruct2 )) == NULL ) {
      fprintf(stderr,"Cannot open %s\n",file2);
      exit(-1);
    }

    if ( soundStruct2.dataFormat != SND_FORMAT_FLOAT ) {
      fprintf(stderr,"%s is not a 32-bit floating point sound file\n", file2);
      exit(-1);
    }


    while ( !eof ) {

      in += D;
      on += I;
      
      eof = floatin( input1, Nw, D, fp1 );
      if (eof)
	(void) floatin( input2, Nw, D, fp2 );
      else
	eof = floatin( input2, Nw, D, fp2 );


      currentThresh = thresh;

      if (useRms) {
	
	float	smoothed;

	rms = 0.;

	for ( i=0; i < Nw; i++ )
	  rms += *(input1+i) * *(input1+i);
	
	rms = sqrt( rms / Nw );

	smoothed = rms;

	for ( i=0; i < smoothWidth; i++ )
	  smoothed += *(smoothRms+i);

	smoothed /= (float) (smoothWidth + 1);
	
	for ( i=smoothWidth-1; i > 0; i-- )
	  *(smoothRms+(i-1)) = *(smoothRms+i);
	  
	*(smoothRms+(smoothWidth-1)) = rms;

	currentThresh *= smoothed;
      }

      fold( input1, Wanal, Nw, buffer1, N, in );
      fold( input2, Wanal, Nw, buffer2, N, in );
      rdft( N, 1, buffer1, bitshuffle, trigland );
      rdft( N, 1, buffer2, bitshuffle, trigland );

      burrowing( buffer1, buffer2, channel1, channel2, N2, mult, currentThresh,
		 smoothWidth, smoothWindow );
		        
      if ( I == 0 ) {
	fwrite( channel1, sizeof(float), N+2, stdout );
	fflush( stdout );
	continue;
      }

      rdft( N, -1, buffer1, bitshuffle, trigland );      
      overlapadd( buffer1, N, Wsyn, output, Nw, on );
      shiftout( output, Nw, I, on );
    }

    exit(0);
}


int usage(woof)
int woof;
{
    fprintf(stderr,
	"burrow:  FFT based spectral ellision\n"
	"%% burrow   [flags] > floatsams\n"
	"	N:	fft length [2^n] (1024)\n"
	"	R:	sampling rate (44100)\n"
	"	M:	window size in samples (N)\n"
	"	D:	decimation factor in samples (M/4)\n"
	"	I:	interpolation factor in samples (D)\n"
	"	f:	source soundfile\n"
	"	F:	filter soundfile\n"
	"	m:	filtering multiplier (.001)\n"
	"	t:	filtering threshold (.001)\n"
	"	r:	adjust threshold by rms of source sound\n"
	"	w:	time smoothing window size [spectral frames] (1)\n"
	"	a:	analysis data output\n"
	"	h:	this cosmic place\n");
    exit(woof);
}