← Back to branch summary

~ubuntu-branches/debian/sid/pgadmin3/sid

« back to all changes in this revision

Viewing changes to pgadmin/pgscript/utilities/m_apm/mapm_fft.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Gerfried Fuchs
  • Date: 2009-07-30 12:27:16 UTC
  • mfrom: (1.1.6 upstream)
  • Revision ID: james.westby@ubuntu.com-20090730122716-fddbh42on721bbs2
Tags: 1.10.0-1
* New upstream release.
* Adjusted watch file to match release candidates.
* Updated to Standards-Version 3.8.2:
  - Moved to Section: database.
  - Add DEB_BUILD_OPTIONS support for parallel building.
  - Move from findstring to filter suggestion for DEB_BUILD_OPTIONS parsing.
* pgagent got split into its own separate source package by upstream.
* Exclude Docs.vcproj from installation.
* Move doc-base.enus from pgadmin3 to pgadmin3-data package, the files are
  in there too.

Show diffs side-by-side

added added

removed removed

Lines of Context:
pgadmin/pgscript/utilities/m_apm/mapm_fft.cpp
 
1
 
 
2
/* 
 
3
 *  M_APM  -  mapm_fft.c
 
4
 *
 
5
 *  This FFT (Fast Fourier Transform) is from Takuya OOURA 
 
6
 *
 
7
 *  Copyright(C) 1996-1999 Takuya OOURA
 
8
 *  email: ooura@mmm.t.u-tokyo.ac.jp
 
9
 *
 
10
 *  See full FFT documentation below ...  (MCR)
 
11
 *
 
12
 *  This software is provided "as is" without express or implied warranty.
 
13
 */
 
14
 
 
15
/*
 
16
 *
 
17
 *      This file contains the FFT based FAST MULTIPLICATION function 
 
18
 *      as well as its support functions.
 
19
 *
 
20
 */
 
21
 
 
22
#include "pgAdmin3.h"
 
23
#include "pgscript/utilities/mapm-lib/m_apm_lc.h"
 
24
 
 
25
#ifndef MM_PI_2
 
26
#define MM_PI_2      1.570796326794896619231321691639751442098584699687
 
27
#endif
 
28
 
 
29
#ifndef WR5000       /* cos(MM_PI_2*0.5000) */
 
30
#define WR5000       0.707106781186547524400844362104849039284835937688
 
31
#endif
 
32
 
 
33
#ifndef RDFT_LOOP_DIV     /* control of the RDFT's speed & tolerance */
 
34
#define RDFT_LOOP_DIV 64
 
35
#endif
 
36
 
 
37
extern void   M_fast_mul_fft(UCHAR *, UCHAR *, UCHAR *, int);
 
38
 
 
39
extern void   M_rdft(int, int, double *);
 
40
extern void   M_bitrv2(int, double *);
 
41
extern void   M_cftfsub(int, double *);
 
42
extern void   M_cftbsub(int, double *);
 
43
extern void   M_rftfsub(int, double *);
 
44
extern void   M_rftbsub(int, double *);
 
45
extern void   M_cft1st(int, double *);
 
46
extern void   M_cftmdl(int, int, double *);
 
47
 
 
48
static double *M_aa_array, *M_bb_array;
 
49
static int    M_size = -1;
 
50
 
 
51
static char   *M_fft_error_msg = (char *)"\'M_fast_mul_fft\', Out of memory";
 
52
 
 
53
/****************************************************************************/
 
54
void    M_free_all_fft()
 
55
{
 
56
if (M_size > 0)
 
57
  {
 
58
   MAPM_FREE(M_aa_array);
 
59
   MAPM_FREE(M_bb_array);
 
60
   M_size = -1;
 
61
  }
 
62
}
 
63
/****************************************************************************/
 
64
/*
 
65
 *      multiply 'uu' by 'vv' with nbytes each
 
66
 *      yielding a 2*nbytes result in 'ww'. 
 
67
 *      each byte contains a base 100 'digit', 
 
68
 *      i.e.: range from 0-99.
 
69
 *
 
70
 *             MSB              LSB
 
71
 *
 
72
 *   uu,vv     [0] [1] [2] ... [N-1]
 
73
 *   ww        [0] [1] [2] ... [2N-1]
 
74
 */
 
75
 
 
76
void    M_fast_mul_fft(UCHAR *ww, UCHAR *uu, UCHAR *vv, int nbytes)
 
