refactor sources to subdir

Signed-off-by: Andy Green <andy.green@linaro.org>

3 files changed, 905 insertions, 0 deletions

diff --git a/lava-fft/Makefile.am b/lava-fft/Makefile.am
new file mode 100644
index 0000000..0021fd7
--- /dev/null
+++ b/lava-fft/Makefile.am
@@ -0,0 +1,7 @@
+bin_PROGRAMS=lava-fft
+lava_fft_SOURCES=lava-fft.c
+lava_fft_CFLAGS=-fPIC -Wall -Werror -D_FORTIFY_SOURCE=2 -fstack-protector -std=gnu99 -pedantic -DINSTALL_DATADIR=\"${datarootdir}\"
+lava_fft_LDFLAGS=-fPIC
+lava_fft_LDADD=-lfftw3 -lm
+
+
diff --git a/lava-fft/compare.gp b/lava-fft/compare.gp
new file mode 100644
index 0000000..8a35cf7
--- /dev/null
+++ b/lava-fft/compare.gp
@@ -0,0 +1,18 @@
+set terminal pngcairo notransparent enhanced font "arial,10" size 760, 600 
+set output 'plot.png'
+
+set title "Test file vs capture comparison" font ",16" 
+#set title "Shunt voltage measurement nonlinearlity Error % (0-100mV)" font ",16" 
+set yrange [  : ] noreverse nowriteback
+set lmargin  9
+set rmargin  9
+
+set ylabel 'magnitude' tc lt 1
+set xlabel 'Frequency (Hz)' tc lt 1
+
+set autoscale xfixmin
+#set arrow from graph 0, 0 to 20000, 20000 nohead ls 1 lc rgb 'red'
+
+plot 'g' using 1:3 title 'Original' with lines lc rgb "#ff0000", \
+	'' using 1:(- $2) title 'Captured' with lines lc rgb "#0000ff";
+
diff --git a/lava-fft/lava-fft.c b/lava-fft/lava-fft.c
new file mode 100644
index 0000000..ed857e2
--- /dev/null
+++ b/lava-fft/lava-fft.c
@@ -0,0 +1,880 @@
+/*
+ * fft.c: utility to generate test wav file on stdout or assess a figure-of
+ * merit from a wav file on stdin compared to the original WAV file.
+ * Assesses in frequency domain.
+ *
+ * Author: Andy Green <andy.green@linaro.org> 
+ * Copyright (C) 2012 Linaro, LTD
+ *
+ * This program 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 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program 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 this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+ * 02110-1301, USA.
+ *
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <unistd.h>
+#include <fftw3.h>
+#include <math.h>
+#include <stdint.h>
+#include <termio.h>
+#include <getopt.h>
+#include <string.h>
+
+enum audio_formats {
+	WAF_PCM = 1,
+	WAF_MULAW = 6,
+	WAF_ALAW = 7,
+	WAF_IBMMULAW = 257,
+	WAF_IBMALAW = 258,
+	WAF_ADPCM = 259,
+};
+
+struct wav_header {
+	char		RIFF[4];
+	uint32_t	chunk_size  __attribute__((packed));
+	char		WAVE[4];
+	char		fmt[4];
+	uint32_t	fmt_chunk_size  __attribute__((packed));
+	uint16_t	audio_format  __attribute__((packed));
+	uint16_t	count_channels  __attribute__((packed));
+	uint32_t	samples_per_second  __attribute__((packed));
+	uint32_t	bytes_per_second  __attribute__((packed));
+	uint16_t	block_align  __attribute__((packed));
+	uint16_t	bits_per_sample  __attribute__((packed));
+	char		data[4];
+	uint32_t	data_length  __attribute__((packed));
+}  __attribute__((packed));
+
+#define DURATION_SAMPLES 65536
+#define RATE 48000
+#define CHANNELS 2
+#define BYTES_PER_SAMPLE 2
+#define MAX_SUPPORTED_CHANNELS 8
+#define MIN_TEST_AMPLITUDE 0x800
+#define SILENCE_FFT_LIMIT 200000
+#define IGNORE_FOR_SILENCE_DETECT 0.05
+
+#define SIN_MODULATION_HZ (RATE / 67)
+#define SIN_CHAN_STEP_HZ 133
+#define SIN_MODULATION_FACTOR 0.6
+#define SQUARE_MODULATION_HZ 1000
+#define SQUARE_CHAN_STEP_HZ 33
+#define SQUARE_MODULATION_FACTOR 0.4
+#define TRI_HZ ((double)RATE / 37)
+
+
+struct wav_header wav = {
+	.RIFF = { 'R', 'I', 'F', 'F' },
+	.WAVE = { 'W', 'A', 'V', 'E' },
+	.fmt = { 'f', 'm', 't', ' ' },
+	.data = { 'd', 'a', 't', 'a' },
+	.audio_format = WAF_PCM,
+	.count_channels = CHANNELS,
+	.samples_per_second = RATE,
+	.bytes_per_second = RATE * CHANNELS * BYTES_PER_SAMPLE,
+	.block_align = CHANNELS * BYTES_PER_SAMPLE,
+	.bits_per_sample = 16,
+	.chunk_size = sizeof(struct wav_header) - 4 +
+			(2 *DURATION_SAMPLES * CHANNELS * BYTES_PER_SAMPLE),
+	.data_length = 2 *DURATION_SAMPLES * CHANNELS * BYTES_PER_SAMPLE,
+	.fmt_chunk_size = 0x10,
+};
+
+int check_header_sanity(struct wav_header *wavh)
+{
+	if (wavh->audio_format != 1) {
+		fprintf(stderr, "Need PCM format WAV\n");
+		return 3;
+	}
+
+	if (
+		wavh->RIFF[0] != 'R' || wavh->RIFF[1] != 'I' ||
+		wavh->RIFF[2] != 'F' || wavh->RIFF[3] != 'F' ||
+		wavh->WAVE[0] != 'W' || wavh->WAVE[1] != 'A' ||
+		wavh->WAVE[2] != 'V' || wavh->WAVE[3] != 'E' ||
+		wavh->fmt[0]  != 'f' || wavh->fmt[1]  != 'm' ||
+		wavh->fmt[2]  != 't' || wavh->fmt[3]  != ' ' ||
+		wavh->data[0] != 'd' || wavh->data[1] != 'a' ||
+		wavh->data[2] != 't' || wavh->data[3] != 'a'
+	) {
+		fprintf(stderr, "problem with WAV header format\n");
+		return 4;
+	}
+
+	return 0;
+}
+
+long sample(int waveform, int fault, double peak_amplitude,
+				unsigned int index, unsigned int channel)
+{
+	long l;
+	double m0, m1;
+
+	/* generate the test waveform */
+
+	switch (waveform) {
+	case 0:
+	default:
+		m0 = (double)wav.samples_per_second /
+			(SQUARE_MODULATION_HZ +
+				(SQUARE_CHAN_STEP_HZ * channel));
+		m1 = ((double)index / m0);
+		m1 -= (long)m1;
+		l = (peak_amplitude * 32767.0 *
+				SQUARE_MODULATION_FACTOR);
+		if (m1 >= 0.5)
+			l = -l;
+		l += (peak_amplitude * 32767.0 *
+			SIN_MODULATION_FACTOR *
+			sin((double)index * ((2 * M_PI) /
+			((double)wav.samples_per_second /
+			(SIN_MODULATION_HZ +
+				    (SIN_CHAN_STEP_HZ * channel))))));
+		break;
+	case 1: /* TRI_HZ = eg, 1297  ^v*/
+		m0 = (double)wav.samples_per_second /
+			      (TRI_HZ + (SIN_CHAN_STEP_HZ * channel));
+		/*
+		 * start at 90 degree phase offset to get a
+		 * good starting detect
+		 */
+		m1 = ((double)index / m0) + 0.25; 
+		m1 -= (long)m1;
+		if (m1 >= 0.5)
+			if (m1 >= 0.75)
+				m1 = (-peak_amplitude *
+					32767.0) *
+				       (0.25 - (m1 - 0.75)) * 4;
+			else
+				m1 = (-peak_amplitude *
+					32767.0) *
+						 (m1 - 0.5) * 4;
+		else
+			if (m1 >= 0.25)
+				m1 = (peak_amplitude *
+					32767.0) *
+				       (0.25 - (m1 - 0.25)) * 4;
+			else
+				m1 = (peak_amplitude *
+					32767.0) * (m1) * 4;
+		l = m1;
+		break;
+	}
+
+	switch (fault) {
+		case 1: /* random samples / 16 .. */
+		case 2:
+		case 3:
+		case 4:
+		case 5: /* 256 */
+		case 6:
+		case 7:
+		case 8:
+		case 9:
+		case 10: /* 8192 */
+		case 11: /* 16384 */
+		case 12: /* 32768 */
+		case 13: /* 64K / per 1.6s @ 48ksps */
+		case 14: /* 128K / per 3.2s @ 48ksps */
+		case 15: /* 256K / per 6.4s @ 48ksps */
+			if ((index & ((1 << (fault + 3)) - 1)) ==
+					  ((1 << (fault + 2)))) {
+				l = (((double)rand() /
+				 (RAND_MAX / 2)) - 1.0) * 32767;
+//				fprintf(stderr, "%ld\n", l);
+			}
+			break;
+		case 0:
+		default:
+			break;
+	}
+
+	return l;
+}
+
+
+
+static struct option options[] = {
+	{ "help",	no_argument,		NULL, 'h' },
+	{ "verbose",	no_argument,		NULL, 'v' }, 
+	{ "rate",	required_argument,	NULL, 'r' },
+	{ "test",	required_argument,	NULL, 't' },
+	{ "noise",	required_argument,	NULL, 'n' },
+	{ "attenuate",	required_argument,	NULL, 'a' },
+	{ "rolloff", 	required_argument,	NULL, 'o' },
+	{ "channels", 	required_argument,	NULL, 'c' },
+	{ "snip", 	required_argument,	NULL, 's' },
+	{ "gnuplot", 	required_argument,	NULL, 'g' },
+	{ "fault", 	required_argument,	NULL, 'f' },
+	{ "waveform", 	required_argument,	NULL, 'w' },
+	{ "ignore-silence", no_argument,	NULL, 'i' },
+	{ "length",	required_argument,	NULL, 'l' },
+	{ NULL, 0, 0, 0 }
+};
+
+int main( int argc, char** argv )
+{
+	fftw_complex    *ref_data[MAX_SUPPORTED_CHANNELS],
+			*data[MAX_SUPPORTED_CHANNELS],
+			*fft_result = NULL,
+			*ref_fft_result[MAX_SUPPORTED_CHANNELS],
+			*fft_background[MAX_SUPPORTED_CHANNELS];
+	fftw_plan       plan_forward[MAX_SUPPORTED_CHANNELS],
+			plan[MAX_SUPPORTED_CHANNELS],
+				bg_plan[MAX_SUPPORTED_CHANNELS];
+	int i, j;
+	int SIZE = DURATION_SAMPLES;
+	int ISSUE_SIZE = DURATION_SAMPLES * 2;
+	int fd = 1; /* stdout */
+	short buf[256];
+	int32_t *buf_32 = (int32_t *)&buf[0];
+	int ring[MAX_SUPPORTED_CHANNELS][DURATION_SAMPLES];
+	int head = 0;
+	int pos = 0;
+	int compare = 0;
+	int test = 3;
+	double noise = 0.0;
+	int n = 1;
+	double r;
+	long l = 0;
+	double peak_amplitude = 0.8;
+	struct wav_header cap_wav;
+	double delta, m0, m1, m2;
+	int skip = 0;
+	double attenuate = 1.0;
+	double ref_largest_ac_mag = 0.0, fdomain_scale, d;
+	int ref_largest_ac_index = 0;
+	double derate_freq = 99999999.0;
+	int derate_index;
+	double derate;
+	double fom = 0;
+	double snip = 0;
+	double leadin;
+	int gnuplot = 0;
+	int gnuplot_channel = 0;
+	double m;
+	int fault = 0;
+	int waveform = 0;
+	int ignore_silence = 0;
+	double length;
+
+	if (sizeof(struct wav_header) != 44) {
+		fprintf(stderr, "Packing problem with struct wav_header\n");
+		return -1;
+	}
+
+	length = (2.0 * (double)DURATION_SAMPLES) / wav.samples_per_second;
+
+	while (n >= 0) {
+		n = getopt_long(argc, argv, "hvr:t:n:a:o:c:s:g:f:w:il:",
+								 options, NULL);
+		if (n < 0)
+			continue;
+		switch (n) {
+
+		case 'l':
+			length = atof(optarg);
+			break;
+
+		case 'i':
+			ignore_silence = 1;
+			break;
+
+		case 's':
+			snip = atof(optarg);
+			break;
+
+		case 'g':
+			gnuplot = 1;
+			gnuplot_channel = atoi(optarg);
+			break;
+
+		case 'c':
+			wav.count_channels = atoi(optarg);
+			wav.bytes_per_second = RATE * wav.count_channels *
+							      BYTES_PER_SAMPLE;
+			wav.block_align = wav.count_channels * BYTES_PER_SAMPLE;
+			wav.chunk_size = sizeof(struct wav_header) - 4 +
+				(DURATION_SAMPLES * wav.count_channels *
+							      BYTES_PER_SAMPLE);
+			wav.data_length = DURATION_SAMPLES *
+					wav.count_channels * BYTES_PER_SAMPLE;
+
+			if (wav.count_channels >= MAX_SUPPORTED_CHANNELS) {
+				fprintf(stderr, "Can't handle %d channels\n",
+							    wav.count_channels);
+				return 9;
+			}
+			break;
+
+		case 'a':
+			attenuate = atof(optarg);
+			break;
+
+		case 'o':
+			derate_freq = atof(optarg);
+			break;
+
+		case 'r':
+			wav.samples_per_second = atoi(optarg);
+			break;
+
+		case 'w':
+			waveform = atoi(optarg);
+			break;
+
+		case 't':
+			test = atoi(optarg);
+			break;
+
+		case 'n':
+			noise = atof(optarg);
+			if (noise < 0 || noise > 1.0) {
+				fprintf(stderr, "mixed noise factor can only "
+						"be a float between 0 and 1\n");
+				return 11;
+			}
+			break;
+
+		case 'f':
+			fault = atoi(optarg);
+			break;
+
+		default:
+		case 'h':
+			fprintf(stderr,
+				"Usage: either... lava-fft > wav_file   or\n"
+					"   cat capture | lava-fft\n"
+				" [--verbose -v] Increase debug on stderr\n"
+				" [--rate -r <sample rate>] set sample rate "
+					"(default 48000)\n"
+				" [--test -t <bitfield>] which channels have "
+					"signal, ch0 = 1, ch1 = 2 etc "
+					"or'd together\n"
+				" [--rolloff -o <frequency> start to linearly "
+					"derate results past this frequency "
+					"down to 0 importance at Nyquist\n"
+				" [--snip - <secs (float)> time to ignore "
+					"after first sound energy seen "
+					"in capture, eg -s 0.25\n"
+				" [ --gnuplot -g <channel index>  issue graph "
+					"data on stdout suitable for gnuplot, "
+					"for channel <channel index>"
+				" [ --fault -f <index>  inject fault...\n"
+				"     1 = random sample every 16\n"
+				"     2 = random sample every 32\n"
+				"     3 = random sample every 64\n"
+				"     4 = random sample every 128\n"
+				"      ...\n"
+				"    13 = random sample every 65536\n"
+				"   |16 = swap channel pairs\n"
+				" [--ignore-silence -i] Do not check channels"
+					" that are disabled in --test "
+					"bitfield... normally they get "
+					"checked to have been silent\n"
+				" [--length -l <secs>] duration of test file\n"
+
+			);
+			return 1;
+		}
+	}
+
+	/* how many samples to jump over after noise detected */
+	snip = snip * wav.samples_per_second;
+
+	/* 100ms + sample duration of silence */
+
+	leadin = ((1.0 / wav.samples_per_second) * DURATION_SAMPLES) + 0.1;
+
+	if (!isatty(0)) {
+		compare = 1;
+
+		if (read(0, &cap_wav, sizeof wav) != sizeof wav) {
+			fprintf(stderr, "Problem reading WAV header\n");
+			return 2;
+		}
+
+		n = check_header_sanity(&wav);
+		if (n)
+			return n;
+
+		if (wav.samples_per_second != cap_wav.samples_per_second) {
+			fprintf(stderr, "Need capture (%u sps) and reference "
+				"(%u sps) WAVs to have the same sample rate\n",
+			    wav.samples_per_second, cap_wav.samples_per_second);
+			return 20;
+		}
+
+		if (wav.count_channels != cap_wav.count_channels) {
+			fprintf(stderr, "Need capture (%u) and reference (%u) "
+				"WAVs to have the same channel count\n",
+				wav.count_channels, cap_wav.count_channels);
+			return 20;
+		}
+
+		switch (cap_wav.bits_per_sample) {
+		case 16:
+			break;
+		case 32:
+			break;
+		default:
+			fprintf(stderr, "Sorry %dbpp not supported\n",
+				cap_wav.bits_per_sample);
+			return 52;
+		}
+	}
+
+	/* allocate FFT array for ref data, which is always needed */
+
+	for (n = 0; n < wav.count_channels; n++)
+		ref_data[n] = (fftw_complex *)fftw_malloc(
+					     sizeof(fftw_complex) * ISSUE_SIZE);
+
+	if (derate_freq > wav.samples_per_second / 2)
+		derate_freq = wav.samples_per_second / 2;
+
+	derate_index = (derate_freq / (wav.samples_per_second / 2)) *
+							DURATION_SAMPLES / 2;
+
+	/*
+	 * always synthesize the test data
+	 *
+	 * waveform 0 is a sine wave modulated with a square wave
+	 * the sine is necessary to detect clipping, since a clipped
+	 * square wave looks pretty much the same ^^
+	 *
+	 * frequencies were selected to be far enough away from DC to have
+	 * reasonable level after capture and separate enough from square
+	 * harmonics to tell them apart
+	 *
+	 * waveform 1 is a triangle wave
+	 *
+	 * we generate a differently pitched version of the selected waveform
+	 * for each channel so they can't be mistaken for each other.
+	 */
+
+	srand(1);
+	for (n = 0; n < wav.count_channels; n++) {
+		for (i = 0 ; i < ISSUE_SIZE ; i++ ) {
+
+			l = sample(waveform, fault, peak_amplitude, i, n);
+
+			if (fault & 16) {
+				ref_data[n ^ 1][i][0] = l;
+				ref_data[n ^ 1][i][1] = 0.0;
+			} else {
+				ref_data[n][i][0] = l;
+				ref_data[n][i][1] = 0.0;
+			}
+		}
+	}
+
+
+	/*
+	 * if we're writing out the synthetic data, just do that and be done
+	 */
+
+	if (!compare) {
+
+		wav.bytes_per_second = wav.samples_per_second *
+					wav.count_channels * BYTES_PER_SAMPLE;
+		wav.block_align = wav.count_channels * BYTES_PER_SAMPLE;
+		wav.chunk_size = sizeof(struct wav_header) - 4 +
+			(2 * (length * wav.samples_per_second) *
+			wav.count_channels * BYTES_PER_SAMPLE);
+		wav.data_length = 2 * (length * wav.samples_per_second) *
+			wav.count_channels * BYTES_PER_SAMPLE;
+
+		if (write(fd, &wav, sizeof wav) != sizeof wav) {
+			fprintf(stderr, "Problem writing output\n");
+			return 40;
+		}
+
+		/*
+		 * write the silent leadin... initial output is funky on some
+		 * soundcards
+		 */
+
+		for (n = 0; n < wav.count_channels; n++)
+			buf[n] = 0;
+
+
+		for (i = 0; i < ((leadin + IGNORE_FOR_SILENCE_DETECT) *
+						wav.samples_per_second); i++)
+			if (write(fd, &buf[0],
+					sizeof buf[0] * wav.count_channels) !=
+					   sizeof buf[0] * wav.count_channels) {
+				fprintf(stderr, "Problem writing to output\n");
+				return 41;
+			}
+
+		srand(1);
+		for (i = 0; i < (length * wav.samples_per_second); i++) {
+
+			for (n = 0; n < wav.count_channels; n++) {
+
+				if (fault & 16)
+					j = n ^ 1;
+				else
+					j = n;
+
+				if (test & (1 << n))
+					buf[pos++] = sample(waveform, fault,
+							  peak_amplitude, i, j);
+				else
+					buf[pos++] = 0;
+			}
+
+			if (pos >= sizeof buf / sizeof buf[0]) {
+				pos = 0;
+				if (write(fd, &buf[0], sizeof buf) !=
+								   sizeof buf) {
+					fprintf(stderr,
+						"Problem writing to output\n");
+					return 41;
+				}
+			}
+		}
+
+		if (pos)
+			if (write(fd, &buf[0],
+				pos * sizeof(buf[0]) * wav.count_channels) != 
+				  (pos * sizeof(buf[0]) * wav.count_channels)) {
+				fprintf(stderr, "Problem writing to output\n");
+				return 42;
+			}
+		if (fd > 2)
+			close(fd);
+
+		goto done;
+	}
+
+	/*
+	 * allocate FFT arrays for analysis
+	 */
+
+	fft_result = (fftw_complex *)fftw_malloc(
+					sizeof(fftw_complex) * SIZE);
+
+	for (i = 0; i < wav.count_channels; i++) {
+		data[i] = (fftw_complex *)fftw_malloc(
+					sizeof(fftw_complex) * SIZE);
+		plan[i] = fftw_plan_dft_1d(SIZE, data[i], fft_result,
+						   FFTW_FORWARD, FFTW_ESTIMATE);
+		fft_background[i] = (fftw_complex *)fftw_malloc(
+					sizeof(fftw_complex) * SIZE);
+		bg_plan[i] = fftw_plan_dft_1d(SIZE, data[i], fft_background[i],
+						   FFTW_FORWARD, FFTW_ESTIMATE);
+
+		ref_fft_result[i] = (fftw_complex *)fftw_malloc(
+					sizeof(fftw_complex) * SIZE);
+
+		plan_forward[i]  = fftw_plan_dft_1d(SIZE, ref_data[i],
+				ref_fft_result[i], FFTW_FORWARD, FFTW_ESTIMATE);
+
+		/* compute synthesized reference FFTs */
+
+		fftw_execute(plan_forward[i]);
+		fftw_destroy_plan(plan_forward[i]);
+	}
+
+	/*
+	 * we are in comparing mode
+	 */
+
+
+	/*
+	 * trim stdin data until we have meaningful amplitude on
+	 * one or the other channel (trim silence)
+	 */
+
+	n = 1;
+	while (n) {
+
+		if (read(0, &buf[0], cap_wav.block_align) !=
+							cap_wav.block_align) {
+			fprintf(stderr, "Unable to find start of test "
+							  "sound before EOF\n");
+			return 30;
+		}
+
+		if (test == 0) {
+			n = 0;
+			continue;
+		}
+
+		for (i = 0; i < cap_wav.count_channels; i++) {
+			switch (cap_wav.bits_per_sample) {
+			case 16:
+				l = buf[i];
+				break;
+			case 32:
+				l = buf_32[i] >> 16;
+				break;
+				
+			}
+			if (l > MIN_TEST_AMPLITUDE ||
+					l < -MIN_TEST_AMPLITUDE)
+				if (skip > (IGNORE_FOR_SILENCE_DETECT *
+						(double)wav.samples_per_second))
+					n = 0;
+		}
+
+		if (n) {
+			skip++;
+			for (i = 0; i < cap_wav.count_channels; i++) {
+				switch (cap_wav.bits_per_sample) {
+				case 16:
+					l = buf[i];
+					break;
+				case 32:
+					l = buf_32[i] >> 16;
+					break;
+				
+				}
+				ring[i][head] = l;
+			}
+			head++;
+			if (head == DURATION_SAMPLES)
+				head = 0;
+		}
+
+	}
+
+	if (skip < SIZE) {
+		fprintf(stderr, "Failed to see enough leadin silence %d %ld\n",
+								       skip, l);
+		return 50;
+	}
+
+
+	/* set up datas with silence from ringbuffer */
+
+	for (n = 0 ; n < SIZE ; n++) {
+
+		for (i = 0; i < cap_wav.count_channels; i++) {
+			data[i][n][0] = ring[i][head];
+			data[i][n][1] = 0.0;
+		}
+
+		head++;
+		if (head == DURATION_SAMPLES)
+				head = 0;
+	}
+
+	/* do fft on "silence" for each channel */
+
+	for (i = 0; i < cap_wav.count_channels; i++)
+		fftw_execute(bg_plan[i]);
+
+
+	/* trim the additional requested snip period */
+
+	while (snip > 0) {
+
+		skip++;
+
+		if (read(0, &buf[0], cap_wav.block_align) !=
+							  cap_wav.block_align) {
+			fprintf(stderr, "Unable to find start of test sound "
+								"before EOF\n");
+			return 30;
+		}
+
+		snip = snip - 1;
+	}
+
+	fprintf(stderr, "Skipping to +%fs in capture\n",
+		(double)skip / (double)wav.samples_per_second);
+
+
+	/* set up datas with stdin content */
+	/* inherit last buf set from snip code */
+
+	for (n = 0 ; n < SIZE ; n++) {
+
+		for (i = 0; i < cap_wav.count_channels; i++) {
+
+			r = ((double)(rand() - (RAND_MAX / 2)) /
+			     (double)(RAND_MAX / 2)) * peak_amplitude * 32767.0;
+
+			switch (cap_wav.bits_per_sample) {
+			case 16:
+				l = (noise * r) + ((1.0 - noise) *
+					(((double)buf[i]) * attenuate));
+				break;
+			case 32:
+				l = (noise * r) + ((1.0 - noise) *
+					(((double)(buf_32[i] >> 16)) *
+								    attenuate));
+				break;
+				
+			}
+
+			data[i][n][0] = l;
+			data[i][n][1] = 0.0;
+		}
+
+		if (n == SIZE - 1)
+			continue;
+
+		if (read(0, &buf[0], cap_wav.block_align) !=
+							cap_wav.block_align)
+			fprintf(stderr, "(Ran out of captured samples)\n");
+	}
+
+
+	if (gnuplot)
+		printf("time\tcaptured\toriginal\n");
+
+	/* go through each channel and assess against synthesized one */
+
+	fom = 0;
+
+	for (i = 0; i < cap_wav.count_channels; i++) {
+
+		if (!(test & (1 << i)) && ignore_silence)
+			continue;
+
+		if (fault & 16)
+			j = i ^ 1;
+		else
+			j = i;
+
+		/* compute captured channel FFT into fft_result */
+
+		fftw_execute(plan[i]);
+
+		/* What's the highest peak that we sent on this channel? */
+
+		ref_largest_ac_mag = 0;
+		ref_largest_ac_index = 0;
+		for (n = 1; n < SIZE / 2; n++) {
+			d = sqrt((ref_fft_result[i][n][0] *
+				  ref_fft_result[i][n][0]) +
+				 (ref_fft_result[i][n][1] *
+				  ref_fft_result[i][n][1]));
+			if (d > ref_largest_ac_mag) {
+				ref_largest_ac_mag = d;
+				ref_largest_ac_index = n;
+			}
+		}
+
+		/*
+		 * looking at the highest freq peak we sent out,
+		 * scale what came in on the frequency domain by the
+		 * ratio of what was seen for the highest energy slot compared
+		 * to what was sent
+		 */
+
+		d = sqrt((fft_result[ref_largest_ac_index][0] *
+				fft_result[ref_largest_ac_index][0]) +
+				(fft_result[ref_largest_ac_index][1] *
+				fft_result[ref_largest_ac_index][1]));
+		fdomain_scale = ref_largest_ac_mag / d;
+
+		/* see how well it compares to template in frequency domain */
+
+		delta = 0;
+		derate = 1.0;
+
+		if (d != 0) {  // because if the biggest magnitude was 0...
+
+			for (n = 1; n < SIZE / 2; n++) {
+
+				if (n > derate_index)
+					derate = 1.0 -
+						((double)(n - derate_index) /
+					 ((double)((SIZE / 2) - derate_index)));
+
+				m0 = sqrt((fft_background[i][n][0] *
+					fft_background[i][n][0]) +
+					(fft_background[i][n][1] *
+					fft_background[i][n][1]));
+
+				m1 = sqrt(((fft_result[n][0] * fdomain_scale) *
+					(fft_result[n][0] * fdomain_scale)) +
+					((fft_result[n][1] * fdomain_scale) *
+					(fft_result[n][1] * fdomain_scale)));
+
+				/*
+				 * compare incoming to either sent template, or
+				 * silence if test disabled
+				 */
+
+				if (test & (1 << i))
+					m2 = derate *
+						sqrt((ref_fft_result[j][n][0] *
+						ref_fft_result[j][n][0]) +
+						(ref_fft_result[j][n][1] *
+						      ref_fft_result[j][n][1]));
+				else
+					m2 = 0;
+
+				if (gnuplot && i == gnuplot_channel)
+					printf("%f\t%f\t%f\n", (double)n *
+					       ((double)wav.samples_per_second /
+							 (double)SIZE), m1, m2);
+
+				m = (derate * (m1 - m0)) - m2;
+
+				/*
+				 * safety net to detect funny "silence"
+				 */
+
+				if (m0 > SILENCE_FFT_LIMIT) {
+					fprintf(stderr, "%f\n", m0);
+					m += m0 - m2;
+				}
+
+				if (m < 0)
+					m = -m;
+				m = m / ref_largest_ac_mag;
+
+				if (m > 0) {
+//					fprintf(stderr, "%f\n", m);
+					delta += m * 100;
+				}
+			}
+		}
+
+		if (delta / SIZE > fom)
+			fom = delta / SIZE;
+
+		fprintf(stderr, "Ch %d: %.6f\n", i, delta / SIZE);
+	}
+
+	if (!gnuplot)
+		printf("%.3f\n", fom);
+
+done:
+
+	if (compare) {
+		for (i = 0; i < wav.count_channels; i++) {
+			fftw_destroy_plan(plan[i]);
+			fftw_destroy_plan(bg_plan[i]);
+			fftw_free(data[i]);
+			fftw_free(fft_background[i]);
+			fftw_free(ref_fft_result[i]);
+		}
+		fftw_free(fft_result);
+	}
+	for (i = 0; i < wav.count_channels; i++)
+		fftw_free(ref_data[i]);
+
+	return 0;
+}
+