-#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 = 0.0;
- double speed = 1.001;
-
- while (!feof(in1) && !feof(in2)) {
- short refs;
-
- // read sample from reference
- if (fread(&refs, 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);
- double 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 = filter(intp_sample - refs);
- speed += 1e-8 * offset;
- fwrite(&speed, sizeof(double), 1, skew);
-
- printf("%f (offset=%f / filt=%f)\n", speed, intp_sample - refs, offset);
- }
-}