/*
* wav8.c  - runs magnetic simulator on raw samples.
* usage: ..\tcc -run wav8.c 
* usage: gcc wav8.c && ./a.out && aplay test.wav
*
* usage: wav8 inputRawfilename outputRawFilename   
* 
* 
*
* idea use fseek to position and update header.
* 1) write header first - final length is unknown for now.
* 2) append data 
* 3) seek END to seek to final length, use ftell to get length
* 4) seek to "data" and update length.


*
Magnetict sticktion only updates if teh input has changed significantly, but is not quantised.

My idea was to have schmitt triggers that fed back and adjusted the bias.

*/

#ifndef MAKE_WAV_H
#define MAKE_WAV_H

void write_wav(char * filename, unsigned long num_samples, short int * data, int s_rate);
    /* open a file named filename, write signed 16-bit values as a
        monoaural WAV file at the specified sampling rate
        and close the file
    */

#endif

/*
http://karplus4arduino.wordpress.com/2011/10/08/making-wav-files-from-c-programs/
karplus4arduino
http://gasstationwithoutpumps.wordpress.com/2011/10/08/making-wav-files-from-c-programs/

2011 October 8

Making WAV files from C programs

Filed under: Digital music — gasstationwithoutpumps @ 21:03
Tags: software, WAV

I’m going to try embedding some short pieces of code in this post, and later work out a better way of distributing larger blocks of code.

The code I want to share today consists of a small C file and the associated header for creating WAV files for producing sounds on a laptop or desktop computer.
This is not a full-featured WAV file module: the format has a lot of options and I didn’t want the complexity of dealing with all of them.
I just wanted to output monophonic sound with a reasonable sampling rate.

*/

/* make_wav.c
 * Creates a WAV file from an array of ints.
 * Output is monophonic, signed 16-bit samples
 * copyright
 * Fri Jun 18 16:36:23 PDT 2010 Kevin Karplus
 * Creative Commons license Attribution-NonCommercial
 *  http://creativecommons.org/licenses/by-nc/3.0/
 */

#include <stdio.h>
#include <assert.h>
#include <math.h>

#define S_RATE  (8000)
#define BUF_SIZE (S_RATE/2) /* 2 second buffer */
short  buffer[BUF_SIZE];

#define MS_NUM_GRAINS	3
#define MS_MAX	100	/* 100 * 100 = 10000 */
#define MS_MIN	-100

short ms_last_positionA[ MS_NUM_GRAINS ] ;
short ms_delta[ MS_NUM_GRAINS ] ;

short ms_state[ MS_NUM_GRAINS ] ;	/* schmitt state 1 or -1 */
short ms_threshold[ MS_NUM_GRAINS ] ;	/* 0 to 1 */	
short ms_output[ MS_NUM_GRAINS ] ;	/* 0 to 1 */
short ms_feedforward[ MS_NUM_GRAINS ] ; /* 0 to 1 */
short ms_feedback[ MS_NUM_GRAINS ] ;	/* 0 to 1 */



#ifndef M_PI
#define M_PI 3.141592654
#endif

int fillBuffer( int t1, int t2, int level  ){
    int i;

    int x1,y1;
    int x2,y2;
    int n;

    float t;
    float amplitude = 32000;
    float freq_Hz = 440;
    float phase=0;

    float freq_radians_per_sample = freq_Hz*2*M_PI/S_RATE;

    x1 = level ;
    y1 = 0;

    x2 = level ;
    y2 = 0;

    /* fill buffer with a sine wave */
    for (i=0; i<BUF_SIZE; i++)
    {
        phase += freq_radians_per_sample;
        //buffer[i] = (int)(amplitude * sin(phase));

        y1 = y1 + x1 / t1;
        x1 = x1 - y1 / t1;
        y2 = y2 + x2 / t2;
        x2 = x2 - y2 / t2;
        buffer[ i ] = x1+x2;
    }
}


