/*
 * moosynth.c
 *
 * Work out how to synthesize a cow
 *
 * (c) 2008 Thomas White <taw27@srcf.ucam.org>
 *
 * This file is part of OpenMooCow - accelerometer moobox simulator
 *
 * OpenMooCow is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * OpenMooCow is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with OpenMooCow.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdlib.h>
#include <stdarg.h>
#include <stdio.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <fftw3.h>
#include <math.h>

int main(int argc, char *argv[]) {

	FILE *fh;
	double *data;
	fftw_complex *ft;
	struct stat statbuf;
	fftw_plan plan;
	int samples;
	int16_t *sdata;
	int i, max_fq, c;
	double max_am, next_max_am, max_ph;
	const int ncomp = 20;
	int components[ncomp];
	
	if ( argc != 2 ) {
		fprintf(stderr, "Syntax: %s <file.pcm>\n", argv[0]);
		return 1;
	}
	
	/* Determine size of data */
	if ( stat(argv[1], &statbuf) == -1 ) {
		fprintf(stderr, "Couldn't stat file '%s'\n", argv[1]);
		return 1;
	}
	samples = statbuf.st_size/2;	/* Two bytes per sample, one channel */
	printf("File size: %lli bytes\n", (long long int)statbuf.st_size);
	
	/* Read the data in */
	fh = fopen(argv[1], "rb");
	if ( fh == NULL ) {
		fprintf(stderr, "Couldn't open file '%s'\n", argv[1]);
		return 1;
	}
	sdata = malloc(samples*2);
	fread(sdata, 2, samples, fh);
	fclose(fh);
	
	data = fftw_malloc(samples*sizeof(double));
	ft = fftw_malloc(samples*sizeof(fftw_complex));
	
	for ( i=0; i<samples; i++ ) {
		data[i] = sdata[i];
	}
	
	plan = fftw_plan_dft_r2c_1d(samples, data, ft, FFTW_ESTIMATE);
	fftw_execute(plan);
	
	printf("DC component: %f %f\n", ft[0][0], ft[0][1]);
	next_max_am = +1000000000.0;
	for ( c=0; c<ncomp-1; c++ ) {
	
		max_am = 0.0;
		max_ph = 0.0;
		max_fq = 0;
		for ( i=1; i<samples/2; i++ ) {
	
			double re, im, am, ph;
	
			re = ft[i][0];
			im = ft[i][1];
			am = sqrt(re*re + im*im);
			ph = atan2(im, re);
		//	printf("Frequency: %i Hz : %f %f\n", (44100/samples)*i, am, 180*(ph/M_PI));
	
			if ( (am > max_am) && (am < next_max_am) ) {
				max_am = am;
				max_fq = i;
				max_ph = ph;
			}
	
		}
		printf("%2i %6i Hz %14.2f %+7.1f\n", c+1, (44100/samples)*max_fq, max_am, 180*(max_ph/M_PI));
		components[c] = max_fq;
		next_max_am = max_am;
		
	}
	
	for ( i=0; i<samples/2; i++ ) {
		int c;
		int found = 0;
		for ( c=0; c<ncomp-1; c++ ) {
			if ( components[c] == i ) found = 1;
		}
		if ( !found ) {
			ft[i][0] = 0.0;
			ft[i][1] = 0.0;
		}
	}
	
	plan = fftw_plan_dft_c2r_1d(samples, ft, data, FFTW_ESTIMATE);
	fftw_execute(plan);
	
	for ( i=0; i<samples; i++ ) {
		sdata[i] = (data[i])/(samples*2);
	}
	fh = fopen("filtered.pcm", "wb");
	fwrite(sdata, 2, samples, fh);
	fclose(fh);
	
	return 0;
	
}