77
{
 
78
int             mflag, i, j, nn2, nn;
 
79
double          carry, nnr, dtemp, *a, *b;
 
80
UCHAR           *w0;
 
81
unsigned long   ul;
 
82
 
 
83
if (M_size < 0)                  /* if first time in, setup working arrays */
 
84
  {
 
85
   if (M_get_sizeof_int() == 2)  /* if still using 16 bit compilers */
 
86
     M_size = 516;
 
87
   else
 
88
     M_size = 8200;
 
89
 
 
90
   M_aa_array = (double *)MAPM_MALLOC(M_size * sizeof(double));
 
91
   M_bb_array = (double *)MAPM_MALLOC(M_size * sizeof(double));
 
92
   
 
93
   if ((M_aa_array == NULL) || (M_bb_array == NULL))
 
94
     {
 
95
      /* fatal, this does not return */
 
96
 
 
97
      M_apm_log_error_msg(M_APM_FATAL, M_fft_error_msg);
 
98
     }
 
99
  }
 
100
 
 
101
nn  = nbytes;
 
102
nn2 = nbytes >> 1;
 
103
 
 
104
if (nn > M_size)
 
105
  {
 
106
   mflag = TRUE;
 
107
 
 
108
   a = (double *)MAPM_MALLOC((nn + 8) * sizeof(double));
 
109
   b = (double *)MAPM_MALLOC((nn + 8) * sizeof(double));
 
110
   
 
111
   if ((a == NULL) || (b == NULL))
 
112
     {
 
113
      /* fatal, this does not return */
 
114
 
 
115
      M_apm_log_error_msg(M_APM_FATAL, M_fft_error_msg);
 
116
     }
 
117
  }
 
118
else
 
119
  {
 
120
   mflag = FALSE;
 
121
 
 
122
   a = M_aa_array;
 
123
   b = M_bb_array;
 
124
  }
 
125
 
 
126
/*
 
127
 *   convert normal base 100 MAPM numbers to base 10000
 
128
 *   for the FFT operation.
 
129
 */
 
130
 
 
131
i = 0;
 
132
for (j=0; j < nn2; j++)
 
133
  {
 
134
   a[j] = (double)((int)uu[i] * 100 + uu[i+1]);
 
135
   b[j] = (double)((int)vv[i] * 100 + vv[i+1]);
 
136
   i += 2;
 
137
  }
 
138
 
 
139
/* zero fill the second half of the arrays */
 
140
 
 
141
for (j=nn2; j < nn; j++)
 
142
  {
 
143
   a[j] = 0.0;
 
144
   b[j] = 0.0;
 
145
  }
 
146
 
 
147
/* perform the forward Fourier transforms for both numbers */
 
148
 
 
149
M_rdft(nn, 1, a);
 
150
M_rdft(nn, 1, b);
 
151
 
 
152
/* perform the convolution ... */
 
153
 
 
154
b[0] *= a[0];
 
155
b[1] *= a[1];
 
156
 
 
157
for (j=3; j <= nn; j += 2)
 
158
  {
 
159
   dtemp  = b[j-1];
 
160
   b[j-1] = dtemp * a[j-1] - b[j] * a[j];
 
161
   b[j]   = dtemp * a[j] + b[j] * a[j-1];
 
162
  }
 
163
 
 
164
/* perform the inverse transform on the result */
 
165
 
 
166
M_rdft(nn, -1, b);
 
167
 
 
168
/* perform a final pass to release all the carries */
 
169
/* we are still in base 10000 at this point        */
 
170
 
 
171
carry = 0.0;
 
172
j     = nn;
 
173
nnr   = 2.0 / (double)nn;
 
174
 
 
175
while (1)
 
176
  {
 
177
   dtemp = b[--j] * nnr + carry + 0.5;
 
178
   ul    = (unsigned long)(dtemp * 1.0E-4);
 
179
   carry = (double)ul;
 
180
   b[j]  = dtemp - carry * 10000.0;
 
181
 
 
182
   if (j == 0)
 
183
     break;
 
184
  }
 
185
 
 
186
/* copy result to our destination after converting back to base 100 */
 
187
 
 
188
w0 = ww;
 
189
M_get_div_rem((int)ul, w0, (w0 + 1));
 
190
 
 
191
for (j=0; j <= (nn - 2); j++)
 
192
  {
 
193
   w0 += 2;
 
194
   M_get_div_rem((int)b[j], w0, (w0 + 1));
 
195
  }
 
196
 
 
197
if (mflag)
 
198
  {
 
199
   MAPM_FREE(b);
 
200
   MAPM_FREE(a);
 
201
  }
 
202
}
 
203
/****************************************************************************/
 
