-
Notifications
You must be signed in to change notification settings - Fork 2
/
test_1d_rodft10.c
111 lines (89 loc) · 3.02 KB
/
test_1d_rodft10.c
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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <fftw3.h>
#include "util.h"
int test_1d_rodft10(int n) {
int N = 2 * n;
double *in = fftw_alloc_real(n);
fftw_complex *in_logical = fftw_alloc_complex(N);
double *out = fftw_alloc_real(n);
fftw_complex *out_logical = fftw_alloc_complex(N);
fftw_plan p = fftw_plan_r2r_1d(n, in, out, FFTW_RODFT10, FFTW_ESTIMATE);
fill_random_1d_real(n, in);
// the first half of the array is identical
for (int i = 0; i < n; ++i) {
in_logical[i] = in[i];
}
// second half is filled according to odd symmetry around (n-0.5)
for (int i = 0; i < n; ++i) {
in_logical[n + i] = -in[n - 1 - i];
}
dump_1d_real("test_1d_rodft10_in.dat", n, in);
dump_1d_cplx("test_1d_rodft10_in_logical.dat", N, in_logical);
// manual implementation of logically-equivalent DFT
double a = 0.5;
double b = 0.0;
dft_1d_cplx(N, in_logical, out_logical, a, b);
fftw_execute(p);
dump_1d_real("test_1d_rodft10_out.dat", n, out);
dump_1d_cplx("test_1d_rodft10_out_logical.dat", N, out_logical);
// check output of logically equivalent DFT
double eps = 1e-12;
int status = 0;
// 1. logically equivalent output should be purely imaginary-valued
for (int i = 0; i < N; ++i) {
if (fabs(creal(out_logical[i])) > eps) {
printf("error: real of [%d] is %g\n", i, creal(out_logical[i]));
status = 1;
} else {
printf("real of [%d] is %g\n", i, creal(out_logical[i]));
}
}
// 2. first output value is zero (?)
if (fabs(cimag(out_logical[0])) > eps) {
printf("error: delta of [%d] is %g\n", 0, -cimag(out_logical[0]));
status = 1;
} else {
printf("match of [%d] (delta=%g)\n", 0, -cimag(out_logical[0]));
}
// 3. next n values should have the output of RODFT10 in the negative imaginary part
double delta;
for (int i = 0; i < n; ++i) {
delta = -cimag(out_logical[i + 1]) - out[i];
if (fabs(delta) > eps) {
printf("error: delta of [%d] is %g\n", i + 1, delta);
status = 1;
} else {
printf("match of [%d] (delta=%g)\n", i + 1, delta);
}
}
// 4. even symmetry around (n-1) of neg. imag. part vs. output of RODFT10
for (int i = 0; i < n - 1; ++i) {
delta = -cimag(out_logical[n + 1 + i]) - out[n - 2 - i];
if (fabs(delta) > eps) {
printf("error: delta of [%d] is %g\n", n + 1 + i, delta);
status = 1;
} else {
printf("match of [%d] (delta=%g)\n", n + 1 + i, delta);
}
}
if (status == 0) {
printf("=> all ok\n");
} else {
printf("=> errors\n");
}
fftw_destroy_plan(p);
fftw_free(in);
fftw_free(in_logical);
fftw_free(out);
fftw_free(out_logical);
return status;
}
int main(int argc, char** argv) {
int status = 0;
status += test_1d_rodft10(4);
status += test_1d_rodft10(5);
return status;
}