//
// hdev1 - Demo program to compute several Allan deviation statistics.
//
//    ADEV - Allan deviation (classic or overlapping)
//    MDEV - Modified Allan deviation
//    HDEV - Hadamard deviation (classic or overlapping)
//
// This was written primarily to explore the accuracy and speed of the
// algorithms. It matches exactly the results from Stable32 (the gold
// standard of commercial clock statistics programs). Input is a file
// of phase (time interval error) measurements. Output is one line or
// a table of Allan deviations (with sample counts in parenthesis).
//
// Several different styles of ADEV tables are possible:
//
//   oct  - octave tau (1, 2, 4, 8, ...)
//   dec  - decade tau (1, 10, 100, ...)
//   124  - 3 steps (1, 2, 4) per decade
//   125  - 3 steps (1, 2, 5) per decade
//   nine - 9 steps (1,2,3,4,5,6,7,8,9) per decade.
//   many - many tau (equal density of tau per decade)
//   all  - all tau (1, 2, 3, 4, 5, 6, ... N)
//
// Build (Win32) with cl /Ox (highly optimized).
//
// 28-Dec-2008 (tvb) Tom Van Baak
//

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#define MIN_ADEV 3  // sanity check for too large tau or too little data.

//
// Calculate Allan deviation (ADEV) from given phase array.
//  if ovlp == 0 then classic calculation (back-to-back intervals).
//  if ovlp == 1 then overlapping intervals.
//

double calc_adev (double data[], long count, long tau, long *terms, int ovlp)
{
    long i, n, stride;
    double sum, v;

    stride = ovlp ? 1 : tau;
    sum = n = 0;
    for (i = 0; (i + 2*tau) < count; i += stride) {
        v = data[i + 2*tau] - 2 * data[i + tau] + data[i];
        sum += v * v;
        n += 1;
    }
    sum /= 2.0;

    *terms = n;
    return (n < MIN_ADEV) ? 0.0 : sqrt(sum / n) / tau;
}

//
// Calculate Hadamard deviation (HDEV) from given phase array.
//  if ovlp == 0 then classic calculation (back-to-back intervals).
//  if ovlp == 1 then overlapping intervals.
//

double calc_hdev (double data[], long count, long tau, long *terms, int ovlp)
{
    long i, n, stride;
    double sum, v;

    stride = ovlp ? 1 : tau;
    sum = n = 0;
    for (i = 0; (i + 3*tau) < count; i += stride) {
        v = data[i + 3*tau] - 3*data[i + 2*tau] + 3*data[i + tau] - data[i];
        sum += v * v;
        n += 1;
    }
    sum /= 6.0;

    *terms = n;
    return (n < MIN_ADEV) ? 0.0 : sqrt(sum / n) / tau;
}

//
// Calculate Modified Allan deviation (MDEV) from given phase array.
// - This is a much faster (unrolled, single-loop) version.
// - Note long tau*tau can exceed 2^31 so it is promoted to double.
//

double calc_mdev_fast (double data[], long count, long tau, long *terms)
{
    long i, n;
    double sum, v;

    sum = v = n = 0;
    for (i = 0; (i + 2*tau) < count && i < tau; i += 1) {
        v += data[i + 2*tau] - 2 * data[i + tau] + data[i];
    }
    sum += v * v;
    n += 1;

    for (i = 0; (i + 3*tau) < count; i += 1) {
        v += data[i + 3*tau] - 3*data[i + 2*tau] + 3*data[i + tau] - data[i];
        sum += v * v;
        n += 1;
    }
    sum /= 2.0 * (double)tau * (double)tau;

    *terms = n;
    return (n < MIN_ADEV) ? 0.0 : sqrt(sum / n) / tau;
}

//
// Perform all three Allan calculations and display results.
//

int overlapping;    // do overlapping calculation
int calc_dp;        // decimal places to display
int n_dp = 6;       // display width of integer counts
int stop;           // hack to stop open-ended table calculations

void sigma_tau (double *data, long count, long tau)
{
    struct { double adev, mdev, hdev; } a;
    struct {  long  adev, mdev, hdev; } n;

    if (stop) return;

    a.adev = calc_adev(data, count, tau, &n.adev, overlapping);
    a.mdev = calc_mdev_fast(data, count, tau, &n.mdev);
    a.hdev = calc_hdev(data, count, tau, &n.hdev, overlapping);

 // stop = (a.adev == 0) || (a.mdev == 0) || (a.hdev == 0);
    stop = (a.adev == 0) && (a.mdev == 0) && (a.hdev == 0);
    if (stop) return;

    // Output in fixed-width table format, or

    if (calc_dp > 0) {
        printf("%*ld", n_dp, tau);
        printf(" %.*le (%*ld)", calc_dp, a.adev, n_dp, n.adev);
        printf(" %.*le (%*ld)", calc_dp, a.mdev, n_dp, n.mdev);
        printf(" %.*le (%*ld)", calc_dp, a.hdev, n_dp, n.hdev);
        printf("\n");
    }

    // Output in compact Excel-friendly CSV format.

    else {
        printf("%ld", tau);
        printf(",%.6le,%ld", a.adev, n.adev);
        printf(",%.6le,%ld", a.mdev, n.mdev);
        printf(",%.6le,%ld", a.hdev, n.hdev);
        printf("\n");
    }
}

