1 /*
2 Fast Fourier/Cosine/Sine Transform
3 dimension :two
4 data length :power of 2
5 decimation :frequency
6 radix :4, 2, row-column
7 data :inplace
8 table :use
9 functions
10 cdft2d: Complex Discrete Fourier Transform
11 rdft2d: Real Discrete Fourier Transform
12 ddct2d: Discrete Cosine Transform
13 ddst2d: Discrete Sine Transform
14 function prototypes
15 void cdft2d(int, int, int, double **, int *, double *);
16 void rdft2d(int, int, int, double **, int *, double *);
17 void ddct2d(int, int, int, double **, double **, int *, double *);
18 void ddst2d(int, int, int, double **, double **, int *, double *);
19
20
21 -------- Complex DFT (Discrete Fourier Transform) --------
22 [definition]
23 <case1>
24 X[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 x[j1][j2] *
25 exp(2*pi*i*j1*k1/n1) *
26 exp(2*pi*i*j2*k2/n2), 0<=k1<n1, 0<=k2<n2
27 <case2>
28 X[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 x[j1][j2] *
29 exp(-2*pi*i*j1*k1/n1) *
30 exp(-2*pi*i*j2*k2/n2), 0<=k1<n1, 0<=k2<n2
31 (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
32 [usage]
33 <case1>
34 ip[0] = 0; // first time only
35 cdft2d(n1, 2*n2, 1, a, ip, w);
36 <case2>
37 ip[0] = 0; // first time only
38 cdft2d(n1, 2*n2, -1, a, ip, w);
39 [parameters]
40 n1 :data length (int)
41 n1 >= 1, n1 = power of 2
42 2*n2 :data length (int)
43 n2 >= 1, n2 = power of 2
44 a[0...n1-1][0...2*n2-1]
45 :input/output data (double **)
46 input data
47 a[j1][2*j2] = Re(x[j1][j2]),
48 a[j1][2*j2+1] = Im(x[j1][j2]),
49 0<=j1<n1, 0<=j2<n2
50 output data
51 a[k1][2*k2] = Re(X[k1][k2]),
52 a[k1][2*k2+1] = Im(X[k1][k2]),
53 0<=k1<n1, 0<=k2<n2
54 ip[0...*]
55 :work area for bit reversal (int *)
56 length of ip >= 2+sqrt(n)
57 (n = max(n1, n2))
58 ip[0],ip[1] are pointers of the cos/sin table.
59 w[0...*]
60 :cos/sin table (double *)
61 length of w >= max(n1/2, n2/2)
62 w[],ip[] are initialized if ip[0] == 0.
63 [remark]
64 Inverse of
65 cdft2d(n1, 2*n2, -1, a, ip, w);
66 is
67 cdft2d(n1, 2*n2, 1, a, ip, w);
68 for (j1 = 0; j1 <= n1 - 1; j1++) {
69 for (j2 = 0; j2 <= 2 * n2 - 1; j2++) {
70 a[j1][j2] *= 1.0 / (n1 * n2);
71 }
72 }
73 .
74
75
76 -------- Real DFT / Inverse of Real DFT --------
77 [definition]
78 <case1> RDFT
79 R[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 a[j1][j2] *
80 cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2),
81 0<=k1<n1, 0<=k2<n2
82 I[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 a[j1][j2] *
83 sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2),
84 0<=k1<n1, 0<=k2<n2
85 <case2> IRDFT (excluding scale)
86 a[k1][k2] = (1/2) * sum_j1=0^n1-1 sum_j2=0^n2-1
87 (R[j1][j2] *
88 cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2) +
89 I[j1][j2] *
90 sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2)),
91 0<=k1<n1, 0<=k2<n2
92 (notes: R[n1-k1][n2-k2] = R[k1][k2],
93 I[n1-k1][n2-k2] = -I[k1][k2],
94 R[n1-k1][0] = R[k1][0],
95 I[n1-k1][0] = -I[k1][0],
96 R[0][n2-k2] = R[0][k2],
97 I[0][n2-k2] = -I[0][k2],
98 0<k1<n1, 0<k2<n2)
99 [usage]
100 <case1>
101 ip[0] = 0; // first time only
102 rdft2d(n1, n2, 1, a, ip, w);
103 <case2>
104 ip[0] = 0; // first time only
105 rdft2d(n1, n2, -1, a, ip, w);
106 [parameters]
107 n1 :data length (int)
108 n1 >= 2, n1 = power of 2
109 n2 :data length (int)
110 n2 >= 2, n2 = power of 2
111 a[0...n1-1][0...n2-1]
112 :input/output data (double **)
113 <case1>
114 output data
115 a[k1][2*k2] = R[k1][k2] = R[n1-k1][n2-k2],
116 a[k1][2*k2+1] = I[k1][k2] = -I[n1-k1][n2-k2],
117 0<k1<n1, 0<k2<n2/2,
118 a[0][2*k2] = R[0][k2] = R[0][n2-k2],
119 a[0][2*k2+1] = I[0][k2] = -I[0][n2-k2],
120 0<k2<n2/2,
121 a[k1][0] = R[k1][0] = R[n1-k1][0],
122 a[k1][1] = I[k1][0] = -I[n1-k1][0],
123 a[n1-k1][1] = R[k1][n2/2] = R[n1-k1][n2/2],
124 a[n1-k1][0] = -I[k1][n2/2] = I[n1-k1][n2/2],
125 0<k1<n1/2,
126 a[0][0] = R[0][0],
127 a[0][1] = R[0][n2/2],
128 a[n1/2][0] = R[n1/2][0],
129 a[n1/2][1] = R[n1/2][n2/2]
130 <case2>
131 input data
132 a[j1][2*j2] = R[j1][j2] = R[n1-j1][n2-j2],
133 a[j1][2*j2+1] = I[j1][j2] = -I[n1-j1][n2-j2],
134 0<j1<n1, 0<j2<n2/2,
135 a[0][2*j2] = R[0][j2] = R[0][n2-j2],
136 a[0][2*j2+1] = I[0][j2] = -I[0][n2-j2],
137 0<j2<n2/2,
138 a[j1][0] = R[j1][0] = R[n1-j1][0],
139 a[j1][1] = I[j1][0] = -I[n1-j1][0],
140 a[n1-j1][1] = R[j1][n2/2] = R[n1-j1][n2/2],
141 a[n1-j1][0] = -I[j1][n2/2] = I[n1-j1][n2/2],
142 0<j1<n1/2,
143 a[0][0] = R[0][0],
144 a[0][1] = R[0][n2/2],
145 a[n1/2][0] = R[n1/2][0],
146 a[n1/2][1] = R[n1/2][n2/2]
147 ip[0...*]
148 :work area for bit reversal (int *)
149 length of ip >= 2+sqrt(n)
150 (n = max(n1, n2/2))
151 ip[0],ip[1] are pointers of the cos/sin table.
152 w[0...*]
153 :cos/sin table (double *)
154 length of w >= max(n1/2, n2/4) + n2/4
155 w[],ip[] are initialized if ip[0] == 0.
156 [remark]
157 Inverse of
158 rdft2d(n1, n2, 1, a, ip, w);
159 is
160 rdft2d(n1, n2, -1, a, ip, w);
161 for (j1 = 0; j1 <= n1 - 1; j1++) {
162 for (j2 = 0; j2 <= n2 - 1; j2++) {
163 a[j1][j2] *= 2.0 / (n1 * n2);
164 }
165 }
166 .
167
168
169 -------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
170 [definition]
171 <case1> IDCT (excluding scale)
172 C[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 a[j1][j2] *
173 cos(pi*j1*(k1+1/2)/n1) *
174 cos(pi*j2*(k2+1/2)/n2),
175 0<=k1<n1, 0<=k2<n2
176 <case2> DCT
177 C[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 a[j1][j2] *
178 cos(pi*(j1+1/2)*k1/n1) *
179 cos(pi*(j2+1/2)*k2/n2),
180 0<=k1<n1, 0<=k2<n2
181 [usage]
182 <case1>
183 ip[0] = 0; // first time only
184 ddct2d(n1, n2, 1, a, t, ip, w);
185 <case2>
186 ip[0] = 0; // first time only
187 ddct2d(n1, n2, -1, a, t, ip, w);
188 [parameters]
189 n1 :data length (int)
190 n1 >= 2, n1 = power of 2
191 n2 :data length (int)
192 n2 >= 2, n2 = power of 2
193 a[0...n1-1][0...n2-1]
194 :input/output data (double **)
195 output data
196 a[k1][k2] = C[k1][k2], 0<=k1<n1, 0<=k2<n2
197 t[0...n1-1][0...n2-1]
198 :work area (double **)
199 ip[0...*]
200 :work area for bit reversal (int *)
201 length of ip >= 2+sqrt(n)
202 (n = max(n1, n2/2))
203 ip[0],ip[1] are pointers of the cos/sin table.
204 w[0...*]
205 :cos/sin table (double *)
206 length of w >= max(n1/2, n2/4) + max(n1, n2)
207 w[],ip[] are initialized if ip[0] == 0.
208 [remark]
209 Inverse of
210 ddct2d(n1, n2, -1, a, t, ip, w);
211 is
212 for (j1 = 0; j1 <= n1 - 1; j1++) {
213 a[j1][0] *= 0.5;
214 }
215 for (j2 = 0; j2 <= n2 - 1; j2++) {
216 a[0][j2] *= 0.5;
217 }
218 ddct2d(n1, n2, 1, a, t, ip, w);
219 for (j1 = 0; j1 <= n1 - 1; j1++) {
220 for (j2 = 0; j2 <= n2 - 1; j2++) {
221 a[j1][j2] *= 4.0 / (n1 * n2);
222 }
223 }
224 .
225
226
227 -------- DST (Discrete Sine Transform) / Inverse of DST --------
228 [definition]
229 <case1> IDST (excluding scale)
230 S[k1][k2] = sum_j1=1^n1 sum_j2=1^n2 A[j1][j2] *
231 sin(pi*j1*(k1+1/2)/n1) *
232 sin(pi*j2*(k2+1/2)/n2),
233 0<=k1<n1, 0<=k2<n2
234 <case2> DST
235 S[k1][k2] = sum_j1=0^n1-1 sum_j2=0^n2-1 a[j1][j2] *
236 sin(pi*(j1+1/2)*k1/n1) *
237 sin(pi*(j2+1/2)*k2/n2),
238 0<k1<=n1, 0<k2<=n2
239 [usage]
240 <case1>
241 ip[0] = 0; // first time only
242 ddst2d(n1, n2, 1, a, t, ip, w);
243 <case2>
244 ip[0] = 0; // first time only
245 ddst2d(n1, n2, -1, a, t, ip, w);
246 [parameters]
247 n1 :data length (int)
248 n1 >= 2, n1 = power of 2
249 n2 :data length (int)
250 n2 >= 2, n2 = power of 2
251 a[0...n1-1][0...n2-1]
252 :input/output data (double **)
253 <case1>
254 input data
255 a[j1][j2] = A[j1][j2], 0<j1<n1, 0<j2<n2,
256 a[j1][0] = A[j1][n2], 0<j1<n1,
257 a[0][j2] = A[n1][j2], 0<j2<n2,
258 a[0][0] = A[n1][n2]
259 (i.e. A[j1][j2] = a[j1 % n1][j2 % n2])
260 output data
261 a[k1][k2] = S[k1][k2], 0<=k1<n1, 0<=k2<n2
262 <case2>
263 output data
264 a[k1][k2] = S[k1][k2], 0<k1<n1, 0<k2<n2,
265 a[k1][0] = S[k1][n2], 0<k1<n1,
266 a[0][k2] = S[n1][k2], 0<k2<n2,
267 a[0][0] = S[n1][n2]
268 (i.e. S[k1][k2] = a[k1 % n1][k2 % n2])
269 t[0...n1-1][0...n2-1]
270 :work area (double **)
271 ip[0...*]
272 :work area for bit reversal (int *)
273 length of ip >= 2+sqrt(n)
274 (n = max(n1, n2/2))
275 ip[0],ip[1] are pointers of the cos/sin table.
276 w[0...*]
277 :cos/sin table (double *)
278 length of w >= max(n1/2, n2/4) + max(n1, n2)
279 w[],ip[] are initialized if ip[0] == 0.
280 [remark]
281 Inverse of
282 ddst2d(n1, n2, -1, a, t, ip, w);
283 is
284 for (j1 = 0; j1 <= n1 - 1; j1++) {
285 a[j1][0] *= 0.5;
286 }
287 for (j2 = 0; j2 <= n2 - 1; j2++) {
288 a[0][j2] *= 0.5;
289 }
290 ddst2d(n1, n2, 1, a, t, ip, w);
291 for (j1 = 0; j1 <= n1 - 1; j1++) {
292 for (j2 = 0; j2 <= n2 - 1; j2++) {
293 a[j1][j2] *= 4.0 / (n1 * n2);
294 }
295 }
296 .
297 */
298
299
cdft2d(int n1,int n2,int isgn,double ** a,int * ip,double * w)300 void cdft2d(int n1, int n2, int isgn, double **a, int *ip, double *w)
301 {
302 void makewt(int nw, int *ip, double *w);
303 void bitrv2col(int n1, int n, int *ip, double **a);
304 void bitrv2row(int n, int n2, int *ip, double **a);
305 void cftbcol(int n1, int n, double **a, double *w);
306 void cftbrow(int n, int n2, double **a, double *w);
307 void cftfcol(int n1, int n, double **a, double *w);
308 void cftfrow(int n, int n2, double **a, double *w);
309 int n;
310
311 n = n1 << 1;
312 if (n < n2) {
313 n = n2;
314 }
315 if (n > (ip[0] << 2)) {
316 makewt(n >> 2, ip, w);
317 }
318 if (n2 > 4) {
319 bitrv2col(n1, n2, ip + 2, a);
320 }
321 if (n1 > 2) {
322 bitrv2row(n1, n2, ip + 2, a);
323 }
324 if (isgn < 0) {
325 cftfcol(n1, n2, a, w);
326 cftfrow(n1, n2, a, w);
327 } else {
328 cftbcol(n1, n2, a, w);
329 cftbrow(n1, n2, a, w);
330 }
331 }
332
333
rdft2d(int n1,int n2,int isgn,double ** a,int * ip,double * w)334 void rdft2d(int n1, int n2, int isgn, double **a, int *ip, double *w)
335 {
336 void makewt(int nw, int *ip, double *w);
337 void makect(int nc, int *ip, double *c);
338 void bitrv2col(int n1, int n, int *ip, double **a);
339 void bitrv2row(int n, int n2, int *ip, double **a);
340 void cftbcol(int n1, int n, double **a, double *w);
341 void cftbrow(int n, int n2, double **a, double *w);
342 void cftfcol(int n1, int n, double **a, double *w);
343 void cftfrow(int n, int n2, double **a, double *w);
344 void rftbcol(int n1, int n, double **a, int nc, double *c);
345 void rftfcol(int n1, int n, double **a, int nc, double *c);
346 int n, nw, nc, n1h, i, j;
347 double xi;
348
349 n = n1 << 1;
350 if (n < n2) {
351 n = n2;
352 }
353 nw = ip[0];
354 if (n > (nw << 2)) {
355 nw = n >> 2;
356 makewt(nw, ip, w);
357 }
358 nc = ip[1];
359 if (n2 > (nc << 2)) {
360 nc = n2 >> 2;
361 makect(nc, ip, w + nw);
362 }
363 n1h = n1 >> 1;
364 if (isgn < 0) {
365 for (i = 1; i <= n1h - 1; i++) {
366 j = n1 - i;
367 xi = a[i][0] - a[j][0];
368 a[i][0] += a[j][0];
369 a[j][0] = xi;
370 xi = a[j][1] - a[i][1];
371 a[i][1] += a[j][1];
372 a[j][1] = xi;
373 }
374 if (n1 > 2) {
375 bitrv2row(n1, n2, ip + 2, a);
376 }
377 cftfrow(n1, n2, a, w);
378 for (i = 0; i <= n1 - 1; i++) {
379 a[i][1] = 0.5 * (a[i][0] - a[i][1]);
380 a[i][0] -= a[i][1];
381 }
382 if (n2 > 4) {
383 rftfcol(n1, n2, a, nc, w + nw);
384 bitrv2col(n1, n2, ip + 2, a);
385 }
386 cftfcol(n1, n2, a, w);
387 } else {
388 if (n2 > 4) {
389 bitrv2col(n1, n2, ip + 2, a);
390 }
391 cftbcol(n1, n2, a, w);
392 if (n2 > 4) {
393 rftbcol(n1, n2, a, nc, w + nw);
394 }
395 for (i = 0; i <= n1 - 1; i++) {
396 xi = a[i][0] - a[i][1];
397 a[i][0] += a[i][1];
398 a[i][1] = xi;
399 }
400 if (n1 > 2) {
401 bitrv2row(n1, n2, ip + 2, a);
402 }
403 cftbrow(n1, n2, a, w);
404 for (i = 1; i <= n1h - 1; i++) {
405 j = n1 - i;
406 a[j][0] = 0.5 * (a[i][0] - a[j][0]);
407 a[i][0] -= a[j][0];
408 a[j][1] = 0.5 * (a[i][1] + a[j][1]);
409 a[i][1] -= a[j][1];
410 }
411 }
412 }
413
414
ddct2d(int n1,int n2,int isgn,double ** a,double ** t,int * ip,double * w)415 void ddct2d(int n1, int n2, int isgn, double **a, double **t,
416 int *ip, double *w)
417 {
418 void makewt(int nw, int *ip, double *w);
419 void makect(int nc, int *ip, double *c);
420 void bitrv2col(int n1, int n, int *ip, double **a);
421 void bitrv2row(int n, int n2, int *ip, double **a);
422 void cftbcol(int n1, int n, double **a, double *w);
423 void cftbrow(int n, int n2, double **a, double *w);
424 void cftfcol(int n1, int n, double **a, double *w);
425 void cftfrow(int n, int n2, double **a, double *w);
426 void rftbcol(int n1, int n, double **a, int nc, double *c);
427 void rftfcol(int n1, int n, double **a, int nc, double *c);
428 void dctbsub(int n1, int n2, double **a, int nc, double *c);
429 void dctfsub(int n1, int n2, double **a, int nc, double *c);
430 int n, nw, nc, n1h, n2h, i, ix, ic, j, jx, jc;
431 double xi;
432
433 n = n1 << 1;
434 if (n < n2) {
435 n = n2;
436 }
437 nw = ip[0];
438 if (n > (nw << 2)) {
439 nw = n >> 2;
440 makewt(nw, ip, w);
441 }
442 nc = ip[1];
443 if (n1 > nc || n2 > nc) {
444 if (n1 > n2) {
445 nc = n1;
446 } else {
447 nc = n2;
448 }
449 makect(nc, ip, w + nw);
450 }
451 n1h = n1 >> 1;
452 n2h = n2 >> 1;
453 if (isgn >= 0) {
454 for (i = 0; i <= n1 - 1; i++) {
455 for (j = 1; j <= n2h - 1; j++) {
456 jx = j << 1;
457 t[i][jx] = a[i][j];
458 t[i][jx + 1] = a[i][n2 - j];
459 }
460 }
461 t[0][0] = a[0][0];
462 t[0][1] = a[0][n2h];
463 t[n1h][0] = a[n1h][0];
464 t[n1h][1] = a[n1h][n2h];
465 for (i = 1; i <= n1h - 1; i++) {
466 ic = n1 - i;
467 t[i][0] = a[i][0];
468 t[ic][1] = a[i][n2h];
469 t[i][1] = a[ic][0];
470 t[ic][0] = a[ic][n2h];
471 }
472 dctfsub(n1, n2, t, nc, w + nw);
473 if (n1 > 2) {
474 bitrv2row(n1, n2, ip + 2, t);
475 }
476 cftfrow(n1, n2, t, w);
477 for (i = 0; i <= n1 - 1; i++) {
478 t[i][1] = 0.5 * (t[i][0] - t[i][1]);
479 t[i][0] -= t[i][1];
480 }
481 if (n2 > 4) {
482 rftfcol(n1, n2, t, nc, w + nw);
483 bitrv2col(n1, n2, ip + 2, t);
484 }
485 cftfcol(n1, n2, t, w);
486 for (i = 0; i <= n1h - 1; i++) {
487 ix = i << 1;
488 ic = n1 - 1 - i;
489 for (j = 0; j <= n2h - 1; j++) {
490 jx = j << 1;
491 jc = n2 - 1 - j;
492 a[ix][jx] = t[i][j];
493 a[ix][jx + 1] = t[i][jc];
494 a[ix + 1][jx] = t[ic][j];
495 a[ix + 1][jx + 1] = t[ic][jc];
496 }
497 }
498 } else {
499 for (i = 0; i <= n1h - 1; i++) {
500 ix = i << 1;
501 ic = n1 - 1 - i;
502 for (j = 0; j <= n2h - 1; j++) {
503 jx = j << 1;
504 jc = n2 - 1 - j;
505 t[i][j] = a[ix][jx];
506 t[i][jc] = a[ix][jx + 1];
507 t[ic][j] = a[ix + 1][jx];
508 t[ic][jc] = a[ix + 1][jx + 1];
509 }
510 }
511 if (n2 > 4) {
512 bitrv2col(n1, n2, ip + 2, t);
513 }
514 cftbcol(n1, n2, t, w);
515 if (n2 > 4) {
516 rftbcol(n1, n2, t, nc, w + nw);
517 }
518 for (i = 0; i <= n1 - 1; i++) {
519 xi = t[i][0] - t[i][1];
520 t[i][0] += t[i][1];
521 t[i][1] = xi;
522 }
523 if (n1 > 2) {
524 bitrv2row(n1, n2, ip + 2, t);
525 }
526 cftbrow(n1, n2, t, w);
527 dctbsub(n1, n2, t, nc, w + nw);
528 for (i = 0; i <= n1 - 1; i++) {
529 for (j = 1; j <= n2h - 1; j++) {
530 jx = j << 1;
531 a[i][j] = t[i][jx];
532 a[i][n2 - j] = t[i][jx + 1];
533 }
534 }
535 a[0][0] = t[0][0];
536 a[0][n2h] = t[0][1];
537 a[n1h][0] = t[n1h][0];
538 a[n1h][n2h] = t[n1h][1];
539 for (i = 1; i <= n1h - 1; i++) {
540 ic = n1 - i;
541 a[i][0] = t[i][0];
542 a[i][n2h] = t[ic][1];
543 a[ic][0] = t[i][1];
544 a[ic][n2h] = t[ic][0];
545 }
546 }
547 }
548
549
ddst2d(int n1,int n2,int isgn,double ** a,double ** t,int * ip,double * w)550 void ddst2d(int n1, int n2, int isgn, double **a, double **t,
551 int *ip, double *w)
552 {
553 void makewt(int nw, int *ip, double *w);
554 void makect(int nc, int *ip, double *c);
555 void bitrv2col(int n1, int n, int *ip, double **a);
556 void bitrv2row(int n, int n2, int *ip, double **a);
557 void cftbcol(int n1, int n, double **a, double *w);
558 void cftbrow(int n, int n2, double **a, double *w);
559 void cftfcol(int n1, int n, double **a, double *w);
560 void cftfrow(int n, int n2, double **a, double *w);
561 void rftbcol(int n1, int n, double **a, int nc, double *c);
562 void rftfcol(int n1, int n, double **a, int nc, double *c);
563 void dstbsub(int n1, int n2, double **a, int nc, double *c);
564 void dstfsub(int n1, int n2, double **a, int nc, double *c);
565 int n, nw, nc, n1h, n2h, i, ix, ic, j, jx, jc;
566 double xi;
567
568 n = n1 << 1;
569 if (n < n2) {
570 n = n2;
571 }
572 nw = ip[0];
573 if (n > (nw << 2)) {
574 nw = n >> 2;
575 makewt(nw, ip, w);
576 }
577 nc = ip[1];
578 if (n1 > nc || n2 > nc) {
579 if (n1 > n2) {
580 nc = n1;
581 } else {
582 nc = n2;
583 }
584 makect(nc, ip, w + nw);
585 }
586 n1h = n1 >> 1;
587 n2h = n2 >> 1;
588 if (isgn >= 0) {
589 for (i = 0; i <= n1 - 1; i++) {
590 for (j = 1; j <= n2h - 1; j++) {
591 jx = j << 1;
592 t[i][jx] = a[i][j];
593 t[i][jx + 1] = a[i][n2 - j];
594 }
595 }
596 t[0][0] = a[0][0];
597 t[0][1] = a[0][n2h];
598 t[n1h][0] = a[n1h][0];
599 t[n1h][1] = a[n1h][n2h];
600 for (i = 1; i <= n1h - 1; i++) {
601 ic = n1 - i;
602 t[i][0] = a[i][0];
603 t[ic][1] = a[i][n2h];
604 t[i][1] = a[ic][0];
605 t[ic][0] = a[ic][n2h];
606 }
607 dstfsub(n1, n2, t, nc, w + nw);
608 if (n1 > 2) {
609 bitrv2row(n1, n2, ip + 2, t);
610 }
611 cftfrow(n1, n2, t, w);
612 for (i = 0; i <= n1 - 1; i++) {
613 t[i][1] = 0.5 * (t[i][0] - t[i][1]);
614 t[i][0] -= t[i][1];
615 }
616 if (n2 > 4) {
617 rftfcol(n1, n2, t, nc, w + nw);
618 bitrv2col(n1, n2, ip + 2, t);
619 }
620 cftfcol(n1, n2, t, w);
621 for (i = 0; i <= n1h - 1; i++) {
622 ix = i << 1;
623 ic = n1 - 1 - i;
624 for (j = 0; j <= n2h - 1; j++) {
625 jx = j << 1;
626 jc = n2 - 1 - j;
627 a[ix][jx] = t[i][j];
628 a[ix][jx + 1] = -t[i][jc];
629 a[ix + 1][jx] = -t[ic][j];
630 a[ix + 1][jx + 1] = t[ic][jc];
631 }
632 }
633 } else {
634 for (i = 0; i <= n1h - 1; i++) {
635 ix = i << 1;
636 ic = n1 - 1 - i;
637 for (j = 0; j <= n2h - 1; j++) {
638 jx = j << 1;
639 jc = n2 - 1 - j;
640 t[i][j] = a[ix][jx];
641 t[i][jc] = -a[ix][jx + 1];
642 t[ic][j] = -a[ix + 1][jx];
643 t[ic][jc] = a[ix + 1][jx + 1];
644 }
645 }
646 if (n2 > 4) {
647 bitrv2col(n1, n2, ip + 2, t);
648 }
649 cftbcol(n1, n2, t, w);
650 if (n2 > 4) {
651 rftbcol(n1, n2, t, nc, w + nw);
652 }
653 for (i = 0; i <= n1 - 1; i++) {
654 xi = t[i][0] - t[i][1];
655 t[i][0] += t[i][1];
656 t[i][1] = xi;
657 }
658 if (n1 > 2) {
659 bitrv2row(n1, n2, ip + 2, t);
660 }
661 cftbrow(n1, n2, t, w);
662 dstbsub(n1, n2, t, nc, w + nw);
663 for (i = 0; i <= n1 - 1; i++) {
664 for (j = 1; j <= n2h - 1; j++) {
665 jx = j << 1;
666 a[i][j] = t[i][jx];
667 a[i][n2 - j] = t[i][jx + 1];
668 }
669 }
670 a[0][0] = t[0][0];
671 a[0][n2h] = t[0][1];
672 a[n1h][0] = t[n1h][0];
673 a[n1h][n2h] = t[n1h][1];
674 for (i = 1; i <= n1h - 1; i++) {
675 ic = n1 - i;
676 a[i][0] = t[i][0];
677 a[i][n2h] = t[ic][1];
678 a[ic][0] = t[i][1];
679 a[ic][n2h] = t[ic][0];
680 }
681 }
682 }
683
684
685 /* -------- initializing routines -------- */
686
687
688 #include <math.h>
689
makewt(int nw,int * ip,double * w)690 void makewt(int nw, int *ip, double *w)
691 {
692 void bitrv2(int n, int *ip, double *a);
693 int nwh, j;
694 double delta, x, y;
695
696 ip[0] = nw;
697 ip[1] = 1;
698 if (nw > 2) {
699 nwh = nw >> 1;
700 delta = atan(1.0) / nwh;
701 w[0] = 1;
702 w[1] = 0;
703 w[nwh] = cos(delta * nwh);
704 w[nwh + 1] = w[nwh];
705 for (j = 2; j <= nwh - 2; j += 2) {
706 x = cos(delta * j);
707 y = sin(delta * j);
708 w[j] = x;
709 w[j + 1] = y;
710 w[nw - j] = y;
711 w[nw - j + 1] = x;
712 }
713 bitrv2(nw, ip + 2, w);
714 }
715 }
716
717
makect(int nc,int * ip,double * c)718 void makect(int nc, int *ip, double *c)
719 {
720 int nch, j;
721 double delta;
722
723 ip[1] = nc;
724 if (nc > 1) {
725 nch = nc >> 1;
726 delta = atan(1.0) / nch;
727 c[0] = 0.5;
728 c[nch] = 0.5 * cos(delta * nch);
729 for (j = 1; j <= nch - 1; j++) {
730 c[j] = 0.5 * cos(delta * j);
731 c[nc - j] = 0.5 * sin(delta * j);
732 }
733 }
734 }
735
736
737 /* -------- child routines -------- */
738
739
bitrv2(int n,int * ip,double * a)740 void bitrv2(int n, int *ip, double *a)
741 {
742 int j, j1, k, k1, l, m, m2;
743 double xr, xi;
744
745 ip[0] = 0;
746 l = n;
747 m = 1;
748 while ((m << 2) < l) {
749 l >>= 1;
750 for (j = 0; j <= m - 1; j++) {
751 ip[m + j] = ip[j] + l;
752 }
753 m <<= 1;
754 }
755 if ((m << 2) > l) {
756 for (k = 1; k <= m - 1; k++) {
757 for (j = 0; j <= k - 1; j++) {
758 j1 = (j << 1) + ip[k];
759 k1 = (k << 1) + ip[j];
760 xr = a[j1];
761 xi = a[j1 + 1];
762 a[j1] = a[k1];
763 a[j1 + 1] = a[k1 + 1];
764 a[k1] = xr;
765 a[k1 + 1] = xi;
766 }
767 }
768 } else {
769 m2 = m << 1;
770 for (k = 1; k <= m - 1; k++) {
771 for (j = 0; j <= k - 1; j++) {
772 j1 = (j << 1) + ip[k];
773 k1 = (k << 1) + ip[j];
774 xr = a[j1];
775 xi = a[j1 + 1];
776 a[j1] = a[k1];
777 a[j1 + 1] = a[k1 + 1];
778 a[k1] = xr;
779 a[k1 + 1] = xi;
780 j1 += m2;
781 k1 += m2;
782 xr = a[j1];
783 xi = a[j1 + 1];
784 a[j1] = a[k1];
785 a[j1 + 1] = a[k1 + 1];
786 a[k1] = xr;
787 a[k1 + 1] = xi;
788 }
789 }
790 }
791 }
792
793
bitrv2col(int n1,int n,int * ip,double ** a)794 void bitrv2col(int n1, int n, int *ip, double **a)
795 {
796 int i, j, j1, k, k1, l, m, m2;
797 double xr, xi;
798
799 ip[0] = 0;
800 l = n;
801 m = 1;
802 while ((m << 2) < l) {
803 l >>= 1;
804 for (j = 0; j <= m - 1; j++) {
805 ip[m + j] = ip[j] + l;
806 }
807 m <<= 1;
808 }
809 if ((m << 2) > l) {
810 for (i = 0; i <= n1 - 1; i++) {
811 for (k = 1; k <= m - 1; k++) {
812 for (j = 0; j <= k - 1; j++) {
813 j1 = (j << 1) + ip[k];
814 k1 = (k << 1) + ip[j];
815 xr = a[i][j1];
816 xi = a[i][j1 + 1];
817 a[i][j1] = a[i][k1];
818 a[i][j1 + 1] = a[i][k1 + 1];
819 a[i][k1] = xr;
820 a[i][k1 + 1] = xi;
821 }
822 }
823 }
824 } else {
825 m2 = m << 1;
826 for (i = 0; i <= n1 - 1; i++) {
827 for (k = 1; k <= m - 1; k++) {
828 for (j = 0; j <= k - 1; j++) {
829 j1 = (j << 1) + ip[k];
830 k1 = (k << 1) + ip[j];
831 xr = a[i][j1];
832 xi = a[i][j1 + 1];
833 a[i][j1] = a[i][k1];
834 a[i][j1 + 1] = a[i][k1 + 1];
835 a[i][k1] = xr;
836 a[i][k1 + 1] = xi;
837 j1 += m2;
838 k1 += m2;
839 xr = a[i][j1];
840 xi = a[i][j1 + 1];
841 a[i][j1] = a[i][k1];
842 a[i][j1 + 1] = a[i][k1 + 1];
843 a[i][k1] = xr;
844 a[i][k1 + 1] = xi;
845 }
846 }
847 }
848 }
849 }
850
851
bitrv2row(int n,int n2,int * ip,double ** a)852 void bitrv2row(int n, int n2, int *ip, double **a)
853 {
854 int i, j, j1, k, k1, l, m;
855 double xr, xi;
856
857 ip[0] = 0;
858 l = n;
859 m = 1;
860 while ((m << 1) < l) {
861 l >>= 1;
862 for (j = 0; j <= m - 1; j++) {
863 ip[m + j] = ip[j] + l;
864 }
865 m <<= 1;
866 }
867 if ((m << 1) > l) {
868 for (k = 1; k <= m - 1; k++) {
869 for (j = 0; j <= k - 1; j++) {
870 j1 = j + ip[k];
871 k1 = k + ip[j];
872 for (i = 0; i <= n2 - 2; i += 2) {
873 xr = a[j1][i];
874 xi = a[j1][i + 1];
875 a[j1][i] = a[k1][i];
876 a[j1][i + 1] = a[k1][i + 1];
877 a[k1][i] = xr;
878 a[k1][i + 1] = xi;
879 }
880 }
881 }
882 } else {
883 for (k = 1; k <= m - 1; k++) {
884 for (j = 0; j <= k - 1; j++) {
885 j1 = j + ip[k];
886 k1 = k + ip[j];
887 for (i = 0; i <= n2 - 2; i += 2) {
888 xr = a[j1][i];
889 xi = a[j1][i + 1];
890 a[j1][i] = a[k1][i];
891 a[j1][i + 1] = a[k1][i + 1];
892 a[k1][i] = xr;
893 a[k1][i + 1] = xi;
894 }
895 j1 += m;
896 k1 += m;
897 for (i = 0; i <= n2 - 2; i += 2) {
898 xr = a[j1][i];
899 xi = a[j1][i + 1];
900 a[j1][i] = a[k1][i];
901 a[j1][i + 1] = a[k1][i + 1];
902 a[k1][i] = xr;
903 a[k1][i + 1] = xi;
904 }
905 }
906 }
907 }
908 }
909
910
cftbcol(int n1,int n,double ** a,double * w)911 void cftbcol(int n1, int n, double **a, double *w)
912 {
913 int i, j, j1, j2, j3, k, k1, ks, l, m;
914 double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
915 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
916
917 for (i = 0; i <= n1 - 1; i++) {
918 l = 2;
919 while ((l << 1) < n) {
920 m = l << 2;
921 for (j = 0; j <= l - 2; j += 2) {
922 j1 = j + l;
923 j2 = j1 + l;
924 j3 = j2 + l;
925 x0r = a[i][j] + a[i][j1];
926 x0i = a[i][j + 1] + a[i][j1 + 1];
927 x1r = a[i][j] - a[i][j1];
928 x1i = a[i][j + 1] - a[i][j1 + 1];
929 x2r = a[i][j2] + a[i][j3];
930 x2i = a[i][j2 + 1] + a[i][j3 + 1];
931 x3r = a[i][j2] - a[i][j3];
932 x3i = a[i][j2 + 1] - a[i][j3 + 1];
933 a[i][j] = x0r + x2r;
934 a[i][j + 1] = x0i + x2i;
935 a[i][j2] = x0r - x2r;
936 a[i][j2 + 1] = x0i - x2i;
937 a[i][j1] = x1r - x3i;
938 a[i][j1 + 1] = x1i + x3r;
939 a[i][j3] = x1r + x3i;
940 a[i][j3 + 1] = x1i - x3r;
941 }
942 if (m < n) {
943 wk1r = w[2];
944 for (j = m; j <= l + m - 2; j += 2) {
945 j1 = j + l;
946 j2 = j1 + l;
947 j3 = j2 + l;
948 x0r = a[i][j] + a[i][j1];
949 x0i = a[i][j + 1] + a[i][j1 + 1];
950 x1r = a[i][j] - a[i][j1];
951 x1i = a[i][j + 1] - a[i][j1 + 1];
952 x2r = a[i][j2] + a[i][j3];
953 x2i = a[i][j2 + 1] + a[i][j3 + 1];
954 x3r = a[i][j2] - a[i][j3];
955 x3i = a[i][j2 + 1] - a[i][j3 + 1];
956 a[i][j] = x0r + x2r;
957 a[i][j + 1] = x0i + x2i;
958 a[i][j2] = x2i - x0i;
959 a[i][j2 + 1] = x0r - x2r;
960 x0r = x1r - x3i;
961 x0i = x1i + x3r;
962 a[i][j1] = wk1r * (x0r - x0i);
963 a[i][j1 + 1] = wk1r * (x0r + x0i);
964 x0r = x3i + x1r;
965 x0i = x3r - x1i;
966 a[i][j3] = wk1r * (x0i - x0r);
967 a[i][j3 + 1] = wk1r * (x0i + x0r);
968 }
969 k1 = 1;
970 ks = -1;
971 for (k = (m << 1); k <= n - m; k += m) {
972 k1++;
973 ks = -ks;
974 wk1r = w[k1 << 1];
975 wk1i = w[(k1 << 1) + 1];
976 wk2r = ks * w[k1];
977 wk2i = w[k1 + ks];
978 wk3r = wk1r - 2 * wk2i * wk1i;
979 wk3i = 2 * wk2i * wk1r - wk1i;
980 for (j = k; j <= l + k - 2; j += 2) {
981 j1 = j + l;
982 j2 = j1 + l;
983 j3 = j2 + l;
984 x0r = a[i][j] + a[i][j1];
985 x0i = a[i][j + 1] + a[i][j1 + 1];
986 x1r = a[i][j] - a[i][j1];
987 x1i = a[i][j + 1] - a[i][j1 + 1];
988 x2r = a[i][j2] + a[i][j3];
989 x2i = a[i][j2 + 1] + a[i][j3 + 1];
990 x3r = a[i][j2] - a[i][j3];
991 x3i = a[i][j2 + 1] - a[i][j3 + 1];
992 a[i][j] = x0r + x2r;
993 a[i][j + 1] = x0i + x2i;
994 x0r -= x2r;
995 x0i -= x2i;
996 a[i][j2] = wk2r * x0r - wk2i * x0i;
997 a[i][j2 + 1] = wk2r * x0i + wk2i * x0r;
998 x0r = x1r - x3i;
999 x0i = x1i + x3r;
1000 a[i][j1] = wk1r * x0r - wk1i * x0i;
1001 a[i][j1 + 1] = wk1r * x0i + wk1i * x0r;
1002 x0r = x1r + x3i;
1003 x0i = x1i - x3r;
1004 a[i][j3] = wk3r * x0r - wk3i * x0i;
1005 a[i][j3 + 1] = wk3r * x0i + wk3i * x0r;
1006 }
1007 }
1008 }
1009 l = m;
1010 }
1011 if (l < n) {
1012 for (j = 0; j <= l - 2; j += 2) {
1013 j1 = j + l;
1014 x0r = a[i][j] - a[i][j1];
1015 x0i = a[i][j + 1] - a[i][j1 + 1];
1016 a[i][j] += a[i][j1];
1017 a[i][j + 1] += a[i][j1 + 1];
1018 a[i][j1] = x0r;
1019 a[i][j1 + 1] = x0i;
1020 }
1021 }
1022 }
1023 }
1024
1025
cftbrow(int n,int n2,double ** a,double * w)1026 void cftbrow(int n, int n2, double **a, double *w)
1027 {
1028 int i, j, j1, j2, j3, k, k1, ks, l, m;
1029 double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
1030 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
1031
1032 l = 1;
1033 while ((l << 1) < n) {
1034 m = l << 2;
1035 for (j = 0; j <= l - 1; j++) {
1036 j1 = j + l;
1037 j2 = j1 + l;
1038 j3 = j2 + l;
1039 for (i = 0; i <= n2 - 2; i += 2) {
1040 x0r = a[j][i] + a[j1][i];
1041 x0i = a[j][i + 1] + a[j1][i + 1];
1042 x1r = a[j][i] - a[j1][i];
1043 x1i = a[j][i + 1] - a[j1][i + 1];
1044 x2r = a[j2][i] + a[j3][i];
1045 x2i = a[j2][i + 1] + a[j3][i + 1];
1046 x3r = a[j2][i] - a[j3][i];
1047 x3i = a[j2][i + 1] - a[j3][i + 1];
1048 a[j][i] = x0r + x2r;
1049 a[j][i + 1] = x0i + x2i;
1050 a[j2][i] = x0r - x2r;
1051 a[j2][i + 1] = x0i - x2i;
1052 a[j1][i] = x1r - x3i;
1053 a[j1][i + 1] = x1i + x3r;
1054 a[j3][i] = x1r + x3i;
1055 a[j3][i + 1] = x1i - x3r;
1056 }
1057 }
1058 if (m < n) {
1059 wk1r = w[2];
1060 for (j = m; j <= l + m - 1; j++) {
1061 j1 = j + l;
1062 j2 = j1 + l;
1063 j3 = j2 + l;
1064 for (i = 0; i <= n2 - 2; i += 2) {
1065 x0r = a[j][i] + a[j1][i];
1066 x0i = a[j][i + 1] + a[j1][i + 1];
1067 x1r = a[j][i] - a[j1][i];
1068 x1i = a[j][i + 1] - a[j1][i + 1];
1069 x2r = a[j2][i] + a[j3][i];
1070 x2i = a[j2][i + 1] + a[j3][i + 1];
1071 x3r = a[j2][i] - a[j3][i];
1072 x3i = a[j2][i + 1] - a[j3][i + 1];
1073 a[j][i] = x0r + x2r;
1074 a[j][i + 1] = x0i + x2i;
1075 a[j2][i] = x2i - x0i;
1076 a[j2][i + 1] = x0r - x2r;
1077 x0r = x1r - x3i;
1078 x0i = x1i + x3r;
1079 a[j1][i] = wk1r * (x0r - x0i);
1080 a[j1][i + 1] = wk1r * (x0r + x0i);
1081 x0r = x3i + x1r;
1082 x0i = x3r - x1i;
1083 a[j3][i] = wk1r * (x0i - x0r);
1084 a[j3][i + 1] = wk1r * (x0i + x0r);
1085 }
1086 }
1087 k1 = 1;
1088 ks = -1;
1089 for (k = (m << 1); k <= n - m; k += m) {
1090 k1++;
1091 ks = -ks;
1092 wk1r = w[k1 << 1];
1093 wk1i = w[(k1 << 1) + 1];
1094 wk2r = ks * w[k1];
1095 wk2i = w[k1 + ks];
1096 wk3r = wk1r - 2 * wk2i * wk1i;
1097 wk3i = 2 * wk2i * wk1r - wk1i;
1098 for (j = k; j <= l + k - 1; j++) {
1099 j1 = j + l;
1100 j2 = j1 + l;
1101 j3 = j2 + l;
1102 for (i = 0; i <= n2 - 2; i += 2) {
1103 x0r = a[j][i] + a[j1][i];
1104 x0i = a[j][i + 1] + a[j1][i + 1];
1105 x1r = a[j][i] - a[j1][i];
1106 x1i = a[j][i + 1] - a[j1][i + 1];
1107 x2r = a[j2][i] + a[j3][i];
1108 x2i = a[j2][i + 1] + a[j3][i + 1];
1109 x3r = a[j2][i] - a[j3][i];
1110 x3i = a[j2][i + 1] - a[j3][i + 1];
1111 a[j][i] = x0r + x2r;
1112 a[j][i + 1] = x0i + x2i;
1113 x0r -= x2r;
1114 x0i -= x2i;
1115 a[j2][i] = wk2r * x0r - wk2i * x0i;
1116 a[j2][i + 1] = wk2r * x0i + wk2i * x0r;
1117 x0r = x1r - x3i;
1118 x0i = x1i + x3r;
1119 a[j1][i] = wk1r * x0r - wk1i * x0i;
1120 a[j1][i + 1] = wk1r * x0i + wk1i * x0r;
1121 x0r = x1r + x3i;
1122 x0i = x1i - x3r;
1123 a[j3][i] = wk3r * x0r - wk3i * x0i;
1124 a[j3][i + 1] = wk3r * x0i + wk3i * x0r;
1125 }
1126 }
1127 }
1128 }
1129 l = m;
1130 }
1131 if (l < n) {
1132 for (j = 0; j <= l - 1; j++) {
1133 j1 = j + l;
1134 for (i = 0; i <= n2 - 2; i += 2) {
1135 x0r = a[j][i] - a[j1][i];
1136 x0i = a[j][i + 1] - a[j1][i + 1];
1137 a[j][i] += a[j1][i];
1138 a[j][i + 1] += a[j1][i + 1];
1139 a[j1][i] = x0r;
1140 a[j1][i + 1] = x0i;
1141 }
1142 }
1143 }
1144 }
1145
1146
cftfcol(int n1,int n,double ** a,double * w)1147 void cftfcol(int n1, int n, double **a, double *w)
1148 {
1149 int i, j, j1, j2, j3, k, k1, ks, l, m;
1150 double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
1151 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
1152
1153 for (i = 0; i <= n1 - 1; i++) {
1154 l = 2;
1155 while ((l << 1) < n) {
1156 m = l << 2;
1157 for (j = 0; j <= l - 2; j += 2) {
1158 j1 = j + l;
1159 j2 = j1 + l;
1160 j3 = j2 + l;
1161 x0r = a[i][j] + a[i][j1];
1162 x0i = a[i][j + 1] + a[i][j1 + 1];
1163 x1r = a[i][j] - a[i][j1];
1164 x1i = a[i][j + 1] - a[i][j1 + 1];
1165 x2r = a[i][j2] + a[i][j3];
1166 x2i = a[i][j2 + 1] + a[i][j3 + 1];
1167 x3r = a[i][j2] - a[i][j3];
1168 x3i = a[i][j2 + 1] - a[i][j3 + 1];
1169 a[i][j] = x0r + x2r;
1170 a[i][j + 1] = x0i + x2i;
1171 a[i][j2] = x0r - x2r;
1172 a[i][j2 + 1] = x0i - x2i;
1173 a[i][j1] = x1r + x3i;
1174 a[i][j1 + 1] = x1i - x3r;
1175 a[i][j3] = x1r - x3i;
1176 a[i][j3 + 1] = x1i + x3r;
1177 }
1178 if (m < n) {
1179 wk1r = w[2];
1180 for (j = m; j <= l + m - 2; j += 2) {
1181 j1 = j + l;
1182 j2 = j1 + l;
1183 j3 = j2 + l;
1184 x0r = a[i][j] + a[i][j1];
1185 x0i = a[i][j + 1] + a[i][j1 + 1];
1186 x1r = a[i][j] - a[i][j1];
1187 x1i = a[i][j + 1] - a[i][j1 + 1];
1188 x2r = a[i][j2] + a[i][j3];
1189 x2i = a[i][j2 + 1] + a[i][j3 + 1];
1190 x3r = a[i][j2] - a[i][j3];
1191 x3i = a[i][j2 + 1] - a[i][j3 + 1];
1192 a[i][j] = x0r + x2r;
1193 a[i][j + 1] = x0i + x2i;
1194 a[i][j2] = x0i - x2i;
1195 a[i][j2 + 1] = x2r - x0r;
1196 x0r = x1r + x3i;
1197 x0i = x1i - x3r;
1198 a[i][j1] = wk1r * (x0i + x0r);
1199 a[i][j1 + 1] = wk1r * (x0i - x0r);
1200 x0r = x3i - x1r;
1201 x0i = x3r + x1i;
1202 a[i][j3] = wk1r * (x0r + x0i);
1203 a[i][j3 + 1] = wk1r * (x0r - x0i);
1204 }
1205 k1 = 1;
1206 ks = -1;
1207 for (k = (m << 1); k <= n - m; k += m) {
1208 k1++;
1209 ks = -ks;
1210 wk1r = w[k1 << 1];
1211 wk1i = w[(k1 << 1) + 1];
1212 wk2r = ks * w[k1];
1213 wk2i = w[k1 + ks];
1214 wk3r = wk1r - 2 * wk2i * wk1i;
1215 wk3i = 2 * wk2i * wk1r - wk1i;
1216 for (j = k; j <= l + k - 2; j += 2) {
1217 j1 = j + l;
1218 j2 = j1 + l;
1219 j3 = j2 + l;
1220 x0r = a[i][j] + a[i][j1];
1221 x0i = a[i][j + 1] + a[i][j1 + 1];
1222 x1r = a[i][j] - a[i][j1];
1223 x1i = a[i][j + 1] - a[i][j1 + 1];
1224 x2r = a[i][j2] + a[i][j3];
1225 x2i = a[i][j2 + 1] + a[i][j3 + 1];
1226 x3r = a[i][j2] - a[i][j3];
1227 x3i = a[i][j2 + 1] - a[i][j3 + 1];
1228 a[i][j] = x0r + x2r;
1229 a[i][j + 1] = x0i + x2i;
1230 x0r -= x2r;
1231 x0i -= x2i;
1232 a[i][j2] = wk2r * x0r + wk2i * x0i;
1233 a[i][j2 + 1] = wk2r * x0i - wk2i * x0r;
1234 x0r = x1r + x3i;
1235 x0i = x1i - x3r;
1236 a[i][j1] = wk1r * x0r + wk1i * x0i;
1237 a[i][j1 + 1] = wk1r * x0i - wk1i * x0r;
1238 x0r = x1r - x3i;
1239 x0i = x1i + x3r;
1240 a[i][j3] = wk3r * x0r + wk3i * x0i;
1241 a[i][j3 + 1] = wk3r * x0i - wk3i * x0r;
1242 }
1243 }
1244 }
1245 l = m;
1246 }
1247 if (l < n) {
1248 for (j = 0; j <= l - 2; j += 2) {
1249 j1 = j + l;
1250 x0r = a[i][j] - a[i][j1];
1251 x0i = a[i][j + 1] - a[i][j1 + 1];
1252 a[i][j] += a[i][j1];
1253 a[i][j + 1] += a[i][j1 + 1];
1254 a[i][j1] = x0r;
1255 a[i][j1 + 1] = x0i;
1256 }
1257 }
1258 }
1259 }
1260
1261
cftfrow(int n,int n2,double ** a,double * w)1262 void cftfrow(int n, int n2, double **a, double *w)
1263 {
1264 int i, j, j1, j2, j3, k, k1, ks, l, m;
1265 double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
1266 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
1267
1268 l = 1;
1269 while ((l << 1) < n) {
1270 m = l << 2;
1271 for (j = 0; j <= l - 1; j++) {
1272 j1 = j + l;
1273 j2 = j1 + l;
1274 j3 = j2 + l;
1275 for (i = 0; i <= n2 - 2; i += 2) {
1276 x0r = a[j][i] + a[j1][i];
1277 x0i = a[j][i + 1] + a[j1][i + 1];
1278 x1r = a[j][i] - a[j1][i];
1279 x1i = a[j][i + 1] - a[j1][i + 1];
1280 x2r = a[j2][i] + a[j3][i];
1281 x2i = a[j2][i + 1] + a[j3][i + 1];
1282 x3r = a[j2][i] - a[j3][i];
1283 x3i = a[j2][i + 1] - a[j3][i + 1];
1284 a[j][i] = x0r + x2r;
1285 a[j][i + 1] = x0i + x2i;
1286 a[j2][i] = x0r - x2r;
1287 a[j2][i + 1] = x0i - x2i;
1288 a[j1][i] = x1r + x3i;
1289 a[j1][i + 1] = x1i - x3r;
1290 a[j3][i] = x1r - x3i;
1291 a[j3][i + 1] = x1i + x3r;
1292 }
1293 }
1294 if (m < n) {
1295 wk1r = w[2];
1296 for (j = m; j <= l + m - 1; j++) {
1297 j1 = j + l;
1298 j2 = j1 + l;
1299 j3 = j2 + l;
1300 for (i = 0; i <= n2 - 2; i += 2) {
1301 x0r = a[j][i] + a[j1][i];
1302 x0i = a[j][i + 1] + a[j1][i + 1];
1303 x1r = a[j][i] - a[j1][i];
1304 x1i = a[j][i + 1] - a[j1][i + 1];
1305 x2r = a[j2][i] + a[j3][i];
1306 x2i = a[j2][i + 1] + a[j3][i + 1];
1307 x3r = a[j2][i] - a[j3][i];
1308 x3i = a[j2][i + 1] - a[j3][i + 1];
1309 a[j][i] = x0r + x2r;
1310 a[j][i + 1] = x0i + x2i;
1311 a[j2][i] = x0i - x2i;
1312 a[j2][i + 1] = x2r - x0r;
1313 x0r = x1r + x3i;
1314 x0i = x1i - x3r;
1315 a[j1][i] = wk1r * (x0i + x0r);
1316 a[j1][i + 1] = wk1r * (x0i - x0r);
1317 x0r = x3i - x1r;
1318 x0i = x3r + x1i;
1319 a[j3][i] = wk1r * (x0r + x0i);
1320 a[j3][i + 1] = wk1r * (x0r - x0i);
1321 }
1322 }
1323 k1 = 1;
1324 ks = -1;
1325 for (k = (m << 1); k <= n - m; k += m) {
1326 k1++;
1327 ks = -ks;
1328 wk1r = w[k1 << 1];
1329 wk1i = w[(k1 << 1) + 1];
1330 wk2r = ks * w[k1];
1331 wk2i = w[k1 + ks];
1332 wk3r = wk1r - 2 * wk2i * wk1i;
1333 wk3i = 2 * wk2i * wk1r - wk1i;
1334 for (j = k; j <= l + k - 1; j++) {
1335 j1 = j + l;
1336 j2 = j1 + l;
1337 j3 = j2 + l;
1338 for (i = 0; i <= n2 - 2; i += 2) {
1339 x0r = a[j][i] + a[j1][i];
1340 x0i = a[j][i + 1] + a[j1][i + 1];
1341 x1r = a[j][i] - a[j1][i];
1342 x1i = a[j][i + 1] - a[j1][i + 1];
1343 x2r = a[j2][i] + a[j3][i];
1344 x2i = a[j2][i + 1] + a[j3][i + 1];
1345 x3r = a[j2][i] - a[j3][i];
1346 x3i = a[j2][i + 1] - a[j3][i + 1];
1347 a[j][i] = x0r + x2r;
1348 a[j][i + 1] = x0i + x2i;
1349 x0r -= x2r;
1350 x0i -= x2i;
1351 a[j2][i] = wk2r * x0r + wk2i * x0i;
1352 a[j2][i + 1] = wk2r * x0i - wk2i * x0r;
1353 x0r = x1r + x3i;
1354 x0i = x1i - x3r;
1355 a[j1][i] = wk1r * x0r + wk1i * x0i;
1356 a[j1][i + 1] = wk1r * x0i - wk1i * x0r;
1357 x0r = x1r - x3i;
1358 x0i = x1i + x3r;
1359 a[j3][i] = wk3r * x0r + wk3i * x0i;
1360 a[j3][i + 1] = wk3r * x0i - wk3i * x0r;
1361 }
1362 }
1363 }
1364 }
1365 l = m;
1366 }
1367 if (l < n) {
1368 for (j = 0; j <= l - 1; j++) {
1369 j1 = j + l;
1370 for (i = 0; i <= n2 - 2; i += 2) {
1371 x0r = a[j][i] - a[j1][i];
1372 x0i = a[j][i + 1] - a[j1][i + 1];
1373 a[j][i] += a[j1][i];
1374 a[j][i + 1] += a[j1][i + 1];
1375 a[j1][i] = x0r;
1376 a[j1][i + 1] = x0i;
1377 }
1378 }
1379 }
1380 }
1381
1382
rftbcol(int n1,int n,double ** a,int nc,double * c)1383 void rftbcol(int n1, int n, double **a, int nc, double *c)
1384 {
1385 int i, j, k, kk, ks;
1386 double wkr, wki, xr, xi, yr, yi;
1387
1388 ks = (nc << 2) / n;
1389 for (i = 0; i <= n1 - 1; i++) {
1390 kk = 0;
1391 for (k = (n >> 1) - 2; k >= 2; k -= 2) {
1392 j = n - k;
1393 kk += ks;
1394 wkr = 0.5 - c[kk];
1395 wki = c[nc - kk];
1396 xr = a[i][k] - a[i][j];
1397 xi = a[i][k + 1] + a[i][j + 1];
1398 yr = wkr * xr - wki * xi;
1399 yi = wkr * xi + wki * xr;
1400 a[i][k] -= yr;
1401 a[i][k + 1] -= yi;
1402 a[i][j] += yr;
1403 a[i][j + 1] -= yi;
1404 }
1405 }
1406 }
1407
1408
rftfcol(int n1,int n,double ** a,int nc,double * c)1409 void rftfcol(int n1, int n, double **a, int nc, double *c)
1410 {
1411 int i, j, k, kk, ks;
1412 double wkr, wki, xr, xi, yr, yi;
1413
1414 ks = (nc << 2) / n;
1415 for (i = 0; i <= n1 - 1; i++) {
1416 kk = 0;
1417 for (k = (n >> 1) - 2; k >= 2; k -= 2) {
1418 j = n - k;
1419 kk += ks;
1420 wkr = 0.5 - c[kk];
1421 wki = c[nc - kk];
1422 xr = a[i][k] - a[i][j];
1423 xi = a[i][k + 1] + a[i][j + 1];
1424 yr = wkr * xr + wki * xi;
1425 yi = wkr * xi - wki * xr;
1426 a[i][k] -= yr;
1427 a[i][k + 1] -= yi;
1428 a[i][j] += yr;
1429 a[i][j + 1] -= yi;
1430 }
1431 }
1432 }
1433
1434
dctbsub(int n1,int n2,double ** a,int nc,double * c)1435 void dctbsub(int n1, int n2, double **a, int nc, double *c)
1436 {
1437 int kk1, kk2, ks1, ks2, n1h, j1, k1, k2;
1438 double w1r, w1i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i;
1439
1440 ks1 = nc / n1;
1441 ks2 = nc / n2;
1442 n1h = n1 >> 1;
1443 kk1 = ks1;
1444 for (k1 = 1; k1 <= n1h - 1; k1++) {
1445 j1 = n1 - k1;
1446 w1r = 2 * c[kk1];
1447 w1i = 2 * c[nc - kk1];
1448 kk1 += ks1;
1449 kk2 = ks2;
1450 for (k2 = 2; k2 <= n2 - 2; k2 += 2) {
1451 x0r = w1r * c[kk2];
1452 x0i = w1i * c[kk2];
1453 x1r = w1r * c[nc - kk2];
1454 x1i = w1i * c[nc - kk2];
1455 wkr = x0r - x1i;
1456 wki = x0i + x1r;
1457 wji = x0r + x1i;
1458 wjr = x0i - x1r;
1459 kk2 += ks2;
1460 x0r = wkr * a[k1][k2] - wki * a[k1][k2 + 1];
1461 x0i = wkr * a[k1][k2 + 1] + wki * a[k1][k2];
1462 x1r = wjr * a[j1][k2] - wji * a[j1][k2 + 1];
1463 x1i = wjr * a[j1][k2 + 1] + wji * a[j1][k2];
1464 a[k1][k2] = x0r + x1i;
1465 a[k1][k2 + 1] = x0i - x1r;
1466 a[j1][k2] = x1r + x0i;
1467 a[j1][k2 + 1] = x1i - x0r;
1468 }
1469 wkr = w1r * 0.5;
1470 wki = w1i * 0.5;
1471 wjr = w1r * c[kk2];
1472 wji = w1i * c[kk2];
1473 x0r = a[k1][0] + a[j1][0];
1474 x0i = a[k1][1] - a[j1][1];
1475 x1r = a[k1][0] - a[j1][0];
1476 x1i = a[k1][1] + a[j1][1];
1477 a[k1][0] = wkr * x0r - wki * x0i;
1478 a[k1][1] = wkr * x0i + wki * x0r;
1479 a[j1][0] = -wjr * x1r + wji * x1i;
1480 a[j1][1] = wjr * x1i + wji * x1r;
1481 }
1482 w1r = 2 * c[kk1];
1483 kk2 = ks2;
1484 for (k2 = 2; k2 <= n2 - 2; k2 += 2) {
1485 wkr = 2 * c[kk2];
1486 wki = 2 * c[nc - kk2];
1487 wjr = w1r * wkr;
1488 wji = w1r * wki;
1489 kk2 += ks2;
1490 x0i = wkr * a[0][k2 + 1] + wki * a[0][k2];
1491 a[0][k2] = wkr * a[0][k2] - wki * a[0][k2 + 1];
1492 a[0][k2 + 1] = x0i;
1493 x0i = wjr * a[n1h][k2 + 1] + wji * a[n1h][k2];
1494 a[n1h][k2] = wjr * a[n1h][k2] - wji * a[n1h][k2 + 1];
1495 a[n1h][k2 + 1] = x0i;
1496 }
1497 a[0][1] *= w1r;
1498 a[n1h][0] *= w1r;
1499 a[n1h][1] *= 0.5;
1500 }
1501
1502
dctfsub(int n1,int n2,double ** a,int nc,double * c)1503 void dctfsub(int n1, int n2, double **a, int nc, double *c)
1504 {
1505 int kk1, kk2, ks1, ks2, n1h, j1, k1, k2;
1506 double w1r, w1i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i;
1507
1508 ks1 = nc / n1;
1509 ks2 = nc / n2;
1510 n1h = n1 >> 1;
1511 kk1 = ks1;
1512 for (k1 = 1; k1 <= n1h - 1; k1++) {
1513 j1 = n1 - k1;
1514 w1r = 2 * c[kk1];
1515 w1i = 2 * c[nc - kk1];
1516 kk1 += ks1;
1517 kk2 = ks2;
1518 for (k2 = 2; k2 <= n2 - 2; k2 += 2) {
1519 x0r = w1r * c[kk2];
1520 x0i = w1i * c[kk2];
1521 x1r = w1r * c[nc - kk2];
1522 x1i = w1i * c[nc - kk2];
1523 wkr = x0r - x1i;
1524 wki = x0i + x1r;
1525 wji = x0r + x1i;
1526 wjr = x0i - x1r;
1527 kk2 += ks2;
1528 x0r = a[k1][k2] - a[j1][k2 + 1];
1529 x0i = a[j1][k2] + a[k1][k2 + 1];
1530 x1r = a[j1][k2] - a[k1][k2 + 1];
1531 x1i = a[k1][k2] + a[j1][k2 + 1];
1532 a[k1][k2] = wkr * x0r + wki * x0i;
1533 a[k1][k2 + 1] = wkr * x0i - wki * x0r;
1534 a[j1][k2] = wjr * x1r + wji * x1i;
1535 a[j1][k2 + 1] = wjr * x1i - wji * x1r;
1536 }
1537 x0r = 2 * c[kk2];
1538 wjr = x0r * w1r;
1539 wji = x0r * w1i;
1540 x0r = w1r * a[k1][0] + w1i * a[k1][1];
1541 x0i = w1r * a[k1][1] - w1i * a[k1][0];
1542 x1r = -wjr * a[j1][0] + wji * a[j1][1];
1543 x1i = wjr * a[j1][1] + wji * a[j1][0];
1544 a[k1][0] = x0r + x1r;
1545 a[k1][1] = x1i + x0i;
1546 a[j1][0] = x0r - x1r;
1547 a[j1][1] = x1i - x0i;
1548 }
1549 w1r = 2 * c[kk1];
1550 kk2 = ks2;
1551 for (k2 = 2; k2 <= n2 - 2; k2 += 2) {
1552 wkr = 2 * c[kk2];
1553 wki = 2 * c[nc - kk2];
1554 wjr = w1r * wkr;
1555 wji = w1r * wki;
1556 kk2 += ks2;
1557 x0i = wkr * a[0][k2 + 1] - wki * a[0][k2];
1558 a[0][k2] = wkr * a[0][k2] + wki * a[0][k2 + 1];
1559 a[0][k2 + 1] = x0i;
1560 x0i = wjr * a[n1h][k2 + 1] - wji * a[n1h][k2];
1561 a[n1h][k2] = wjr * a[n1h][k2] + wji * a[n1h][k2 + 1];
1562 a[n1h][k2 + 1] = x0i;
1563 }
1564 w1r *= 2;
1565 a[0][0] *= 2;
1566 a[0][1] *= w1r;
1567 a[n1h][0] *= w1r;
1568 }
1569
1570
dstbsub(int n1,int n2,double ** a,int nc,double * c)1571 void dstbsub(int n1, int n2, double **a, int nc, double *c)
1572 {
1573 int kk1, kk2, ks1, ks2, n1h, j1, k1, k2;
1574 double w1r, w1i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i;
1575
1576 ks1 = nc / n1;
1577 ks2 = nc / n2;
1578 n1h = n1 >> 1;
1579 kk1 = ks1;
1580 for (k1 = 1; k1 <= n1h - 1; k1++) {
1581 j1 = n1 - k1;
1582 w1r = 2 * c[kk1];
1583 w1i = 2 * c[nc - kk1];
1584 kk1 += ks1;
1585 kk2 = ks2;
1586 for (k2 = 2; k2 <= n2 - 2; k2 += 2) {
1587 x0r = w1r * c[kk2];
1588 x0i = w1i * c[kk2];
1589 x1r = w1r * c[nc - kk2];
1590 x1i = w1i * c[nc - kk2];
1591 wkr = x0r - x1i;
1592 wki = x0i + x1r;
1593 wji = x0r + x1i;
1594 wjr = x0i - x1r;
1595 kk2 += ks2;
1596 x0r = wkr * a[k1][k2] - wki * a[k1][k2 + 1];
1597 x0i = wkr * a[k1][k2 + 1] + wki * a[k1][k2];
1598 x1r = wjr * a[j1][k2] - wji * a[j1][k2 + 1];
1599 x1i = wjr * a[j1][k2 + 1] + wji * a[j1][k2];
1600 a[k1][k2] = x1i - x0r;
1601 a[k1][k2 + 1] = x1r + x0i;
1602 a[j1][k2] = x0i - x1r;
1603 a[j1][k2 + 1] = x0r + x1i;
1604 }
1605 wkr = w1r * 0.5;
1606 wki = w1i * 0.5;
1607 wjr = w1r * c[kk2];
1608 wji = w1i * c[kk2];
1609 x0r = a[k1][0] + a[j1][0];
1610 x0i = a[k1][1] - a[j1][1];
1611 x1r = a[k1][0] - a[j1][0];
1612 x1i = a[k1][1] + a[j1][1];
1613 a[k1][1] = wkr * x0r - wki * x0i;
1614 a[k1][0] = wkr * x0i + wki * x0r;
1615 a[j1][1] = -wjr * x1r + wji * x1i;
1616 a[j1][0] = wjr * x1i + wji * x1r;
1617 }
1618 w1r = 2 * c[kk1];
1619 kk2 = ks2;
1620 for (k2 = 2; k2 <= n2 - 2; k2 += 2) {
1621 wkr = 2 * c[kk2];
1622 wki = 2 * c[nc - kk2];
1623 wjr = w1r * wkr;
1624 wji = w1r * wki;
1625 kk2 += ks2;
1626 x0i = wkr * a[0][k2 + 1] + wki * a[0][k2];
1627 a[0][k2 + 1] = wkr * a[0][k2] - wki * a[0][k2 + 1];
1628 a[0][k2] = x0i;
1629 x0i = wjr * a[n1h][k2 + 1] + wji * a[n1h][k2];
1630 a[n1h][k2 + 1] = wjr * a[n1h][k2] - wji * a[n1h][k2 + 1];
1631 a[n1h][k2] = x0i;
1632 }
1633 a[0][1] *= w1r;
1634 a[n1h][0] *= w1r;
1635 a[n1h][1] *= 0.5;
1636 }
1637
1638
dstfsub(int n1,int n2,double ** a,int nc,double * c)1639 void dstfsub(int n1, int n2, double **a, int nc, double *c)
1640 {
1641 int kk1, kk2, ks1, ks2, n1h, j1, k1, k2;
1642 double w1r, w1i, wkr, wki, wjr, wji, x0r, x0i, x1r, x1i;
1643
1644 ks1 = nc / n1;
1645 ks2 = nc / n2;
1646 n1h = n1 >> 1;
1647 kk1 = ks1;
1648 for (k1 = 1; k1 <= n1h - 1; k1++) {
1649 j1 = n1 - k1;
1650 w1r = 2 * c[kk1];
1651 w1i = 2 * c[nc - kk1];
1652 kk1 += ks1;
1653 kk2 = ks2;
1654 for (k2 = 2; k2 <= n2 - 2; k2 += 2) {
1655 x0r = w1r * c[kk2];
1656 x0i = w1i * c[kk2];
1657 x1r = w1r * c[nc - kk2];
1658 x1i = w1i * c[nc - kk2];
1659 wkr = x0r - x1i;
1660 wki = x0i + x1r;
1661 wji = x0r + x1i;
1662 wjr = x0i - x1r;
1663 kk2 += ks2;
1664 x0r = a[j1][k2 + 1] - a[k1][k2];
1665 x0i = a[k1][k2 + 1] + a[j1][k2];
1666 x1r = a[k1][k2 + 1] - a[j1][k2];
1667 x1i = a[j1][k2 + 1] + a[k1][k2];
1668 a[k1][k2] = wkr * x0r + wki * x0i;
1669 a[k1][k2 + 1] = wkr * x0i - wki * x0r;
1670 a[j1][k2] = wjr * x1r + wji * x1i;
1671 a[j1][k2 + 1] = wjr * x1i - wji * x1r;
1672 }
1673 x0r = 2 * c[kk2];
1674 wjr = x0r * w1r;
1675 wji = x0r * w1i;
1676 x0r = w1r * a[k1][1] + w1i * a[k1][0];
1677 x0i = w1r * a[k1][0] - w1i * a[k1][1];
1678 x1r = -wjr * a[j1][1] + wji * a[j1][0];
1679 x1i = wjr * a[j1][0] + wji * a[j1][1];
1680 a[k1][0] = x0r + x1r;
1681 a[k1][1] = x1i + x0i;
1682 a[j1][0] = x0r - x1r;
1683 a[j1][1] = x1i - x0i;
1684 }
1685 w1r = 2 * c[kk1];
1686 kk2 = ks2;
1687 for (k2 = 2; k2 <= n2 - 2; k2 += 2) {
1688 wkr = 2 * c[kk2];
1689 wki = 2 * c[nc - kk2];
1690 wjr = w1r * wkr;
1691 wji = w1r * wki;
1692 kk2 += ks2;
1693 x0i = wkr * a[0][k2] - wki * a[0][k2 + 1];
1694 a[0][k2] = wkr * a[0][k2 + 1] + wki * a[0][k2];
1695 a[0][k2 + 1] = x0i;
1696 x0i = wjr * a[n1h][k2] - wji * a[n1h][k2 + 1];
1697 a[n1h][k2] = wjr * a[n1h][k2 + 1] + wji * a[n1h][k2];
1698 a[n1h][k2 + 1] = x0i;
1699 }
1700 w1r *= 2;
1701 a[0][0] *= 2;
1702 a[0][1] *= w1r;
1703 a[n1h][0] *= w1r;
1704 }
1705
1706