home *** CD-ROM | disk | FTP | other *** search
/ nisttime.carsoncity.k12.mi.us / nisttime.carsoncity.k12.mi.us.tar / nisttime.carsoncity.k12.mi.us / pub / autolock / adev.c next >
C/C++ Source or Header  |  2003-12-09  |  2KB  |  48 lines
  1. #include <math.h>
  2. double adev(double y[],int n,int lag)
  3. {
  4. /*
  5.     this subroutine computes the approximate Allan
  6.     deviation of the frequency data in array y
  7.     (as produced by autolock, for example).  There
  8.     are n points in the array, and the deviation is
  9.     computed for lag.  The computed value is returned
  10.     through the function.
  11.     this subroutine assumes that the input data are
  12.     equally spaced in time, but it has no way of 
  13.     checking that this is true.
  14.     the input array y represents the frequency averages
  15.     for a constant averaging time -- basically one unit
  16.     of lag. thus this subroutine is called with a value
  17.     of lag >= 1.
  18.     if the value of lag is > 1 then the definition of
  19.     adev requires frequency estimates that are made using
  20.     an averaging time of lag values.  if the input array
  21.     contain times, then it is a simple matter to compute
  22.     the second difference of the input array using points
  23.     spaced lag apart. This can't be done for frequency
  24.     values however since the input values are always 
  25.     computed as frequency estimates with unit lag.
  26. */
  27. double ysum;    /*running sum of freq. differences*/
  28. double y1,y2;     /*the average frequency over the first and second intervals*/
  29. double yy;    /*the difference y2 - y1*/
  30. int npt;    /*twice the number of points in ysum*/
  31. int i;
  32. int j;
  33.     if(lag > n/4) return(0);    /*too big a lag for the array*/
  34.     if(lag <=  0) return(0);    /*error */
  35.     ysum=0;
  36.     npt= 0;
  37.     for(i=0; i < (n-2*lag); i++)      /*the index i points to the 1st point*/
  38.        {
  39.        y1=y2=0;
  40.        for(j=i; j< i+lag; j++) y1 += y[j];
  41.        for(j=i+lag; j< i+2*lag; j++) y2 += y[j];
  42.        yy= (y2 - y1)/lag;
  43.        ysum += yy*yy;
  44.        npt += 2;    /*add 2 for every point - thus twice number of points*/
  45.     }
  46.     return(sqrt(ysum/npt));
  47. }
  48.