int main (int argc, char *argv[])
{
    long count, tau, decade, clock_start, clock_end;
    double *phase, elapsed;
    double *read_data (FILE *f, long *count);
    char *user_tau;

    // 1. User specifies one tau, or keyword for table of sigma(tau).

    user_tau = (argc > 1) ? argv[1] : "124";
    fprintf(stderr, "(tau = '%s')\n", user_tau);

    // 2. User sets output precision (zero means CSV format).

    calc_dp = (argc > 2) ? atoi(argv[2]) : 4;
    if (calc_dp < 0) calc_dp = 0;
    fprintf(stderr, "(output %d decimal places)\n", calc_dp);

    if (argc > 2 && strchr(argv[2], ',')) {
        n_dp = atoi(strchr(argv[2], ',') + 1);
        fprintf(stderr, "(term counts %d digits)\n", n_dp);
    }

    // 3. Overlapping option

    overlapping = (argc > 3) ? atoi(argv[3]) : 1;
    fprintf(stderr, "(overlapping %s)\n", overlapping ? "yes" : "no");

    // Read phase data from stdin.

    phase = read_data(stdin, &count);
    fprintf(stderr, "(read %ld samples)\n", count);

    // Display heading and start timing the calculations.

    fprintf(stderr, "Displaying Tau, ADEV, MDEV, HDEV (each with term counts)\n");
    clock_start = clock();
    stop = 0;

    // Table: octaves of tau.

    if (strncmp(user_tau, "oct", 3) == 0) {
        for (tau = 1; !stop; tau *= 2) {
            sigma_tau(phase, count, tau);
        }
    } else

    // Table: decades of tau.

    if (strncmp(user_tau, "dec", 3) == 0) {
        for (tau = 1; !stop; tau *= 10) {
            sigma_tau(phase, count, tau);
        }
    } else

    // Table: all possible tau.

    if (strncmp(user_tau, "all", 3) == 0) {
        for (tau = 1; !stop; tau += 1) {
            sigma_tau(phase, count, tau);
        }
    } else

    // Table: many tau (equal density of tau per decade); e.g., m50

    if (strncmp(user_tau, "m", 1) == 0) {
        long i, otau, density;

        density = atoi(&user_tau[1]);
        if (density < 1) density = 1;
        otau = 0;
        for (decade = 1; !stop; decade *= 10) {
            for (i = 0; i < density; i += 1) {
                tau = (long) (decade * pow(10.0, (double) i / density));
                if (tau > otau) {
                    sigma_tau(phase, count, tau);
                    otau = tau;
                }
            }
        }
    } else

    // Table: nine tau (each digit) per decade.

    if (strncmp(user_tau, "nin", 3) == 0) {
        for (decade = 1; !stop; decade *= 10) {
            for (tau = 1; tau < 10; tau += 1) {
                sigma_tau(phase, count, tau * decade);
            }
        }
    } else

    // Table: 1-2-4 steps per decade.

    if (strcmp(user_tau, "124") == 0) {
        for (decade = 1; !stop; decade *= 10) {
            sigma_tau(phase, count, 1 * decade);
            sigma_tau(phase, count, 2 * decade);
            sigma_tau(phase, count, 4 * decade);
        }
    } else

    // Table: 1-2-5 steps per decade.

    if (strcmp(user_tau, "125") == 0) {
        for (decade = 1; !stop; decade *= 10) {
            sigma_tau(phase, count, 1 * decade);
            sigma_tau(phase, count, 2 * decade);
            sigma_tau(phase, count, 5 * decade);
        }
    } else

    // One user-specified tau number.

    if (user_tau[0] > '0' && user_tau[0] <= '9') {
        tau = atoi(user_tau);
        sigma_tau(phase, count, tau);
    } else

    {
        fprintf(stderr, "invalid tau: %s\n", user_tau);
        exit(1);
    }

    // Display elapsed time.

    clock_end = clock();
    elapsed = (double)(clock_end - clock_start) / CLOCKS_PER_SEC;
    fprintf(stderr, "(elapsed time %.2lf seconds)\n", elapsed);

    free(phase);
    return 0;   // exit(0);
}

//
// Read file into memory.
// - One time interval reading per line, pre-scaled to seconds units.
// - Sample data (e.g., from a time interval counter on GPS receiver):
//   1.6862E-08
//   2.6083E-08
//   5.1941E-09
//   2.0985E-08
//   9.2359E-09
//   ...
//

double *read_data (FILE *f, long *count)
{
    double *data;
    long index, maxcount, size;
    char line[1000];

    data = NULL;
    index = maxcount = 0;
    while (fgets(line, sizeof line, f) != NULL) {
        if (maxcount == 0) {
            maxcount = 1024;
            size = maxcount * sizeof data[0];
            data = malloc(size);

        } else if (index == maxcount) {
            maxcount *= 2;
            size = maxcount * sizeof data[0];
            data = realloc(data, size);
        }
        if (data == NULL) {
            fprintf(stderr, "malloc failed (%ld bytes)\n", size);
            exit(1);
        }
        data[index] = atof(line);
        index += 1;
    }
    *count = index;
    return data;
}