204
 
 
205
/*
 
206
 *    The following info is from Takuya OOURA's documentation : 
 
207
 *
 
208
 *    NOTE : MAPM only uses the 'RDFT' function (as well as the 
 
209
 *           functions RDFT calls). All the code from here down 
 
210
 *           in this file is from Takuya OOURA. The only change I
 
211
 *           made was to add 'M_' in front of all the functions
 
212
 *           I used. This was to guard against any possible 
 
213
 *           name collisions in the future.
 
214
 *
 
215
 *    MCR  06 July 2000
 
216
 *
 
217
 *
 
218
 *    General Purpose FFT (Fast Fourier/Cosine/Sine Transform) Package
 
219
 *    
 
220
 *    Description:
 
221
 *        A package to calculate Discrete Fourier/Cosine/Sine Transforms of 
 
222
 *        1-dimensional sequences of length 2^N.
 
223
 *    
 
224
 *        fft4g_h.c  : FFT Package in C       - Simple Version I   (radix 4,2)
 
225
 *        
 
226
 *        rdft: Real Discrete Fourier Transform
 
227
 *    
 
228
 *    Method:
 
229
 *        -------- rdft --------
 
230
 *        A method with a following butterfly operation appended to "cdft".
 
231
 *        In forward transform :
 
232
 *            A[k] = sum_j=0^n-1 a[j]*W(n)^(j*k), 0<=k<=n/2, 
 
233
 *                W(n) = exp(2*pi*i/n), 
 
234
 *        this routine makes an array x[] :
 
235
 *            x[j] = a[2*j] + i*a[2*j+1], 0<=j<n/2
 
236
 *        and calls "cdft" of length n/2 :
 
237
 *            X[k] = sum_j=0^n/2-1 x[j] * W(n/2)^(j*k), 0<=k<n.
 
238
 *        The result A[k] are :
 
239
 *            A[k]     = X[k]     - (1+i*W(n)^k)/2 * (X[k]-conjg(X[n/2-k])),
 
240
 *            A[n/2-k] = X[n/2-k] + 
 
241
 *                            conjg((1+i*W(n)^k)/2 * (X[k]-conjg(X[n/2-k]))),
 
242
 *                0<=k<=n/2
 
243
 *            (notes: conjg() is a complex conjugate, X[n/2]=X[0]).
 
244
 *        ----------------------
 
245
 *    
 
246
 *    Reference:
 
247
 *        * Masatake MORI, Makoto NATORI, Tatuo TORII: Suchikeisan, 
 
248
 *          Iwanamikouzajyouhoukagaku18, Iwanami, 1982 (Japanese)
 
249
 *        * Henri J. Nussbaumer: Fast Fourier Transform and Convolution 
 
250
 *          Algorithms, Springer Verlag, 1982
 
251
 *        * C. S. Burrus, Notes on the FFT (with large FFT paper list)
 
252
 *          http://www-dsp.rice.edu/research/fft/fftnote.asc
 
253
 *    
 
254
 *    Copyright:
 
255
 *        Copyright(C) 1996-1999 Takuya OOURA
 
256
 *        email: ooura@mmm.t.u-tokyo.ac.jp
 
257
 *        download: http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html
 
258
 *        You may use, copy, modify this code for any purpose and 
 
259
 *        without fee. You may distribute this ORIGINAL package.
 
260
 */
 
261
 
 
262
/*
 
263
 
 
264
functions
 
265
    rdft: Real Discrete Fourier Transform
 
266
 
 
267
function prototypes
 
268
    void rdft(int, int, double *);
 
269
 
 
270
-------- Real DFT / Inverse of Real DFT --------
 
271
    [definition]
 
272
        <case1> RDFT
 
273
            R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
 
274
            I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
 
275
        <case2> IRDFT (excluding scale)
 
276
            a[k] = (R[0] + R[n/2]*cos(pi*k))/2 + 
 
277
                   sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) + 
 
278
                   sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
 
279
    [usage]
 
280
        <case1>
 
281
            rdft(n, 1, a);
 
282
        <case2>
 
283
            rdft(n, -1, a);
 
284
    [parameters]
 
285
        n              :data length (int)
 
286
                        n >= 2, n = power of 2
 
287
        a[0...n-1]     :input/output data (double *)
 
288
                        <case1>
 
289
                            output data
 
290
                                a[2*k] = R[k], 0<=k<n/2
 
291
                                a[2*k+1] = I[k], 0<k<n/2
 
292
                                a[1] = R[n/2]
 
293
                        <case2>
 
294
                            input data
 
295
                                a[2*j] = R[j], 0<=j<n/2
 
296
                                a[2*j+1] = I[j], 0<j<n/2
 
297
                                a[1] = R[n/2]
 
298
    [remark]
 
299
        Inverse of 
 
300
            rdft(n, 1, a);
 
301
        is 
 
302
            rdft(n, -1, a);
 
303
            for (j = 0; j <= n - 1; j++) {
 
304
                a[j] *= 2.0 / n;
 
305
            }
 
306
*/
 
307
 
 
308
 
 
309
void    M_rdft(int n, int isgn, double *a)
 
310
{
 
311
    double xi;
 
312
 
 
313
    if (isgn >= 0) {
 
314
        if (n > 4) {
 
315
            M_bitrv2(n, a);
 
316
            M_cftfsub(n, a);
 
317
            M_rftfsub(n, a);
 
318
        } else if (n == 4) {
 
319
            M_cftfsub(n, a);
 
320
        }
 
321
        xi = a[0] - a[1];
 
322
        a[0] += a[1];
 
323
        a[1] = xi;
 
324
    } else {
 
325
        a[1] = 0.5 * (a[0] - a[1]);
 
326
        a[0] -= a[1];
 
327
        if (n > 4) {
 
328
            M_rftbsub(n, a);
 
329
            M_bitrv2(n, a);
 
330
            M_cftbsub(n, a);
 
331
        } else if (n == 4) {
 
332
            M_cftfsub(n, a);
 
333
        }
 
334
    }
 
335
}
 
336
 
 
337
 
 
338
 
 
339
void    M_bitrv2(int n, double *a)
 
