From f0d94272ccbb005ef20ddf9140fa3b4fa4fcf7c9 Mon Sep 17 00:00:00 2001
From: "sgunderson@bigfoot.com" <>
Date: Thu, 28 Dec 2006 12:41:48 +0100
Subject: [PATCH] Add the skew estimator.

---
 estimate-skew.c | 88 +++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 88 insertions(+)
 create mode 100644 estimate-skew.c

diff --git a/estimate-skew.c b/estimate-skew.c
new file mode 100644
index 0000000..1c44f86
--- /dev/null
+++ b/estimate-skew.c
@@ -0,0 +1,88 @@
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <stdlib.h>
+
+double frac(double x)
+{
+	return x - (int)x;
+}
+
+double filter(double x)
+{
+	static double l = 0.0;
+	l = 0.001 * x + 0.999 * l;
+	return l;
+}
+
+int main(int argc, char **argv)
+{
+	FILE *in1, *in2, *skew;
+
+	// open input1 (reference)
+	if (strcmp(argv[1], "-") == 0) {
+		in1 = stdin;
+	} else {
+		in1 = fopen(argv[1], "rb");
+		if (in1 == NULL) {
+			perror(argv[1]);
+			exit(1);
+		}
+	}
+
+	// open input2
+	if (strcmp(argv[2], "-") == 0) {
+		in2 = stdin;
+	} else {
+		in2 = fopen(argv[2], "rb");
+		if (in2 == NULL) {
+			perror(argv[2]);
+			exit(1);
+		}
+	}
+
+	// open (estimated) skew
+	if (strcmp(argv[3], "-") == 0) {
+		skew = stderr;
+	} else {
+		skew = fopen(argv[3], "wb");
+		if (skew == NULL) {
+			perror(argv[3]);
+			exit(1);
+		}
+	}
+
+	short prev_sample, sample = 0;
+	int in_pos = -1;
+	double p = -1.0;
+	double speed = 1.0;
+
+	while (!feof(in1) && !feof(in2)) {
+		// read sample from reference
+		if (fread(&sample, sizeof(short), 1, in1) != 1) {
+			exit(0);
+		}
+
+		// read sample from in2 (at current speed)
+		p += speed;
+		while ((int)(ceil(p)) > in_pos) {
+			prev_sample = sample;
+			if (fread(&sample, sizeof(short), 1, in2) != 1) {
+				exit(0);
+			}
+			++in_pos;
+		}
+
+		// linear interpolation (works well since delta_p varies so slowly)
+		double t = frac(p);
+		short intp_sample = prev_sample * (1.0 - t) + sample * t;
+
+		// subtract the two samples and pull it through a filter. we assume
+		// the sines are normalized for now (and that there's no dc skew)
+		double offset = -0.001 * filter(intp_sample - sample);
+		speed += offset;
+		fwrite(&speed, sizeof(double), 1, skew);
+
+		printf("%f (offset=%f)\n", speed, offset);
+	}
+}
-- 
2.39.2