John Miles wrote:
If you'd like to make the C source available, I'll look at building an
incremental version...?
If you can accept the current quick-hack state of the code (yes, I am
not too proud of the cleanness, so no comments about coding style or the
general uglyness), find it attached. The main goal of this code is
correctness. Find it attached.
I'll hacking away further.
An incremental version is a trivial extention for MTOTDEV as the outer
loop steps the 3*m sample window over the data-samples. So you need to
keep a 3*m big backlog of samples for the largest m you need. For
simplicity I made the expantion into the 9*m series explicit and with a
3*m offset from the index of the original data only for clarity while
debugging the code. Some of that debugging was not in the core
algorithm, but in the missing 0.73 factor, which is the bias factor for
white FM noise. A good quality measure would adjust the bias factor
based on the dominant noise-source, but that comes as later refinement.
Anyway, hope you find the code usefull. Will hack away to make more
values calculated. Currently the Hadamard Total Deviation is not
correct, as a quick-check with SP1065 page 118 would tell you. But that
was a quick-hack anyway.
Oh, yes... among the particular uglyness issues is that the routines
know just how many samples the time or frequency sample series contains.
Ugly as hell, but for the intended purpose (known test-series) that is ok.
Cheers,
Magnus
#include <math.h>
#include <stdio.h>
double td_f[1000];
double td_t[1001];
double td_mtot_t[1601];
// Expected series of l is
// l0 = 1234567890
// l1 = 395529916
// l2 = 1209410747
// l3 = 633705974
void ts_gen()
{
int i;
long long l;
l = 1234567890;
td_t[0] = 0;
//printf("%12.10f\n", td_t[0]);
for (i = 0; i < 1000; i++)
{
//printf("%i %li\n", i, l);
td_f[i] = (double)l / 2147483647;
td_t[i+1] = td_t[i] + td_f[i];
//printf("%14.10f\n", td_t[i+1]);
l = (16807 * l) % 2147483647;
}
}
void ts_gen_drift(double D, double tau0)
{
int i;
td_t[0] = 0;
for (i = 0; i < 1000; i++)
{
td_f[i] = D*i*tau0;
td_t[i+1] = td_t[i] + td_f[i];
}
}
double get_avg(int i, int m)
{
int j;
double a;
a = 0;
for (j = 0; j < m; j++)
a += td_f[i+j];
return a / m;
}
double calc_max(int m)
{
int i;
double n,o;
for (i = 0; i < 1000; i+=m)
{
n = get_avg(i,m);
if (i == 0)
o = n;
else if (n > o)
o = n;
}
return o;
}
double calc_min(int m)
{
int i;
double n,o;
for (i = 0; i < 1000; i+=m)
{
n = get_avg(i,m);
if (i == 0)
o = n;
else if (n < o)
o = n;
}
return o;
}
double calc_avg(int m)
{
int i;
double a;
a = 0;
for (i = 0; i < 1000; i+=m)
a += get_avg(i,m);
return a / (1000/m);
}
double calc_sdev(int m)
{
int i;
double a,d;
a = calc_avg(m);
d = 0;
for (i = 0; i < 1000; i+=m)
d += pow(get_avg(i,m) - a, 2);
return sqrt(d / (1000/m - 1));
}
double calc_nadev(int m)
{
int i;
double d;
d = 0;
for (i = 0; i + m < 1000; i+=m)
d += pow(get_avg(i+m,m) - get_avg(i,m), 2);
return sqrt(d / (2 * (1000/m - 1)));
}
double calc_osadev(int m)
{
int j,i,c;
double d,e;
d = 0;
c = 1000;
for (i = 0; i < c -2*m+1; i++)
{
e = 0;
for (j = 0; j < m; j++)
e += td_f[i+j+m] - td_f[i+j];
d += pow(e,2);
}
return sqrt(d / (2 *m*m * (c - 2*m+1)));
}
double calc_oadev(int m)
{
int i,c;
double d,e;
d = 0;
c = 1000;
for (i = 0; i < c -2*m+1; i++)
{
e = td_t[i+2*m] - 2 * td_t[i+m] + td_t[i];
d += pow(e,2);
}
return sqrt(d / (2 *m*m * (c - 2*m+1)));
}
double calc_madev(int m)
{
int j,i,k,c;
double d,e,md;
d = 0;
c = 1000;
for (j = 1; j <= c -3*m+2; j++)
{
e = 0;
for (i = j; i <= j+m-1; i++)
for (k = i; k <= i+m-1; k++)
e += td_f[k+m-1] - td_f[k-1];
d += e*e;
}
md = m;
return sqrt(d / (2 *md*md*md*md * (c - 3*md+2)));
}
double calc_tdev(int m)
{
int j,i,k,c;
double d,e,md;
d = 0;
c = 1000;
for (j = 1; j <= c -3*m+2; j++)
{
e = 0;
for (i = j; i <= j+m-1; i++)
for (k = i; k <= i+m-1; k++)
e += td_f[k+m-1] - td_f[k-1];
d += e*e;
}
md = m;
return sqrt(d / (6 *md*md * (c - 3*md+2)));
}
double calc_hdev(int m)
{
int i;
double d;
d = 0;
for (i = 0; i + 2*m < 1000; i+=m)
d += pow(get_avg(i+2*m,m) - 2*get_avg(i+m,m) + get_avg(i,m), 2);
return sqrt(d / (6 * (1000/m - 2)));
}
double calc_ohdev(int m)
{
int i,c;
double d,e,md;
d = 0;
c = 1001;
for (i = 0; i < c -3*m; i++)
{
e = td_t[i+3*m] - 3 * td_t[i+2*m] + 3 * td_t[i+m] - td_t[i];
d += e*e;
}
md = m;
return sqrt(d / (6 *md*md * (c - 3*md)));
}
double calc_mhdev(int m)
{
int i,j,c;
double d,e,md;
d = 0;
c = 1001;
for (j = 0; j < c -4*m; j++)
{
e = 0;
for (i=j; i <= j+m-1; i++)
e -= td_t[i+3*m] - 3 * td_t[i+2*m] + 3 * td_t[i+m] - td_t[i];
d += e*e;
}
md = m;
return sqrt(d / (6 *md*md*md*md * (c - 4*md+1)));
}
double get_wt(int i, int N)
{
i = i - 1;
if (i < 0)
{
i = -i;
return 2*td_t[0] - td_t[i];
}
else if (i > N - 1)
{
i = 2*(N-1) - i;
return 2*td_t[N-1] - td_t[i];
}
return td_t[i];
}
double calc_totdev(int m)
{
int c,i;
double d,e,md;
d = 0;
c = 1001;
for (i = 2; i <= c-1; i++)
{
e = get_wt(i-m,c) -2*get_wt(i,c) + get_wt(i+m,c);
d += e*e;
}
md = m;
return sqrt(d / (2 *md*md* (c - 2)));
}
/**
* 1 <= n <= N
* 1 <= l <= 3m
* i = l - 1
* x0(n+o) = td_mtot_t(o)
* x0#(n-l) = x0(n+l-1) = td_mtot_t(l-1) = td_mtot_t(i)
* x0#(n+l) = x0(n+l) = td_mtot_t(l); = td_mtot_t(i+1)
* x0#(n+3m+l-1) = x0(n+3m-l) = td_mtot_t(3m-l) = td_mtot_t(3m-i+1)
*
* i = n-l => l = n-i
* i = n+3m+l-1 => l = i - n - 3m +1
*/
double get_mtot_wt(int i, int n, int m)
{
int l;
if (i < n+1)
{
l = n-i;
i = n+l-1;
return td_mtot_t[i-(n-1)];
}
else if (i > n + 3*m)
{
l = i - n - 3*m + 1;
i = n+3*m - l;
return td_mtot_t[i-(n-1)];
}
return td_mtot_t[i-(n-1)];
}
double get_mtot_wtavg(int i, int n, int m)
{
int j;
double a,md;
a = 0;
for (j = 0; j < m; j++)
a+=get_mtot_wt(i+j,n,m);
md = m;
return a/md;
}
double get_mtot_wtavg2(int i, int m)
{
int j;
double a,md;
a = 0;
for (j = 0; j < m; j++)
a+=td_mtot_t[i+j];
md = m;
return a/md;
}
double calc_mtotdev(int m)
{
int n,i,c,j,x;
double d,e,x1,x2,f,md;
d = 0;
c = 1001;
md = m;
for (n = 0; n < c -3*m+1; n++)
{
// Calculate average frequency
x = (3*m+1)/2;
for (x1 = 0, i = 0, j = 0; i < 3*m/2; i++, j++)
x1 += (td_t[n+i+x]-td_t[n+i])/x;
f = x1/j;
// Compensate for average frequency
for (i = 0; i <= 3*m-1; i++)
{
int l;
l = i + 1;
e = td_t[n+i] - f*i;
//td_mtot_t[3*m-1-i] = e;
//td_mtot_t[3*m +i] = e;
//td_mtot_t[9*m-1-i] = e;
td_mtot_t[300+n-l] = e;
td_mtot_t[300+n+l-1] = e;
// i = 3*m - h => h = 3*m - i
// n + 3*m + h - 1 = n + 6*m - i - 1 = n + 6*m - l
td_mtot_t[300+n-l+6*m] = e;
}
//printf("%e\n", (td_mtot_t[3*m-1]-td_mtot_t[0])/(3*m-1));
e = 0;
//for (i = 0; i <= 6*m-1; i++)
for (i = n-3*m; i <= n+3*m-1; i++)
{
//f = get_mtot_wtavg2(i,m) -2.0*get_mtot_wtavg2(i+m,m) + get_mtot_wtavg2(i+2*m,m);
f = get_mtot_wtavg2(300+i,m) -2.0*get_mtot_wtavg2(300+i+m,m) + get_mtot_wtavg2(300+i+2*m,m);
e += f*f;
}
d += e/(6*md);
}
return sqrt(d / (0.73 * 2 *md*md * (c - 3*md+1)));
}
double calc_ttotdev(int m)
{
double md;
md = m;
return calc_mtotdev(m)*sqrt(md*md/3);
}
double calc_htotdev(int m)
{
int n,i,c,j,x;
double d,e,x1,x2,f,md;
d = 0;
c = 1000;
md = m;
for (n = 0; n < c -3*m+1; n++)
{
// Calculate average frequency
x = (3*m+1)/2;
for (x1 = 0, i = 0, j = 0; i < 3*m/2; i++, j++)
x1 += (td_f[n+i+x]-td_f[n+i])/x;
f = x1/j;
//printf("m %i n %i x1 %e x2 %e x2-x1 %e f %e ", m, n, x1, x2, x2-x1, f);
// Compensate for average frequency
for (i = 0; i <= 3*m-1; i++)
{
int l;
l = i + 1;
e = td_f[n+i] - f*i;
//td_mtot_t[3*m-1-i] = e;
//td_mtot_t[3*m +i] = e;
//td_mtot_t[9*m-1-i] = e;
td_mtot_t[300+n-l] = e;
td_mtot_t[300+n+l-1] = e;
// i = 3*m - h => h = 3*m - i
// n + 3*m + h - 1 = n + 6*m - i - 1 = n + 6*m - l
td_mtot_t[300+n-l+6*m] = e;
}
//printf("%e\n", (td_mtot_t[3*m-1]-td_mtot_t[0])/(3*m-1));
e = 0;
//for (i = 0; i <= 6*m-1; i++)
for (i = n-3*m; i <= n+3*m-1; i++)
{
//f = get_mtot_wtavg2(i,m) -2.0*get_mtot_wtavg2(i+m,m) + get_mtot_wtavg2(i+2*m,m);
f = get_mtot_wtavg2(300+i,m) -2.0*get_mtot_wtavg2(300+i+m,m) + get_mtot_wtavg2(300+i+2*m,m);
e += f*f;
}
d += e/(6*md);
}
return sqrt(d / (6.0 * (c - 3*md+1)));
}
int main()
{
ts_gen();
//ts_gen_drift(1.0E-4, 1.0);
printf(" m = 1 m = 10 m = 100\n");
printf("Max %12.10f %12.10f %12.10f\n", calc_max(1), calc_max(10), calc_max(100));
printf("Min %12.10f %12.10f %12.10f\n", calc_min(1), calc_min(10), calc_min(100));
printf("Avg %12.10f %12.10f %12.10f\n", calc_avg(1), calc_avg(10), calc_avg(100));
printf("Sdev %12.10f %12.10f %12.10f\n", calc_sdev(1), calc_sdev(10), calc_sdev(100));
printf("NAdev %12.10f %12.10f %12.10f\n", calc_nadev(1), calc_nadev(10), calc_nadev(100));
printf("OSAdev %12.10f %12.10f %12.10f\n", calc_osadev(1), calc_osadev(10), calc_osadev(100));
printf("OAdev %12.10f %12.10f %12.10f\n", calc_oadev(1), calc_oadev(10), calc_oadev(100));
printf("MAdev %12.10f %12.10f %12.10f\n", calc_madev(1), calc_madev(10), calc_madev(100));
printf("Tdev %12.10f %12.10f %12.10f\n", calc_tdev(1), calc_tdev(10), calc_tdev(100));
printf("Hdev %12.10f %12.10f %12.10f\n", calc_hdev(1), calc_hdev(10), calc_hdev(100));
printf("OHdev %12.10f %12.10f %12.10f\n", calc_ohdev(1), calc_ohdev(10), calc_ohdev(100));
printf("MHdev %12.10f %12.10f %12.10f\n", calc_mhdev(1), calc_mhdev(10), calc_mhdev(100));
printf("HTOTdev %12.10f %12.10f %12.10f\n", calc_htotdev(1), calc_htotdev(10), calc_htotdev(100));
printf("TOTdev %12.10f %12.10f %12.10f\n", calc_totdev(1), calc_totdev(10), calc_totdev(100));
printf("MTOTdev %12.10f %12.10f %12.10f\n", calc_mtotdev(1), calc_mtotdev(10), calc_mtotdev(100));
printf("TTOTdev %12.10f %12.10f %12.10f\n", calc_ttotdev(1), calc_ttotdev(10), calc_ttotdev(100));
return -1;
}
_______________________________________________
time-nuts mailing list -- [email protected]
To unsubscribe, go to https://www.febo.com/cgi-bin/mailman/listinfo/time-nuts
and follow the instructions there.