void write_little_endian(unsigned int word, int num_bytes, FILE *wav_file)
{
    unsigned buf;
    while(num_bytes>0)
    {   buf = word & 0xff;
        fwrite(&buf, 1,1, wav_file);
        num_bytes--;
    word >>= 8;
    }
}

/* information about the WAV file format from
http://ccrma.stanford.edu/courses/422/projects/WaveFormat/
*/

void write_wav(char * filename, unsigned long num_samples, short int * data, int s_rate)
{
    FILE* wav_file;
    unsigned int sample_rate;
    unsigned int num_channels;
    unsigned int bytes_per_sample;
    unsigned int byte_rate;
    unsigned long i;    /* counter for samples */

    num_channels = 1;   /* monoaural */
    bytes_per_sample = 2;

    if (s_rate<=0) sample_rate = S_RATE;
    else sample_rate = (unsigned int) s_rate;

    byte_rate = sample_rate*num_channels*bytes_per_sample;

    wav_file = fopen(filename, "wb");
    assert(wav_file);   /* make sure it opened */

    /* write RIFF header */
    fwrite("RIFF", 1, 4, wav_file);
    write_little_endian(36 + bytes_per_sample* num_samples*num_channels, 4, wav_file);
    fwrite("WAVE", 1, 4, wav_file);

    /* write fmt  subchunk */
    fwrite("fmt ", 1, 4, wav_file);
    write_little_endian(16, 4, wav_file);   /* SubChunk1Size is 16 */
    write_little_endian(1, 2, wav_file);    /* PCM is format 1 */
    write_little_endian(num_channels, 2, wav_file);
    write_little_endian(sample_rate, 4, wav_file);
    write_little_endian(byte_rate, 4, wav_file);
    write_little_endian(num_channels*bytes_per_sample, 2, wav_file);  /* block align */
    write_little_endian(8*bytes_per_sample, 2, wav_file);  /* bits/sample */

    /* write data subchunk */
    fwrite("data", 1, 4, wav_file);
    i = bytes_per_sample* num_samples*num_channels;
    i = 0x5555;
    
    printf("%ld\n",i);
    	
    write_little_endian(bytes_per_sample* num_samples*num_channels, 4, wav_file);
    for (i=0; i< num_samples; i++)
    {   
	write_little_endian((unsigned int)(data[i]),bytes_per_sample, wav_file);
    }
    fclose(wav_file);
}


void write_wav_header(char * filename, unsigned long num_samples, short int * data, int s_rate)
{
    FILE* wav_file;
    unsigned int sample_rate;
    unsigned int num_channels;
    unsigned int bytes_per_sample;
    unsigned int byte_rate;
    unsigned long i;    /* counter for samples */

    num_channels = 1;   /* monoaural */
    bytes_per_sample = 2;

    if (s_rate<=0) sample_rate = S_RATE;
    else sample_rate = (unsigned int) s_rate;

    byte_rate = sample_rate*num_channels*bytes_per_sample;

    wav_file = fopen(filename, "wb+");
    assert(wav_file);   /* make sure it opened */
    fseek( wav_file,0,SEEK_SET);

    /* write RIFF header */
    fwrite("RIFF", 1, 4, wav_file);
    write_little_endian(36 + bytes_per_sample* num_samples*num_channels, 4, wav_file);
    fwrite("WAVE", 1, 4, wav_file);

    /* write fmt  subchunk */
    fwrite("fmt ", 1, 4, wav_file);
    write_little_endian(16, 4, wav_file);   /* SubChunk1Size is 16 */
    write_little_endian(1, 2, wav_file);    /* PCM is format 1 */
    write_little_endian(num_channels, 2, wav_file);
    write_little_endian(sample_rate, 4, wav_file);
    write_little_endian(byte_rate, 4, wav_file);
    write_little_endian(num_channels*bytes_per_sample, 2, wav_file);  /* block align */
    write_little_endian(8*bytes_per_sample, 2, wav_file);  /* bits/sample */

    /* write data subchunk */
    fwrite("data", 1, 4, wav_file);
    
    i = bytes_per_sample* num_samples*num_channels;
    i = 0x55555555;
    i = 0x10;
        
    //printf("write_wav_header: %ld\n",i);

    write_little_endian( i , 4, wav_file);
/*
    for (i=0; i< num_samples; i++)
    {   
	write_little_endian((unsigned int)(data[i]),bytes_per_sample, wav_file);
    }
*/
    fclose(wav_file);
}