340
{
 
341
    int j0, k0, j1, k1, l, m, i, j, k;
 
342
    double xr, xi, yr, yi;
 
343
    
 
344
    l = n >> 2;
 
345
    m = 2;
 
346
    while (m < l) {
 
347
        l >>= 1;
 
348
        m <<= 1;
 
349
    }
 
350
    if (m == l) {
 
351
        j0 = 0;
 
352
        for (k0 = 0; k0 < m; k0 += 2) {
 
353
            k = k0;
 
354
            for (j = j0; j < j0 + k0; j += 2) {
 
355
                xr = a[j];
 
356
                xi = a[j + 1];
 
357
                yr = a[k];
 
358
                yi = a[k + 1];
 
359
                a[j] = yr;
 
360
                a[j + 1] = yi;
 
361
                a[k] = xr;
 
362
                a[k + 1] = xi;
 
363
                j1 = j + m;
 
364
                k1 = k + 2 * m;
 
365
                xr = a[j1];
 
366
                xi = a[j1 + 1];
 
367
                yr = a[k1];
 
368
                yi = a[k1 + 1];
 
369
                a[j1] = yr;
 
370
                a[j1 + 1] = yi;
 
371
                a[k1] = xr;
 
372
                a[k1 + 1] = xi;
 
373
                j1 += m;
 
374
                k1 -= m;
 
375
                xr = a[j1];
 
376
                xi = a[j1 + 1];
 
377
                yr = a[k1];
 
378
                yi = a[k1 + 1];
 
379
                a[j1] = yr;
 
380
                a[j1 + 1] = yi;
 
381
                a[k1] = xr;
 
382
                a[k1 + 1] = xi;
 
383
                j1 += m;
 
384
                k1 += 2 * m;
 
385
                xr = a[j1];
 
386
                xi = a[j1 + 1];
 
387
                yr = a[k1];
 
388
                yi = a[k1 + 1];
 
389
                a[j1] = yr;
 
390
                a[j1 + 1] = yi;
 
391
                a[k1] = xr;
 
392
                a[k1 + 1] = xi;
 
393
                for (i = n >> 1; i > (k ^= i); i >>= 1);
 
394
            }
 
395
            j1 = j0 + k0 + m;
 
396
            k1 = j1 + m;
 
397
            xr = a[j1];
 
398
            xi = a[j1 + 1];
 
399
            yr = a[k1];
 
400
            yi = a[k1 + 1];
 
401
            a[j1] = yr;
 
402
            a[j1 + 1] = yi;
 
403
            a[k1] = xr;
 
404
            a[k1 + 1] = xi;
 
405
            for (i = n >> 1; i > (j0 ^= i); i >>= 1);
 
406
        }
 
407
    } else {
 
408
        j0 = 0;
 
409
        for (k0 = 2; k0 < m; k0 += 2) {
 
410
            for (i = n >> 1; i > (j0 ^= i); i >>= 1);
 
411
            k = k0;
 
412
            for (j = j0; j < j0 + k0; j += 2) {
 
413
                xr = a[j];
 
414
                xi = a[j + 1];
 
415
                yr = a[k];
 
416
                yi = a[k + 1];
 
417
                a[j] = yr;
 
418
                a[j + 1] = yi;
 
419
                a[k] = xr;
 
420
                a[k + 1] = xi;
 
421
                j1 = j + m;
 
422
                k1 = k + m;
 
423
                xr = a[j1];
 
424
                xi = a[j1 + 1];
 
425
                yr = a[k1];
 
426
                yi = a[k1 + 1];
 
427
                a[j1] = yr;
 
428
                a[j1 + 1] = yi;
 
429
                a[k1] = xr;
 
430
                a[k1 + 1] = xi;
 
431
                for (i = n >> 1; i > (k ^= i); i >>= 1);
 
432
            }
 
433
        }
 
434
    }
 
435
}
 
436
 
 
437
 
 
438
 
 
439
void    M_cftfsub(int n, double *a)
 
440
{
 
441
    int j, j1, j2, j3, l;
 
442
    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
 
443
    
 
444
    l = 2;
 
445
    if (n > 8) {
 
446
        M_cft1st(n, a);
 
447
        l = 8;
 
448
        while ((l << 2) < n) {
 
449
            M_cftmdl(n, l, a);
 
450
            l <<= 2;
 
451
        }
 
452
    }
 
453
    if ((l << 2) == n) {
 
454
        for (j = 0; j < l; j += 2) {
 
455
            j1 = j + l;
 
456
            j2 = j1 + l;
 
457
            j3 = j2 + l;
 
458
            x0r = a[j] + a[j1];
 
459
            x0i = a[j + 1] + a[j1 + 1];
 
460
            x1r = a[j] - a[j1];
 
461
            x1i = a[j + 1] - a[j1 + 1];
 
462
            x2r = a[j2] + a[j3];
 
463
            x2i = a[j2 + 1] + a[j3 + 1];
 
464
            x3r = a[j2] - a[j3];
 
465
            x3i = a[j2 + 1] - a[j3 + 1];
 
466
            a[j] = x0r + x2r;
 
467
            a[j + 1] = x0i + x2i;
 
468
            a[j2] = x0r - x2r;
 
469
            a[j2 + 1] = x0i - x2i;
 
470
            a[j1] = x1r - x3i;
 
471
            a[j1 + 1] = x1i + x3r;
 
472
            a[j3] = x1r + x3i;
 
473
            a[j3 + 1] = x1i - x3r;
 
474
        }
 
475
    } else {
 
476
        for (j = 0; j < l; j += 2) {
 
477
            j1 = j + l;
 
478
            x0r = a[j] - a[j1];
 
479
            x0i = a[j + 1] - a[j1 + 1];
 
480
            a[j] += a[j1];
 
481
            a[j + 1] += a[j1 + 1];
 
482
            a[j1] = x0r;
 
483
            a[j1 + 1] = x0i;
 
484
        }
 
485
    }
 
486
}
 
