#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>

static unsigned int FREQUENCY = 8192;

static double (*waveForm)( double angle );

enum {
    BYTES_PER_SAMPLE = 2
};

enum {
    MAX_SAMPLE = (~(~0 << (8*BYTES_PER_SAMPLE-1)))
};

static int
readWidth( unsigned int* width, FILE* fp )
{
    int headerNumbersRead = 0;
    int inNumber = 0;
    int number = 0;
    int inHeader = 1;
    int ch;

    if ( fp == 0 ) {
	return __LINE__;
    }

    while ( inHeader && ( ch = fgetc( fp ) ) != EOF ) {
	if ( isdigit( ch ) ) {
	    if ( !inNumber ) {
		number = 0;
		inNumber = 1;
	    }

	    number *= 10;
	    number += ( ch - '0' );
	} else {
	    if ( inNumber ) {
		switch ( ++headerNumbersRead ) {
		case 1:
		    *width = number;
		    inHeader = 0;
		    break;
		}
		inNumber = 0;
	    }

	    if ( ch == '#' ) {
		while ( ( ch = fgetc( fp ) ) != EOF && ch != '\n' ) {
		    /* do nothing */;
		}
	    } else if ( !isspace( ch ) ) {
		fprintf( stderr, "Illegal header character <%c>\n", ch );
	    }
	}
    }

    if ( ch != '\n' ) {
	while ( ( ch = fgetc( fp ) ) != EOF && ch != '\n' ) {
	    /* do nothing */;
	}
    }

    if ( *width == 0 ) {
	fprintf( stderr, "Invalid width\n" );
	fclose( fp );
	return __LINE__;
    }

    return 0;
}

static double intensityMap[ 256 ];

static int
initIntensityMap( double gamma )
{
    unsigned int ii;

    for ( ii=0; ii < 256; ++ii ) {
	intensityMap[ ii ] = pow( (double)ii / 255.0, 1.0 / gamma );
    }

    return 0;
}

static double
remapIntensity( unsigned int index )
{
    return intensityMap[ index ];
}

static double
squareWave( double angle )
{
    double waveNumber;
    double proportion;

    proportion = modf( angle / ( 2.0 * M_PI ), &waveNumber );

    if ( proportion < 0.5 ) {
	return 1.0;
    } else {
	return -1.0;
    }
}

static double
triangleWave( double angle )
{
    double waveNumber;
    double proportion;

    proportion = modf( angle / ( 2.0 * M_PI ), &waveNumber );

    if ( proportion < 0.5 ) {
	return 1.0 - 4.0 * proportion;
    } else {
	return 4.0 * proportion - 3.0;
    }
}

static double
sawtoothWave( double angle )
{
    double waveNumber;
    double proportion;

    proportion = modf( angle / ( 2.0 * M_PI ), &waveNumber );

    return 2.0 * proportion - 1.0;
}

static double
toothsawWave( double angle )
{
    double waveNumber;
    double proportion;

    proportion = modf( angle / ( 2.0 * M_PI ), &waveNumber );

    return 1.0 - 2.0 * proportion;
}

static double (*waveForm)( double angle );

static double
getSample( double intensity, unsigned int sampleNumber, double hz )
{
    return intensity
	* (*waveForm)( sampleNumber * hz * 2.0 * M_PI / (double)FREQUENCY );
}

static void
emitIntegerSample( double samp )
{
    int val = (int)( MAX_SAMPLE * samp );

    if ( val < -MAX_SAMPLE ) {
	val = -MAX_SAMPLE;
    } else if ( val > MAX_SAMPLE ) {
	val = MAX_SAMPLE;
    }

    fputc( ( val >> 8 ) & 0xFF, stdout );
    fputc( ( val >> 0 ) & 0xFF, stdout );
}

static void
emitFloatingSample( double samp )
{
    fwrite( (const char*)&samp, sizeof(samp), 1, stdout );
}

