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

#define Max_Table_Size 131073

typedef struct _pickme {

  int		bin;
  float		value;

} pickme;

int sortIncreasing( a, b )
pickme *a, *b;
{

  if ( a->value > b->value )
    return 1;

  if ( a->value < b->value )
    return -1;

  return 0;
}     

int sortDecreasing( a, b )
pickme *a, *b;
{

  if ( a->value < b->value )
    return 1;

  if ( a->value > b->value )
    return -1;

  return 0;
}     


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,
                real, imag,
                amp, phase,
		index,
		*bitshuffle,
		sortDecreasing(),
		sortIncreasing();

    float       ctime = 0.,
                tinc,
                fpos = 0.,
                scalar = 1.,
                a1,a2,
                b1,b2,
                starttime = -1.,
                endtime = -1.,
		increasing = 1,
		*trigland,
                *Hwin,
                *Wanal,
                *Wsyn,
                *input1,
                *input2,
                *buffer1,
                *buffer2,
                *channel1,
                *channel2,
                *output,
                prand();

    double      rmult = 1.;

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

    FILE        *fp1,
                *fp2;

    Tabl        *ftabl1;

    pickme	*picks,
		*mirror;

    MrSoundStruct       soundStruct1,
			soundStruct2;


    synt = 0.0;

    while( (ch= crack( argc, argv, "R|N|M|D|I|f|F|t|T|s|e|d|S|X|h", 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':   rmult = atof(arg_option);
                        break;
	case 't':   strcpy(table1, arg_option);
                        break;
	case 's':   starttime = atof(arg_option);
                        break;
	case 'e':   endtime = atof(arg_option);
                        break;
	case 'S':   srrand(atoi(arg_option));
                        break;
	case 'X':   scalar = atof(arg_option);
                        break;
	case 'd':   increasing = 0;
                        break;
	case 'h':   usage(1);
        }
      }


    if ( starttime < 0. || endtime < 0. ) {
      fprintf(stderr,"Please specify start and end times!\n");
      usage(1);
    }

    if (Nw == 0)
        Nw = N;

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

    if (I == 0)
      I = D;

    myPI = 4.*atan(1.);
    myTWOPI = 8.*atan(1.);
    N2 = N>>1;
    Nw2 = Nw>>1;
    tinc = (float) D / (float) R;

    if ( strcmp(file1, file2) == 0 ) {
      fprintf(stderr,"lame and boring file selection\n");
      usage(1);
    }

    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) );
    Wanal = (float *) space( Nw, sizeof(float) );
    Wsyn = (float *) space( Nw, sizeof(float) );
    Hwin = (float *) space( Nw, sizeof(float) );
    input1 = (float *) space( Nw, sizeof(float) );
    input2 = (float *) space( Nw, sizeof(float) );
    output = (float *) space( Nw, sizeof(float) );

    picks = (pickme *) space( N2+1, sizeof(pickme) );
    mirror = (pickme *) space( N2+1, sizeof(pickme) );

    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) (N>>1)), sizeof(int) );
  
    /* initialize FFT needs */

    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);
    }