487
 
 
488
 
 
489
 
 
490
void    M_cftbsub(int n, double *a)
 
491
{
 
492
    int j, j1, j2, j3, l;
 
493
    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
 
494
    
 
495
    l = 2;
 
496
    if (n > 8) {
 
497
        M_cft1st(n, a);
 
498
        l = 8;
 
499
        while ((l << 2) < n) {
 
500
            M_cftmdl(n, l, a);
 
501
            l <<= 2;
 
502
        }
 
503
    }
 
504
    if ((l << 2) == n) {
 
505
        for (j = 0; j < l; j += 2) {
 
506
            j1 = j + l;
 
507
            j2 = j1 + l;
 
508
            j3 = j2 + l;
 
509
            x0r = a[j] + a[j1];
 
510
            x0i = -a[j + 1] - a[j1 + 1];
 
511
            x1r = a[j] - a[j1];
 
512
            x1i = -a[j + 1] + a[j1 + 1];
 
513
            x2r = a[j2] + a[j3];
 
514
            x2i = a[j2 + 1] + a[j3 + 1];
 
515
            x3r = a[j2] - a[j3];
 
516
            x3i = a[j2 + 1] - a[j3 + 1];
 
517
            a[j] = x0r + x2r;
 
518
            a[j + 1] = x0i - x2i;
 
519
            a[j2] = x0r - x2r;
 
520
            a[j2 + 1] = x0i + x2i;
 
521
            a[j1] = x1r - x3i;
 
522
            a[j1 + 1] = x1i - x3r;
 
523
            a[j3] = x1r + x3i;
 
524
            a[j3 + 1] = x1i + x3r;
 
525
        }
 
526
    } else {
 
527
        for (j = 0; j < l; j += 2) {
 
528
            j1 = j + l;
 
529
            x0r = a[j] - a[j1];
 
530
            x0i = -a[j + 1] + a[j1 + 1];
 
531
            a[j] += a[j1];
 
532
            a[j + 1] = -a[j + 1] - a[j1 + 1];
 
533
            a[j1] = x0r;
 
534
            a[j1 + 1] = x0i;
 
535
        }
 
536
    }
 
537
}
 
538
 
 
539
 
 
540
 
 
541
void    M_cft1st(int n, double *a)
 
