1 /* test of fft4f2d.c */
2
3 #include <math.h>
4 #include <stdio.h>
5 #include "alloc.h"
6 #define MAX(x,y) ((x) > (y) ? (x) : (y))
7
8 /* random number generator, 0 <= RND < 1 */
9 #define RND(p) ((*(p) = (*(p) * 7141 + 54773) % 259200) * (1.0 / 259200))
10
11
main()12 int main()
13 {
14 void cdft2d(int, int, int, double **, int *, double *);
15 void rdft2d(int, int, int, double **, int *, double *);
16 void ddct2d(int, int, int, double **, double **, int *, double *);
17 void ddst2d(int, int, int, double **, double **, int *, double *);
18 void putdata2d(int n1, int n2, double **a);
19 double errorcheck2d(int n1, int n2, double scale, double **a);
20 int *ip, n1, n2, n, i;
21 double **a, **t, *w, err;
22
23 printf("data length n1=? (n1 = power of 2) \n");
24 scanf("%d", &n1);
25 printf("data length n2=? (n2 = power of 2) \n");
26 scanf("%d", &n2);
27
28 a = alloc_2d_double(n1, n2);
29 t = alloc_2d_double(n1, n2);
30 n = MAX(n1, n2 / 2);
31 ip = alloc_1d_int(2 + (int) sqrt(n + 0.5));
32 n = MAX(n1 / 2, n2 / 4) + MAX(n1, n2);
33 w = alloc_1d_double(n);
34 ip[0] = 0;
35
36 /* check of CDFT */
37 putdata2d(n1, n2, a);
38 cdft2d(n1, n2, 1, a, ip, w);
39 cdft2d(n1, n2, -1, a, ip, w);
40 err = errorcheck2d(n1, n2, 2.0 / n1 / n2, a);
41 printf("cdft2d err= %g \n", err);
42
43 /* check of RDFT */
44 putdata2d(n1, n2, a);
45 rdft2d(n1, n2, 1, a, ip, w);
46 rdft2d(n1, n2, -1, a, ip, w);
47 err = errorcheck2d(n1, n2, 2.0 / n1 / n2, a);
48 printf("rdft2d err= %g \n", err);
49
50 /* check of DDCT */
51 putdata2d(n1, n2, a);
52 ddct2d(n1, n2, 1, a, t, ip, w);
53 ddct2d(n1, n2, -1, a, t, ip, w);
54 for (i = 0; i <= n1 - 1; i++) {
55 a[i][0] *= 0.5;
56 }
57 for (i = 0; i <= n2 - 1; i++) {
58 a[0][i] *= 0.5;
59 }
60 err = errorcheck2d(n1, n2, 4.0 / n1 / n2, a);
61 printf("ddct2d err= %g \n", err);
62
63 /* check of DDST */
64 putdata2d(n1, n2, a);
65 ddst2d(n1, n2, 1, a, t, ip, w);
66 ddst2d(n1, n2, -1, a, t, ip, w);
67 for (i = 0; i <= n1 - 1; i++) {
68 a[i][0] *= 0.5;
69 }
70 for (i = 0; i <= n2 - 1; i++) {
71 a[0][i] *= 0.5;
72 }
73 err = errorcheck2d(n1, n2, 4.0 / n1 / n2, a);
74 printf("ddst2d err= %g \n", err);
75
76 free_1d_double(w);
77 free_1d_int(ip);
78 free_2d_double(t);
79 free_2d_double(a);
80 return 0;
81 }
82
83
putdata2d(int n1,int n2,double ** a)84 void putdata2d(int n1, int n2, double **a)
85 {
86 int j1, j2, seed = 0;
87
88 for (j1 = 0; j1 <= n1 - 1; j1++) {
89 for (j2 = 0; j2 <= n2 - 1; j2++) {
90 a[j1][j2] = RND(&seed);
91 }
92 }
93 }
94
95
errorcheck2d(int n1,int n2,double scale,double ** a)96 double errorcheck2d(int n1, int n2, double scale, double **a)
97 {
98 int j1, j2, seed = 0;
99 double err = 0, e;
100
101 for (j1 = 0; j1 <= n1 - 1; j1++) {
102 for (j2 = 0; j2 <= n2 - 1; j2++) {
103 e = RND(&seed) - a[j1][j2] * scale;
104 err = MAX(err, fabs(e));
105 }
106 }
107 return err;
108 }
109
110