void write_wav_length(char * filename, unsigned long num_samples, short int * data, int s_rate)
{
    FILE* wav_file;
    unsigned long i;    /* counter for samples */
    
    wav_file = fopen( filename, "rb+");
    assert(wav_file);   /* make sure it opened */

    /* find the end of the file in bytes */	
    fseek( wav_file,0,SEEK_END);
    i = ftell(wav_file);

    /* now update data length at the start of the file */
    fseek( wav_file,36,SEEK_SET);
    /* write data subchunk */
    fwrite("data", 1, 4, wav_file);
    //printf("write_wav_header: %ld 0x%lX\n",i,i);
    write_little_endian( i , 4, wav_file);
    
    fseek( wav_file,0,SEEK_END);
    fclose(wav_file);
}


void write_wav_data(char * filename, unsigned long num_samples, short int * data, int s_rate)
{


    FILE* wav_file;
    unsigned int sample_rate;
    unsigned int num_channels;
    unsigned int bytes_per_sample;
    unsigned int byte_rate;
    unsigned long i;    /* counter for samples */

    num_channels = 1;   /* monoaural */
    bytes_per_sample = 2;

    if (s_rate<=0) sample_rate = S_RATE;
    else sample_rate = (unsigned int) s_rate;

    byte_rate = sample_rate*num_channels*bytes_per_sample;

    /* apend data samples onto the end of the file, but it does not update length */	    	
    wav_file = fopen(filename, "ab+");
    assert(wav_file);   /* make sure it opened */
    /* if we use "ab+" we append, so we do not need to seek */	
//    fseek( wav_file,36,SEEK_SET);


    for (i=0; i< num_samples; i++)
    {   
	write_little_endian((unsigned int)(data[i]),bytes_per_sample, wav_file);
    }
    fclose(wav_file);
}



void copy_wav_data(char * filename,char * filename_in )
{
	short buf[100];
    FILE* wav_file;
    FILE* wav_file_in;
    size_t got;
	short sample;
    short * samplesP;

	
	
    /* apend data samples onto the end of the file, but it does not update length */	    	
    wav_file    = fopen(filename,    "ab+");
    wav_file_in = fopen(filename_in, "rb" );
	if ( ! wav_file    ){ printf( "error: output file ");    }
	if ( ! wav_file_in ){ printf( "error: Input file ");     }

	// printf("%p %p",wav_file,wav_file_in );

	if ( wav_file_in ){	
		while(1){
			/* 
			local call recording is 44100, mono. 
			44100/8000 = 5.5125  
			*/		
			got = fread( &buf[0], sizeof( short ),1, wav_file_in );
			//got = fread( &buf[0], sizeof( short ),11, wav_file_in );
			if( feof(wav_file_in) )
			{ 
				break ;
			}
			/*  *0,1,2,3,4,5, *6,7,8,9,10*/
			if ( wav_file ){ fwrite(&buf[0], sizeof( short ),1, wav_file); }
			//if ( wav_file ){ fwrite(&buf[6], sizeof( short ),1, wav_file); }
			//printf( "%d", got );
		};	
		if ( wav_file    ){ fclose( wav_file );    }
		if ( wav_file_in ){ fclose( wav_file_in ); }
	}
}