542
{
 
543
    int j, kj, kr;
 
544
    double ew, wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
 
545
    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
 
546
    
 
547
    x0r = a[0] + a[2];
 
548
    x0i = a[1] + a[3];
 
549
    x1r = a[0] - a[2];
 
550
    x1i = a[1] - a[3];
 
551
    x2r = a[4] + a[6];
 
552
    x2i = a[5] + a[7];
 
553
    x3r = a[4] - a[6];
 
554
    x3i = a[5] - a[7];
 
555
    a[0] = x0r + x2r;
 
556
    a[1] = x0i + x2i;
 
557
    a[4] = x0r - x2r;
 
558
    a[5] = x0i - x2i;
 
559
    a[2] = x1r - x3i;
 
560
    a[3] = x1i + x3r;
 
561
    a[6] = x1r + x3i;
 
562
    a[7] = x1i - x3r;
 
563
    wn4r = WR5000;
 
564
    x0r = a[8] + a[10];
 
565
    x0i = a[9] + a[11];
 
566
    x1r = a[8] - a[10];
 
567
    x1i = a[9] - a[11];
 
568
    x2r = a[12] + a[14];
 
569
    x2i = a[13] + a[15];
 
570
    x3r = a[12] - a[14];
 
571
    x3i = a[13] - a[15];
 
572
    a[8] = x0r + x2r;
 
573
    a[9] = x0i + x2i;
 
574
    a[12] = x2i - x0i;
 
575
    a[13] = x0r - x2r;
 
576
    x0r = x1r - x3i;
 
577
    x0i = x1i + x3r;
 
578
    a[10] = wn4r * (x0r - x0i);
 
579
    a[11] = wn4r * (x0r + x0i);
 
580
    x0r = x3i + x1r;
 
581
    x0i = x3r - x1i;
 
582
    a[14] = wn4r * (x0i - x0r);
 
583
    a[15] = wn4r * (x0i + x0r);
 
584
    ew = MM_PI_2 / n;
 
585
    kr = 0;
 
586
    for (j = 16; j < n; j += 16) {
 
587
        for (kj = n >> 2; kj > (kr ^= kj); kj >>= 1);
 
588
        wk1r = cos(ew * kr);
 
589
        wk1i = sin(ew * kr);
 
590
        wk2r = 1 - 2 * wk1i * wk1i;
 
591
        wk2i = 2 * wk1i * wk1r;
 
592
        wk3r = wk1r - 2 * wk2i * wk1i;
 
593
        wk3i = 2 * wk2i * wk1r - wk1i;
 
594
        x0r = a[j] + a[j + 2];
 
595
        x0i = a[j + 1] + a[j + 3];
 
596
        x1r = a[j] - a[j + 2];
 
597
        x1i = a[j + 1] - a[j + 3];
 
598
        x2r = a[j + 4] + a[j + 6];
 
599
        x2i = a[j + 5] + a[j + 7];
 
600
        x3r = a[j + 4] - a[j + 6];
 
601
        x3i = a[j + 5] - a[j + 7];
 
602
        a[j] = x0r + x2r;
 
603
        a[j + 1] = x0i + x2i;
 
604
        x0r -= x2r;
 
605
        x0i -= x2i;
 
606
        a[j + 4] = wk2r * x0r - wk2i * x0i;
 
607
        a[j + 5] = wk2r * x0i + wk2i * x0r;
 
608
        x0r = x1r - x3i;
 
609
        x0i = x1i + x3r;
 
610
        a[j + 2] = wk1r * x0r - wk1i * x0i;
 
611
        a[j + 3] = wk1r * x0i + wk1i * x0r;
 
612
        x0r = x1r + x3i;
 
613
        x0i = x1i - x3r;
 
614
        a[j + 6] = wk3r * x0r - wk3i * x0i;
 
615
        a[j + 7] = wk3r * x0i + wk3i * x0r;
 
616
        x0r = wn4r * (wk1r - wk1i);
 
617
        wk1i = wn4r * (wk1r + wk1i);
 
618
        wk1r = x0r;
 
619
        wk3r = wk1r - 2 * wk2r * wk1i;
 
620
        wk3i = 2 * wk2r * wk1r - wk1i;
 
621
        x0r = a[j + 8] + a[j + 10];
 
622
        x0i = a[j + 9] + a[j + 11];
 
623
        x1r = a[j + 8] - a[j + 10];
 
624
        x1i = a[j + 9] - a[j + 11];
 
625
        x2r = a[j + 12] + a[j + 14];
 
626
        x2i = a[j + 13] + a[j + 15];
 
627
        x3r = a[j + 12] - a[j + 14];
 
628
        x3i = a[j + 13] - a[j + 15];
 
629
        a[j + 8] = x0r + x2r;
 
630
        a[j + 9] = x0i + x2i;
 
631
        x0r -= x2r;
 
632
        x0i -= x2i;
 
633
        a[j + 12] = -wk2i * x0r - wk2r * x0i;
 
634
        a[j + 13] = -wk2i * x0i + wk2r * x0r;
 
635
        x0r = x1r - x3i;
 
636
        x0i = x1i + x3r;
 
637
        a[j + 10] = wk1r * x0r - wk1i * x0i;
 
638
        a[j + 11] = wk1r * x0i + wk1i * x0r;
 
639
        x0r = x1r + x3i;
 
640
        x0i = x1i - x3r;
 
641
        a[j + 14] = wk3r * x0r - wk3i * x0i;
 
642
        a[j + 15] = wk3r * x0i + wk3i * x0r;
 
643
    }
 
644
}
 
645
 
 
646
 
 
647
 
 
648
void    M_cftmdl(int n, int l, double *a)
 