/* read in traversal function table */

    ftabl1 = loadtabl( table1 );

    while ( !eof ) {

      in += D;
      on += I;

      eof = floatin( input1, Nw, D, fp1 );
      if (eof)
        floatin( input2, Nw, D, fp2 );
      else
        eof = floatin( input2, Nw, D, fp2 );
      
      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 );

      if (ctime >= starttime) {

        if (ctime >= endtime) {

          for ( i = 0 ; i <= N2 ; i++ ) {
            imag = phase = ( real = amp = i<<1 ) + 1 ;

/*
            a1 = ( i == N2 ? *(buffer1+1) : *(buffer1+real) ) ;
            b1 = ( i == 0 || i == N2 ? 0. : *(buffer1+imag) ) ;
*/

	    a2 = ( i == N2 ? *(buffer2+1) : *(buffer2+real) ) ;
            b2 = ( i == 0 || i == N2 ? 0. : *(buffer2+imag) ) ;

/*
            *(channel1+phase) = -atan2( b1, a1 );
*/

            *(channel2+amp) = hypot( a2, b2 );
            *(channel2+phase) = -atan2( b2, a2 );
          }

          for ( i=0; i <= N2; i++ ) {
            phase = (amp = i<<1) + 1;

              *(channel1+amp) = *(channel2+amp);

/*	    if ( *(channel1+phase) == 0. ) */
              *(channel1+phase) = *(channel2+phase);
          }

          for ( i = 0 ; i <= N2 ; i++ ) {
            imag = phase = ( real = amp = i<<1 ) + 1 ;
            *(buffer1+real) = *(channel1+amp) * cos( *(channel1+phase) );
            if ( i != N2 ) 
              *(buffer1+imag) = -*(channel1+amp) * sin( *(channel1+phase) );
          }
        }

        else {

          fpos = (ctime - starttime) / (endtime - starttime) * 
                ((float) ftabl1->elms);

          for ( i = 0 ; i <= N2 ; i++ ) {
            imag = phase = ( real = amp = i<<1 ) + 1 ;
            a1 = ( i == N2 ? *(buffer1+1) : *(buffer1+real) ) ;
            b1 = ( i == 0 || i == N2 ? 0. : *(buffer1+imag) ) ;
            a2 = ( i == N2 ? *(buffer2+1) : *(buffer2+real) ) ;
            b2 = ( i == 0 || i == N2 ? 0. : *(buffer2+imag) ) ;

            *(channel1+amp) = hypot( a1, b1 );
            *(channel1+phase) = -atan2( b1, a1 );
            *(channel2+amp) = hypot( a2, b2 );
            *(channel2+phase) = -atan2( b2, a2 );

	    /* find amplitude differences between home and visitors */

	    (picks+i)->value = fabs( *(channel1+amp) - *(channel2+amp) );
	    (picks+i)->bin = i;
          }
	  

	  if (increasing)
	    qsort( (void *) picks, N2+1, sizeof(pickme), sortIncreasing );

	  else
	    qsort( (void *) picks, N2+1, sizeof(pickme), sortDecreasing );


	  /* now we create an effective mirror of the sorted distribution.
	     we will assure that the initial transition will be made from
	     small spectral differences (when the sort behavior is increasing)
	     and the ending transition will also be made from small spectral
             differences */ 

	  for ( i=0; i <= N2; i += 2 ) {
	    (mirror+(i/2))->bin = (picks+i)->bin;
	    (mirror+(i/2))->value = (picks+i)->value;
	  }

	  for ( i=1; i <= N2; i += 2 ) {
	    (mirror+(N2-(i/2)))->bin = (picks+i)->bin;
	    (mirror+(N2-(i/2)))->value = (picks+i)->value;
	  }

	  /*

          for ( i=0; i <= N2; i++ )
	    fprintf(stderr, "%d		%d	%f\n", i,
				 (mirror+i)->bin, (mirror+i)->value);

	  fprintf(stderr, "\n");

	  */

	  index = (*(ftabl1->data+(int)(fpos + .5))) * (N2+1);

	  /* choose the bins that are least different first */

          for ( i=0; i <= index; i++ ) {

	    int		amp = ((mirror+i)->bin)<<1,
			phase = (((mirror+i)->bin)<<1) + 1;	

	    *(channel1+amp) = *(channel2+amp);
	    *(channel1+phase) = *(channel2+phase);
          }


          for ( i = 0 ; i <= N2 ; i++ ) {
            imag = phase = ( real = amp = i<<1 ) + 1 ;
            *(buffer1+real) = *(channel1+amp) * cos( *(channel1+phase) );

            if ( i != N2 ) 
              *(buffer1+imag) = -*(channel1+amp) * sin( *(channel1+phase) );
          }
        }
      }

      else {

	for ( i = 0 ; i <= N2 ; i++ ) {
	  imag = phase = ( real = amp = i<<1 ) + 1 ;
	  a1 = ( i == N2 ? *(buffer1+1) : *(buffer1+real) ) ;
	  b1 = ( i == 0 || i == N2 ? 0. : *(buffer1+imag) ) ;
	  a2 = ( i == N2 ? *(buffer2+1) : *(buffer2+real) ) ;
	  b2 = ( i == 0 || i == N2 ? 0. : *(buffer2+imag) ) ;
	  *(channel1+amp) = hypot( a1, b1 );
	  *(channel1+phase) = -atan2( b1, a1 );
	  *(channel2+amp) = hypot( a2, b2 );
	  *(channel2+phase) = -atan2( b2, a2 );
	}
	        
        for ( i = 0 ; i <= N2 ; i++ ) {
          imag = phase = ( real = amp = i<<1 ) + 1 ;
          *(buffer1+real) = *(channel1+amp) * cos( *(channel1+phase) );
          if ( i != N2 ) 
            *(buffer1+imag) = -*(channel1+amp) * sin( *(channel1+phase) );
        }
      }

      ctime += tinc;

      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,
        "morphine:  adaptive spectral morphing\n"
        "%% morphine   [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/8)\n"
        "       I:      interpolation factor in samples (D)\n"
        "       f:      departure soundfile (unspecified)\n"
        "       F:      arrival soundfile (unspecified)\n"
        "       t:      transition mapping function (unspecified)\n"
        "       s:      begin time of interpolation (unspecified)\n"
        "       e:      end time of interpolation (unspecified)\n"
        "       d:      use decreasing sorting behavior\n"
        "       h:      this bouncy place\n");
    exit(woof);
}