/* 
* CD quality is 44100 samples per second. 
* We want 8000 samples per second
* 44100/8000 = 5.5125
* 441  /80	 = 5.5125
* 
* The samples seem to be mono, so I suggested taking every 11 samples. 
* This would be okay for 44000, but we sample at 44100. 
* You need a leap sample evey 80 CD samples.
* I used Audacity to resample
*
* convert mono, 44100 sample/second to 8000 samples/second 
*
* Every 11 samples take sample 0 and sample 6 and 
*/

void copy_wav_data_cd_8k_doesNotWork(char * filename,char * filename_in )
{
	short buf[100];
    FILE* wav_file;
    FILE* wav_file_in;
    size_t got;

	short sample40;
	
	printf("\n===============copy_wav_data_cd_8k()==================================================================\n");
	
    /* apend data samples onto the end of the file, but it does not update length */	    	
    wav_file    = fopen(filename,    "wb+");
    wav_file_in = fopen(filename_in, "rb" );
	if ( ! wav_file    ){ printf( "error: output file ");    }
	if ( ! wav_file_in ){ printf( "error: Input file ");     }

	// printf("%p %p",wav_file,wav_file_in );
	sample40 = 1;
	
	if ( wav_file_in ){	
		while(1){
			/* 
			local call recording is 44100, mono. 
			44100/8000 = 5.5125  
			*/		
			//got = fread( &buf[0], sizeof( short ),1, wav_file_in );
			sample40 += 11;
			if ( sample40 > 80){
			  /* every 40 samples throw one away this is for the 0.1 in 44.1ksamples */	
			  got = fread( &buf[0], sizeof( short ),1, wav_file_in );
              sample40 = sample40 -80 ;		
	          printf("\n=== %d ==\n", sample40 );
			  
			}
			got = fread( &buf[0], sizeof( short ),11, wav_file_in );
			if( feof(wav_file_in) )
			{ 
				break ;
			}
			/*  *0,1,2,3,4,5, *6,7,8,9,10 */
			if ( wav_file ){ fwrite(&buf[0], sizeof( short ),1, wav_file); }
			if ( wav_file ){ fwrite(&buf[6], sizeof( short ),1, wav_file); }
			//printf( "%d", got );
		};	
		if ( wav_file    ){ fclose( wav_file );    }
		if ( wav_file_in ){ fclose( wav_file_in ); }
	}
}



/*
*
* The 44100 samples need to be subsampled down to 8000.
* Use the method of adding absolute time
*
* for each CD_sample add 1/41000 to CD time
* when CD_time > 125us output a sample and subtract 125us
*
*/

void copy_wav_data_cd_8k_works(char * filename,char * filename_in )
{
	short buf[100];
    FILE* wav_file;
    FILE* wav_file_in;
    size_t got;
	
	float cd_time;
	
	/* 
	for each CD_sample add 1/41000 to CD time

	short sample40;
	
	printf("\n===============copy_wav_data_cd_8k()==================================================================\n");
	
    /* apend data samples onto the end of the file, but it does not update length */	    	
    wav_file    = fopen(filename,    "wb+");
    wav_file_in = fopen(filename_in, "rb" );
	if ( ! wav_file    ){ printf( "error: output file ");    }
	if ( ! wav_file_in ){ printf( "error: Input file ");     }

	// printf("%p %p",wav_file,wav_file_in );
	cd_time=0.0;
	if ( wav_file_in ){	
		while(1){
			/* 
			local call recording is 44100, mono. 
			44100/8000 = 5.5125  
			*/		
			//got = fread( &buf[0], sizeof( short ),1, wav_file_in );
			got = fread( &buf[0], sizeof( short ),1, wav_file_in );
			cd_time +=1/44100.0;
			
			if( feof(wav_file_in) )
			{ 
				break ;
			}
			if ( cd_time > 1/8000.0){	
			  if ( wav_file ){ fwrite(&buf[0], sizeof( short ),1, wav_file); }
			  //printf( "%d", got );
			  cd_time -= 1/8000.0;
			}  
		};	
		if ( wav_file    ){ fclose( wav_file );    }
		if ( wav_file_in ){ fclose( wav_file_in ); }
	}
}




