1 files changed, 140 insertions, 0 deletions
diff --git a/src/moosynth.c b/src/moosynth.c
new file mode 100644
index 0000000..6d78415
--- /dev/null
+++ b/src/moosynth.c
@@ -0,0 +1,140 @@
+/*
+ * 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;
+	
+}
+