Add the skew estimator.

diff --git a/estimate-skew.c b/estimate-skew.c
new file mode 100644 (file)
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);
+       }
+}