int
main( int argc, const char* argv[] )
{
    unsigned int help = 0;
    double linesPerSecond = 64.0;
    double gamma = 0.44;
    double minHz = 128.0;
    double maxHz = 4096.0;
    double logMinHz;
    double logMaxHz;
    double* hzTable;
    double gain = 1.00;
    unsigned int ii;
    unsigned int yy = 0;
    unsigned int width;
    unsigned char* buf[2];
    unsigned int primaryBuf = 0;
    void (*emitSample)( double ) = emitIntegerSample;

    for ( ii=1; ii < argc; ++ii ) {
	if ( !strcmp( argv[ ii ], "-lps" ) ) {
	    linesPerSecond = strtod( argv[ ++ii ], 0 );
	} else if ( !strcmp( argv[ ii ], "-gamma" ) ) {
	    gamma = strtod( argv[ ++ii ], 0 );
	} else if ( !strcmp( argv[ ii ], "-minHz" ) ) {
	    minHz = strtod( argv[ ++ii ], 0 );
	} else if ( !strcmp( argv[ ii ], "-maxHz" ) ) {
	    maxHz = strtod( argv[ ++ii ], 0 );
	} else if ( !strcmp( argv[ ii ], "-gain" ) ) {
	    gain = strtod( argv[ ++ii ], 0 );
	} else if ( !strcmp( argv[ ii ], "-freq" ) ) {
	    FREQUENCY = strtol( argv[ ++ii ], 0, 0 );
	} else if ( !strcmp( argv[ ii ], "-wave" ) ) {
	    ++ii;
	    if ( !strcmp( argv[ ii ], "sine" ) ) {
		waveForm = sin;
	    } else if ( !strcmp( argv[ ii ], "square" ) ) {
		waveForm = squareWave;
	    } else if ( !strcmp( argv[ ii ], "triangle" ) ) {
		waveForm = triangleWave;
	    } else if ( !strcmp( argv[ ii ], "sawtooth" ) ) {
		waveForm = sawtoothWave;
	    } else if ( !strcmp( argv[ ii ], "toothsaw" ) ) {
		waveForm = toothsawWave;
	    } else {
		fprintf( stderr, "Unknown waveform: %s\n", argv[ii] );
		help = 1;
	    }
	} else if ( !strcmp( argv[ ii ], "-float" ) ) {
	    emitSample = emitFloatingSample;
	} else {
	    help = 1;
	}
    }

    if ( FREQUENCY == 0 || gain < 0.0000001 || linesPerSecond < 0.0000001 ) {
	fprintf( stderr, "Bad parameters\n" );
	help = 1;
    }

    if ( help ) {
	fprintf( stderr, "Usage: %s\n", argv[0] );
	fprintf( stderr, "\t[-lps linesPerSecond]\n" );
	fprintf( stderr, "\t[-gamma gammaCorrection]\n" );
	fprintf( stderr, "\t[-minHz minimumFrequency]\n" );
	fprintf( stderr, "\t[-maxHz maximumFrequency]\n" );
	fprintf( stderr, "\t[-gain gain]\n" );
	fprintf( stderr, "\t[-freq outputSamplesPerSecond]\n" );
	fprintf( stderr,
		"\t[-wave sine | triangle | square | sawtooth | toothsaw ]\n"
	    );
	fprintf( stderr, "\t[-float]\n" );
	fprintf( stderr, "\t[-help]\n" );
	return __LINE__;
    }

    if ( minHz < 1.0 ) {
	minHz = 1.0;
    }

    if ( maxHz < 1.0 ) {
	maxHz = 1.0;
    }

    if ( readWidth( &width, stdin ) != 0 ) {
	fprintf( stderr, "Couldn't get width\n" );
	return __LINE__;
    }

    buf[ 0 ] = (unsigned char*)malloc( width );
    buf[ 1 ] = (unsigned char*)malloc( width );

    if ( buf[ 0 ] == 0 || buf[ 1 ] == 0 ) {
	fprintf( stderr, "Couldn't get memory for buffers\n" );
	return __LINE__;
    }

    initIntensityMap( gamma );

    hzTable = (double*)malloc( width * sizeof( double ) );

    if ( hzTable == 0 ) {
	fprintf( stderr, "Couldn't get hzTable memory\n" );
	return __LINE__;
    }

    logMinHz = log( minHz );
    logMaxHz = log( maxHz );

    for ( ii=0; ii < width; ++ii ) {
	double horz = (double) ii / (double)width;
	double logHz = logMinHz + horz * ( logMaxHz - logMinHz );

	hzTable[ ii ] = exp( logHz );
    }

    if ( waveForm == 0 ) {
	waveForm = sin;
    }

	/* adjust gain to avoid clipping */
    gain /= (double) width;

    fread( buf[ 0 ], 1, width, stdin );

    while ( !feof( stdin ) ) {
	static int previousRow = 0;
	double vert = ( (double) yy * linesPerSecond ) / (double)FREQUENCY;
	unsigned int kk;
	double samp = 0.0;
	int val;
	double realRow;
	double proportion;
	int row;
	int otherRow;
	unsigned char* aPtr;
	unsigned char* bPtr;

	proportion = modf( vert, &realRow );
	otherRow = ( (int)realRow ) % 2;
	row = ( otherRow + 1 ) % 2;

	if ( row != previousRow ) {
	    previousRow = row;
	    if ( fread( buf[ row ], 1, width, stdin ) != width ) {
		break;
	    }
	}

	aPtr = buf[ row ];
	bPtr = buf[ otherRow ];

	for ( kk=0; kk < width; ++kk ) {
	    double partA = remapIntensity( *aPtr++ );
	    double partB = remapIntensity( *bPtr++ );
	    double intensity;

	    intensity = proportion * partA + ( 1.0 - proportion ) * partB;

	    samp += getSample( intensity, yy, hzTable[kk] );
	}

	(*emitSample)( samp * gain );

	++yy;
    }

    return 0;
}