649
{
 
650
    int j, j1, j2, j3, k, kj, kr, m, m2;
 
651
    double ew, wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
 
652
    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
 
653
    
 
654
    m = l << 2;
 
655
    for (j = 0; j < l; j += 2) {
 
656
        j1 = j + l;
 
657
        j2 = j1 + l;
 
658
        j3 = j2 + l;
 
659
        x0r = a[j] + a[j1];
 
660
        x0i = a[j + 1] + a[j1 + 1];
 
661
        x1r = a[j] - a[j1];
 
662
        x1i = a[j + 1] - a[j1 + 1];
 
663
        x2r = a[j2] + a[j3];
 
664
        x2i = a[j2 + 1] + a[j3 + 1];
 
665
        x3r = a[j2] - a[j3];
 
666
        x3i = a[j2 + 1] - a[j3 + 1];
 
667
        a[j] = x0r + x2r;
 
668
        a[j + 1] = x0i + x2i;
 
669
        a[j2] = x0r - x2r;
 
670
        a[j2 + 1] = x0i - x2i;
 
671
        a[j1] = x1r - x3i;
 
672
        a[j1 + 1] = x1i + x3r;
 
673
        a[j3] = x1r + x3i;
 
674
        a[j3 + 1] = x1i - x3r;
 
675
    }
 
676
    wn4r = WR5000;
 
677
    for (j = m; j < l + m; j += 2) {
 
678
        j1 = j + l;
 
679
        j2 = j1 + l;
 
680
        j3 = j2 + l;
 
681
        x0r = a[j] + a[j1];
 
682
        x0i = a[j + 1] + a[j1 + 1];
 
683
        x1r = a[j] - a[j1];
 
684
        x1i = a[j + 1] - a[j1 + 1];
 
685
        x2r = a[j2] + a[j3];
 
686
        x2i = a[j2 + 1] + a[j3 + 1];
 
687
        x3r = a[j2] - a[j3];
 
688
        x3i = a[j2 + 1] - a[j3 + 1];
 
689
        a[j] = x0r + x2r;
 
690
        a[j + 1] = x0i + x2i;
 
691
        a[j2] = x2i - x0i;
 
692
        a[j2 + 1] = x0r - x2r;
 
693
        x0r = x1r - x3i;
 
694
        x0i = x1i + x3r;
 
695
        a[j1] = wn4r * (x0r - x0i);
 
696
        a[j1 + 1] = wn4r * (x0r + x0i);
 
697
        x0r = x3i + x1r;
 
698
        x0i = x3r - x1i;
 
699
        a[j3] = wn4r * (x0i - x0r);
 
700
        a[j3 + 1] = wn4r * (x0i + x0r);
 
701
    }
 
702
    ew = MM_PI_2 / n;
 
703
    kr = 0;
 
704
    m2 = 2 * m;
 
705
    for (k = m2; k < n; k += m2) {
 
706
        for (kj = n >> 2; kj > (kr ^= kj); kj >>= 1);
 
707
        wk1r = cos(ew * kr);
 
708
        wk1i = sin(ew * kr);
 
709
        wk2r = 1 - 2 * wk1i * wk1i;
 
710
        wk2i = 2 * wk1i * wk1r;
 
711
        wk3r = wk1r - 2 * wk2i * wk1i;
 
712
        wk3i = 2 * wk2i * wk1r - wk1i;
 
713
        for (j = k; j < l + k; j += 2) {
 
714
            j1 = j + l;
 
715
            j2 = j1 + l;
 
716
            j3 = j2 + l;
 
717
            x0r = a[j] + a[j1];
 
718
            x0i = a[j + 1] + a[j1 + 1];
 
719
            x1r = a[j] - a[j1];
 
720
            x1i = a[j + 1] - a[j1 + 1];
 
721
            x2r = a[j2] + a[j3];
 
722
            x2i = a[j2 + 1] + a[j3 + 1];
 
723
            x3r = a[j2] - a[j3];
 
724
            x3i = a[j2 + 1] - a[j3 + 1];
 
725
            a[j] = x0r + x2r;
 
726
            a[j + 1] = x0i + x2i;
 
727
            x0r -= x2r;
 
728
            x0i -= x2i;
 
729
            a[j2] = wk2r * x0r - wk2i * x0i;
 
730
            a[j2 + 1] = wk2r * x0i + wk2i * x0r;
 
731
            x0r = x1r - x3i;
 
732
            x0i = x1i + x3r;
 
733
            a[j1] = wk1r * x0r - wk1i * x0i;
 
734
            a[j1 + 1] = wk1r * x0i + wk1i * x0r;
 
735
            x0r = x1r + x3i;
 
736
            x0i = x1i - x3r;
 
737
            a[j3] = wk3r * x0r - wk3i * x0i;
 
738
            a[j3 + 1] = wk3r * x0i + wk3i * x0r;
 
739
        }
 
740
        x0r = wn4r * (wk1r - wk1i);
 
741
        wk1i = wn4r * (wk1r + wk1i);
 
742
        wk1r = x0r;
 
743
        wk3r = wk1r - 2 * wk2r * wk1i;
 
744
        wk3i = 2 * wk2r * wk1r - wk1i;
 
745
        for (j = k + m; j < l + (k + m); j += 2) {
 
746
            j1 = j + l;
 
747
            j2 = j1 + l;
 
748
            j3 = j2 + l;
 
749
            x0r = a[j] + a[j1];
 
750
            x0i = a[j + 1] + a[j1 + 1];
 
751
            x1r = a[j] - a[j1];
 
752
            x1i = a[j + 1] - a[j1 + 1];
 
753
            x2r = a[j2] + a[j3];
 
754
            x2i = a[j2 + 1] + a[j3 + 1];
 
755
            x3r = a[j2] - a[j3];
 
756
            x3i = a[j2 + 1] - a[j3 + 1];
 
757
            a[j] = x0r + x2r;
 
758
            a[j + 1] = x0i + x2i;
 
759
            x0r -= x2r;
 
760
            x0i -= x2i;
 
761
            a[j2] = -wk2i * x0r - wk2r * x0i;
 
762
            a[j2 + 1] = -wk2i * x0i + wk2r * x0r;
 
763
            x0r = x1r - x3i;
 
764
            x0i = x1i + x3r;
 
765
            a[j1] = wk1r * x0r - wk1i * x0i;
 
766
            a[j1 + 1] = wk1r * x0i + wk1i * x0r;
 
767
            x0r = x1r + x3i;
 
768
            x0i = x1i - x3r;
 
769
            a[j3] = wk3r * x0r - wk3i * x0i;
 
770
            a[j3 + 1] = wk3r * x0i + wk3i * x0r;
 
771
        }
 
772
    }
 
773
}
 
774
 
 
775
 
 
776
 
 
777
void    M_rftfsub(int n, double *a)
 