/*
*
* The 44100 samples need to be subsampled down to 8000.
* Use the method of adding absolute time
*
* for each CD_sample add 1/41000 to CD time
* when CD_time > 125us output a sample and subtract 125us
*
*/

/*
*
*  
* 
* only updae if input has changed significantly!
* 
* #define MS_NUM_GRAINS	10
* short int ms_last_positionA[ NUM_GRAINS ]
* short int ms_delta[ NUM_GRAINS ]
*
*
*   if (( input-last_position )> ms_delta){
*     lasp_positionA = input
*   }
*
*
*
*/
void magnetic_simulate_sticktion(char * filename,char * filename_in )
{
	short buf[100];
   	FILE* wav_file;
   	FILE* wav_file_in;
    	size_t got;
	
	float cd_time;
	
	/* 
	for each CD_sample add 1/41000 to CD time
	*/
	short sample40;
	short input;
	short middle;
	short output;
	int count;

	printf("\n===============copy_wav_data_cd_8k()==================================================================\n");
	
	for( count = 0 ; count < MS_NUM_GRAINS ; count++){
          ms_delta[ 0 ] = 0;
          ms_last_positionA[ count  ] = 0; 
        }

	ms_delta[ 0 ] = 5000;
	ms_delta[ 1 ] = 3000;
	ms_delta[ 2 ] = 5000;
	ms_delta[ 3 ] = 7000;
	ms_delta[ 4 ] = 1100;

    /* apend data samples onto the end of the file, but it does not update length */	    	
    wav_file    = fopen(filename,    "wb+");
    wav_file_in = fopen(filename_in, "rb" );
	if ( ! wav_file    ){ printf( "error: output file ");    }
	if ( ! wav_file_in ){ printf( "error: Input file ");     }

	// printf("%p %p",wav_file,wav_file_in );
	cd_time=0.0;
	if ( wav_file_in ){	
		while(1){
			//got = fread( &buf[0], sizeof( short ),1, wav_file_in );
			got = fread( &buf[0], sizeof( short ),1, wav_file_in );

			if( feof(wav_file_in) )
			{ 
				break ;
			}
			input = buf[0];		
			for( count = 0 ; count < MS_NUM_GRAINS ; count++){

			  /* stiction simulation */
			  if ( abs( input - ms_last_positionA[ count  ] ) > ms_delta[ count ] ){
                            ms_last_positionA[ count  ] = input ;  
                          }
			  input = ms_last_positionA[ count  ];
                        }
			input = (input-ms_last_positionA[ count  ]);

	                buf[0] = input;
			if ( wav_file ){ fwrite(&buf[0], sizeof( short ),1, wav_file); }
			  
		};	
		if ( wav_file    ){ fclose( wav_file );    }
		if ( wav_file_in ){ fclose( wav_file_in ); }
	}
}









