1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
|
/* Test smoothness scaling behaviour for 1D */
#include <stdio.h>
#include <math.h>
double trans(double *v, int luord, double vv);
double f(double x) {
double fp[5] = { +1.0, 0.7, -0.3, 0.0, 0.0 };
double y;
y = trans(fp, 5, x);
return y;
}
int main() {
double min = 0.0;
double max = 1.0;
int i, res = 4;
int di = 1;
#define LOC(xx) (min + (max-min) * (xx)/(res-1.0))
/* For each resolution */
for (i = 0; i < 10; i++, res *= 2) {
double tse = 0.0; /* Total squared error */
int j;
/* For each grid point with neigbors */
for (j = 1; j < (res-1); j++) {
double err;
double y1, y2, y3;
y1 = f(LOC(j-1));
y2 = f(LOC(j));
y3 = f(LOC(j+1));
err = 0.5 * (y3 + y1) - y2;
tse += err * err;
// tse += fabs(err);
}
/* Apply adjustments and corrections to error squared */
tse *= pow((res-1.0), 4.0); /* Aprox. geometric resolution factor */
tse /= pow((res-2.0),(double)di); /* Average squared non-smoothness */
// tse /= (di * pow((res-2.0),(double)di)); /* Average squared non-smoothness */
printf("Res %d, tse = %f\n",res,tse);
}
return 0;
}
/* Transfer function */
double trans(
double *v, /* Pointer to first parameter */
int luord, /* Number of parameters */
double vv /* Source of value */
) {
double g;
int ord;
for (ord = 0; ord < luord; ord++) {
int nsec; /* Number of sections */
double sec; /* Section */
g = v[ord]; /* Parameter */
nsec = ord + 1; /* Increase sections for each order */
vv *= (double)nsec;
sec = floor(vv);
if (((int)sec) & 1)
g = -g; /* Alternate action in each section */
vv -= sec;
if (g >= 0.0) {
vv = vv/(g - g * vv + 1.0);
} else {
vv = (vv - g * vv)/(1.0 - g * vv);
}
vv += sec;
vv /= (double)nsec;
}
return vv;
}
|
