Jamoma API  0.6.0.a19
fft4f2d.cpp
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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