/*
*
*  
* 
* only updae if input has changed significantly!
* 
* #define MS_NUM_GRAINS	10
* short int ms_last_positionA[ NUM_GRAINS ]
* short int ms_delta[ NUM_GRAINS ]
*
*
*   if (( input-last_position )> ms_delta){
*     lasp_positionA = input
*   }
*
*
*
*/
void magnetic_simulate_schmitt(char * filename,char * filename_in )
{
	short buf[100];
   	FILE* wav_file;
   	FILE* wav_file_in;
    	size_t got;
	
	float cd_time;
	
	/* 
	for each CD_sample add 1/41000 to CD time
	*/
	short sample40;
	short input;
	short output;
	short middle;
	short output;
	int count;

	printf("\n===============copy_wav_data_cd_8k()==================================================================\n");
	
	for( count = 0 ; count < MS_NUM_GRAINS ; count++){
          ms_delta[ count ] = 0;
          ms_last_positionA[ count  ] = 0; 

	 /* iput += 16bit signed, so MS_MAX = 100 for now */
  	  ms_state[ count ] 	  = MS_MAX ;	/* schmitt state 1 or -1 */
  	  ms_threshold[ count ]   = MS_MAX ;	/* 0 to 1 */	
  	  ms_output[ count ]      = MS_MAX ;	/* 0 to 1 */
  	  ms_feedforward[ count ] = MS_MAX ;    /* 0 to 1 */
  	  ms_feedback[ count ] 	  = MS_MAX ;	/* 0 to 1 */
        }


    /* apend data samples onto the end of the file, but it does not update length */	    	
    wav_file    = fopen(filename,    "wb+");
    wav_file_in = fopen(filename_in, "rb" );
	if ( ! wav_file    ){ printf( "error: output file ");    }
	if ( ! wav_file_in ){ printf( "error: Input file ");     }

	// printf("%p %p",wav_file,wav_file_in );
	cd_time=0.0;
	if ( wav_file_in ){	
		while(1){
			//got = fread( &buf[0], sizeof( short ),1, wav_file_in );
			got = fread( &buf[0], sizeof( short ),1, wav_file_in );

			if( feof(wav_file_in) )
			{ 
				break ;
			}

			input = buf[0];		
			for( count = 0 ; count < MS_NUM_GRAINS ; count++){
				input  = input - ms_state[ count ] * ms_feedback[ count ];
				output = ms_state[ count ] * ms_feedforward[ count ];
                        }
			for( count = 0 ; count < MS_NUM_GRAINS ; count++){
			  if ( input > ms_threshold[ count ] ){
			    ms_state[ count ] = MS_MAX; 
                          } 
			  if ( input < ms_threshold[ count ] ){
			    ms_state[ count ] = MS_MIN; 
			  }
                        }

	                buf[0] = output;
			if ( wav_file ){ fwrite(&buf[0], sizeof( short ),1, wav_file); }
			  
		};	
		if ( wav_file    ){ fclose( wav_file );    }
		if ( wav_file_in ){ fclose( wav_file_in ); }
	}
}

/*
Note that the above code is not a full program, just a tiny library routine that can be used to generate WAV files. Here is a test program to generate sine waves (inefficiently) to see if the program works. Once you have the basic idea, you can write arbitrary programs to generate WAV files and play them back with QuickTime Player, RealPlayer, Audacity, or any of several other sound-playing programs.  This is really old-school computer music, from the days when computers were so slow that quality sounds could not be generated in real time.  It provides a good way to get started in learning the algorithms, though, as you don’t need to worry about efficiency or all the details needed to compute samples and provide them to a DAC at exactly the right rate.
*/

/* test_make_wav.c
 * Fri Jun 18 17:13:19 PDT 2010 Kevin Karplus
 * Test program for the write_wav function in make_wav.c
 */
 
int main(int argc, char ** argv)
{ 
	int i;
    if( argc < 2 ){ 
      printf("\nThis puts a wav header onto raw audio in 8000ms mono little endian\n usage: wav5 inputRawfilename outputWavFilename \n");  
    }
    char * ipFileName = "sine1.raw";
    char * opFileName = "sine1__schmitt_mag.raw";
    if( argc > 1 ){ ipFileName = argv[1]; }
    if( argc > 2 ){ opFileName = argv[2]; }  
    printf("Using: %s > %s \n",ipFileName,opFileName  );
	
#ifdef wav	
   /* wave file
	HEADER
	"data"+length
	data
	if the data length is not known the length needs to be updated. 
	*/
    //write_wav_header( "pcmtx_dhrdlin30.wav",  BUF_SIZE, buffer, S_RATE );
    write_wav_header( opFileName,  BUF_SIZE, buffer, S_RATE );
    /* NOTE: the function below needs to append not write op file */	

    copy_wav_data(opFileName,ipFileName); /* copy output input */

    /* update length */	
    write_wav_length( opFileName,    BUF_SIZE, buffer, S_RATE );
#endif

	
    magnetic_simulate_schmitt( opFileName, ipFileName );
    return 0;
}

