1 #include <stdlib.h>
2 #include <stdio.h>
3
4 #include <fftw3.h>
5
6 #if DCT_TEST_PRECISION == 1
7 typedef float float_prec;
8 #define PF "%.7f"
9 #define FFTW_PLAN fftwf_plan
10 #define FFTW_MALLOC fftwf_malloc
11 #define FFTW_FREE fftwf_free
12 #define FFTW_PLAN_CREATE fftwf_plan_r2r_1d
13 #define FFTW_EXECUTE fftwf_execute
14 #define FFTW_DESTROY_PLAN fftwf_destroy_plan
15 #define FFTW_CLEANUP fftwf_cleanup
16 #elif DCT_TEST_PRECISION == 2
17 typedef double float_prec;
18 #define PF "%.18f"
19 #define FFTW_PLAN fftw_plan
20 #define FFTW_MALLOC fftw_malloc
21 #define FFTW_FREE fftw_free
22 #define FFTW_PLAN_CREATE fftw_plan_r2r_1d
23 #define FFTW_EXECUTE fftw_execute
24 #define FFTW_DESTROY_PLAN fftw_destroy_plan
25 #define FFTW_CLEANUP fftw_cleanup
26 #elif DCT_TEST_PRECISION == 3
27 typedef long double float_prec;
28 #define PF "%.18Lf"
29 #define FFTW_PLAN fftwl_plan
30 #define FFTW_MALLOC fftwl_malloc
31 #define FFTW_FREE fftwl_free
32 #define FFTW_PLAN_CREATE fftwl_plan_r2r_1d
33 #define FFTW_EXECUTE fftwl_execute
34 #define FFTW_DESTROY_PLAN fftwl_destroy_plan
35 #define FFTW_CLEANUP fftwl_cleanup
36 #else
37 #error DCT_TEST_PRECISION must be a number 1-3
38 #endif
39
40
41 enum type {
42 DCT_I = 1,
43 DCT_II = 2,
44 DCT_III = 3,
45 DCT_IV = 4,
46 DST_I = 5,
47 DST_II = 6,
48 DST_III = 7,
49 DST_IV = 8,
50 };
51
gen(int type,int sz)52 int gen(int type, int sz)
53 {
54 float_prec *a, *b;
55 FFTW_PLAN p;
56 int i, tp;
57
58 a = FFTW_MALLOC(sizeof(*a) * sz);
59 if (a == NULL) {
60 fprintf(stderr, "failure\n");
61 exit(EXIT_FAILURE);
62 }
63 b = FFTW_MALLOC(sizeof(*b) * sz);
64 if (b == NULL) {
65 fprintf(stderr, "failure\n");
66 exit(EXIT_FAILURE);
67 }
68
69 switch(type) {
70 case DCT_I:
71 tp = FFTW_REDFT00;
72 break;
73 case DCT_II:
74 tp = FFTW_REDFT10;
75 break;
76 case DCT_III:
77 tp = FFTW_REDFT01;
78 break;
79 case DCT_IV:
80 tp = FFTW_REDFT11;
81 break;
82 case DST_I:
83 tp = FFTW_RODFT00;
84 break;
85 case DST_II:
86 tp = FFTW_RODFT10;
87 break;
88 case DST_III:
89 tp = FFTW_RODFT01;
90 break;
91 case DST_IV:
92 tp = FFTW_RODFT11;
93 break;
94 default:
95 fprintf(stderr, "unknown type\n");
96 exit(EXIT_FAILURE);
97 }
98
99 switch(type) {
100 case DCT_I:
101 case DCT_II:
102 case DCT_III:
103 case DCT_IV:
104 for(i=0; i < sz; ++i) {
105 a[i] = i;
106 }
107 break;
108 case DST_I:
109 case DST_II:
110 case DST_III:
111 case DST_IV:
112 /* TODO: what should we do for dst's?*/
113 for(i=0; i < sz; ++i) {
114 a[i] = i;
115 }
116 break;
117 default:
118 fprintf(stderr, "unknown type\n");
119 exit(EXIT_FAILURE);
120 }
121
122 p = FFTW_PLAN_CREATE(sz, a, b, tp, FFTW_ESTIMATE);
123 FFTW_EXECUTE(p);
124 FFTW_DESTROY_PLAN(p);
125
126 for(i=0; i < sz; ++i) {
127 printf(PF"\n", b[i]);
128 }
129 FFTW_FREE(b);
130 FFTW_FREE(a);
131
132 return 0;
133 }
134
main(int argc,char * argv[])135 int main(int argc, char* argv[])
136 {
137 int n, tp;
138
139 if (argc < 3) {
140 fprintf(stderr, "missing argument: program type n\n");
141 exit(EXIT_FAILURE);
142 }
143 tp = atoi(argv[1]);
144 n = atoi(argv[2]);
145
146 gen(tp, n);
147 FFTW_CLEANUP();
148
149 return 0;
150 }
151