778
{
 
779
    int i, i0, j, k;
 
780
    double ec, w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi;
 
781
    
 
782
    ec = 2 * MM_PI_2 / n;
 
783
    wkr = 0;
 
784
    wki = 0;
 
785
    wdi = cos(ec);
 
786
    wdr = sin(ec);
 
787
    wdi *= wdr;
 
788
    wdr *= wdr;
 
789
    w1r = 1 - 2 * wdr;
 
790
    w1i = 2 * wdi;
 
791
    ss = 2 * w1i;
 
792
    i = n >> 1;
 
793
    while (1) {
 
794
        i0 = i - 4 * RDFT_LOOP_DIV;
 
795
        if (i0 < 4) {
 
796
            i0 = 4;
 
797
        }
 
798
        for (j = i - 4; j >= i0; j -= 4) {
 
799
            k = n - j;
 
800
            xr = a[j + 2] - a[k - 2];
 
801
            xi = a[j + 3] + a[k - 1];
 
802
            yr = wdr * xr - wdi * xi;
 
803
            yi = wdr * xi + wdi * xr;
 
804
            a[j + 2] -= yr;
 
805
            a[j + 3] -= yi;
 
806
            a[k - 2] += yr;
 
807
            a[k - 1] -= yi;
 
808
            wkr += ss * wdi;
 
809
            wki += ss * (0.5 - wdr);
 
810
            xr = a[j] - a[k];
 
811
            xi = a[j + 1] + a[k + 1];
 
812
            yr = wkr * xr - wki * xi;
 
813
            yi = wkr * xi + wki * xr;
 
814
            a[j] -= yr;
 
815
            a[j + 1] -= yi;
 
816
            a[k] += yr;
 
817
            a[k + 1] -= yi;
 
818
            wdr += ss * wki;
 
819
            wdi += ss * (0.5 - wkr);
 
820
        }
 
821
        if (i0 == 4) {
 
822
            break;
 
823
        }
 
824
        wkr = 0.5 * sin(ec * i0);
 
825
        wki = 0.5 * cos(ec * i0);
 
826
        wdr = 0.5 - (wkr * w1r - wki * w1i);
 
827
        wdi = wkr * w1i + wki * w1r;
 
828
        wkr = 0.5 - wkr;
 
829
        i = i0;
 
830
    }
 
831
    xr = a[2] - a[n - 2];
 
832
    xi = a[3] + a[n - 1];
 
833
    yr = wdr * xr - wdi * xi;
 
834
    yi = wdr * xi + wdi * xr;
 
835
    a[2] -= yr;
 
836
    a[3] -= yi;
 
837
    a[n - 2] += yr;
 
838
    a[n - 1] -= yi;
 
839
}
 
840
 
 
841
 
 
842
 
 
843
void    M_rftbsub(int n, double *a)
 
844
{
 
845
    int i, i0, j, k;
 
846
    double ec, w1r, w1i, wkr, wki, wdr, wdi, ss, xr, xi, yr, yi;
 
847
    
 
848
    ec = 2 * MM_PI_2 / n;
 
849
    wkr = 0;
 
850
    wki = 0;
 
851
    wdi = cos(ec);
 
852
    wdr = sin(ec);
 
853
    wdi *= wdr;
 
854
    wdr *= wdr;
 
855
    w1r = 1 - 2 * wdr;
 
856
    w1i = 2 * wdi;
 
857
    ss = 2 * w1i;
 
858
    i = n >> 1;
 
859
    a[i + 1] = -a[i + 1];
 
860
    while (1) {
 
861
        i0 = i - 4 * RDFT_LOOP_DIV;
 
862
        if (i0 < 4) {
 
863
            i0 = 4;
 
864
        }
 
865
        for (j = i - 4; j >= i0; j -= 4) {
 
866
            k = n - j;
 
867
            xr = a[j + 2] - a[k - 2];
 
868
            xi = a[j + 3] + a[k - 1];
 
869
            yr = wdr * xr + wdi * xi;
 
870
            yi = wdr * xi - wdi * xr;
 
871
            a[j + 2] -= yr;
 
872
            a[j + 3] = yi - a[j + 3];
 
873
            a[k - 2] += yr;
 
874
            a[k - 1] = yi - a[k - 1];
 
875
            wkr += ss * wdi;
 
876
            wki += ss * (0.5 - wdr);
 
877
            xr = a[j] - a[k];
 
878
            xi = a[j + 1] + a[k + 1];
 
879
            yr = wkr * xr + wki * xi;
 
880
            yi = wkr * xi - wki * xr;
 
881
            a[j] -= yr;
 
882
            a[j + 1] = yi - a[j + 1];
 
883
            a[k] += yr;
 
884
            a[k + 1] = yi - a[k + 1];
 
885
            wdr += ss * wki;
 
886
            wdi += ss * (0.5 - wkr);
 
887
        }
 
888
        if (i0 == 4) {
 
889
            break;
 
890
        }
 
891
        wkr = 0.5 * sin(ec * i0);
 
892
        wki = 0.5 * cos(ec * i0);
 
893
        wdr = 0.5 - (wkr * w1r - wki * w1i);
 
894
        wdi = wkr * w1i + wki * w1r;
 
895
        wkr = 0.5 - wkr;
 
896
        i = i0;
 
897
    }
 
898
    xr = a[2] - a[n - 2];
 
899
    xi = a[3] + a[n - 1];
 
900
    yr = wdr * xr + wdi * xi;
 
901
    yi = wdr * xi - wdi * xr;
 
902
    a[2] -= yr;
 
903
    a[3] = yi - a[3];
 
904
    a[n - 2] += yr;
 
905
    a[n - 1] = yi - a[n - 1];
 
906
    a[1] = -a[1];
 
907
}
 
908