← Back to branch summary

~ubuntu-branches/ubuntu/quantal/linpsk/quantal

« back to all changes in this revision

Viewing changes to linpsk/fft.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Bruce Walker
  • Date: 2002-02-06 11:43:38 UTC
  • Revision ID: james.westby@ubuntu.com-20020206114338-xqmjmhh01lpjm0g4
Tags: upstream-0.6.2
Import upstream version 0.6.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
linpsk/fft.cpp
 
1
// fft.cpp: implementation of the Cfft class.
 
2
//  This is a slightly modified version of Takuya OOURA's
 
3
//     original radix 4 FFT package.
 
4
//Copyright(C) 1996-1998 Takuya OOURA
 
5
//    (email: ooura@mmm.t.u-tokyo.ac.jp).
 
6
//////////////////////////////////////////////////////////////////////
 
7
 
 
8
 
 
9
#include <math.h>
 
10
#include "fft.h"
 
11
 
 
12
//////////////////////////////////////////////////////////////////////
 
13
// Local Defines
 
14
//////////////////////////////////////////////////////////////////////
 
15
#define SQRT_FFT_SIZE 46//sqrt(2048)
 
16
#define K_2PI 6.283185307
 
17
 
 
18
//////////////////////////////////////////////////////////////////////
 
19
// A pure input sin wave ... Asin(wt)... will produce an fft output 
 
20
//   peak of (N*A/4)^2  where N is FFT_SIZE.
 
21
// To convert to a Power dB range:
 
22
//   PdBmax = 10*log10( (N*A/4)^2 + K_C ) + K_B
 
23
//   PdBmin = 10*log10( 0 + K_C ) + K_B
 
24
//  if (N*A/4)^2 >> K_C 
 
25
//  Then K_B = PdBmax - 20*log10( N*A/4 )
 
26
//       K_C = 10 ^ ( (PdBmin-K_B)/10 )
 
27
//  for power range of 0 to 100 dB with input(A) of 32767 and N=4096
 
28
//                      K_B = -50.51473275  and K_C = 1.1258312e5
 
29
// To eliminate the multiply by 10, divide by 10 so for an output
 
30
//              range of 0 to 100dB the stored value is 0.0 to 10.0
 
31
//   so final constant K_B = -5.051473275
 
32
#define K_B (-5.051473275)
 
33
#define K_C (1.1258312e5)
 
34
//////////////////////////////////////////////////////////////////////
 
35
// Construction/Destruction
 
36
//////////////////////////////////////////////////////////////////////
 
37
Cfft::Cfft()
 
38
{
 
39
        InBuf   = new double[FFT_SIZE];
 
40
        WindowTbl = new double[FFT_SIZE];
 
41
        SinCosTbl = new double[FFT_SIZE/2];
 
42
        WorkArea = new int[SQRT_FFT_SIZE+2];
 
43
        pFFTAveBuf = new double[FFT_SIZE];      //Even index's hold average
 
44
        WorkArea[0] = 0;
 
45
        makewt(FFT_SIZE/4, WorkArea, SinCosTbl);
 
46
    makect(FFT_SIZE/4, WorkArea, SinCosTbl + WorkArea[0]);
 
47
        for(int i=0; i<FFT_SIZE; i++)
 
48
        {
 
49
                pFFTAveBuf[i] = 0.0;
 
50
// Pick a data windowing function:
 
51
//              WindowTbl[i] = 1.0;                                                                             //rectangle
 
52
//              WindowTbl[i] = .54 - .46*cos( (K_2PI*i)/(FFT_SIZE-1) ); //Hamming
 
53
                WindowTbl[i] = (.5 - .5*cos( (K_2PI*i)/(FFT_SIZE-1) )); //Hanning
 
54
        }
 
55
        m_AveSize = 1;
 
56
        m_AverageOn = false;
 
57
}
 
58
 
 
59
Cfft::~Cfft()
 
60
{                                                       // free all resources
 
61
        if(WorkArea)
 
62
        {
 
63
                delete WorkArea;
 
64
                WorkArea = NULL;
 
65
        }
 
66
        if(SinCosTbl)
 
67
        {
 
68
                delete SinCosTbl;
 
69
                SinCosTbl = NULL;
 
70
        }
 
71
        if(WindowTbl)
 
72
        {
 
73
                delete WindowTbl;
 
74
                WindowTbl = NULL;
 
75
        }
 
76
        if(pFFTAveBuf)
 
77
        {
 
78
                delete pFFTAveBuf;
 
79
                pFFTAveBuf = NULL;
 
80
        }
 
81
        if(InBuf)
 
82
                {
 
83
                        delete InBuf;
 
84
                        InBuf = NULL;
 
85
                }       
 
86
 
 
87
}
 
88
 
 
89
//////////////////////////////////////////////////////////////////////
 
90
// "InBuf[]" is first multiplied by a window function and then
 
91
//  calculates an "FFT_SIZE" point FFT on "InBuf[]".
 
92
//  The result is converted to dB, and the range "start" to "stop"
 
93
//  is multiplied by "gain" and copied into "OutBuf[]".
 
94
//  If "Ave" is > 1, a LP smoothing filter 
 
95
//  is calculated on the output.
 
96
//     Execution times for 4096 point fft:
 
97
//                       ~820nSec/sample with 266MHz PII
 
98
// The function returns true if the input is overloaded
 
99
//////////////////////////////////////////////////////////////////////
 
100
bool Cfft::CalcFFT(double * InPut, int start, int stop,double gain, int Ave,int* OutBuf)
 
101
{
 
102
int i;
 
103
bool Ovrld = false;
 
104
        if( Ave>1 )
 
105
                m_AverageOn = true;
 
106
        else
 
107
                m_AverageOn = false;
 
108
        m_AveSize = Ave;
 
109
        for(i=0; i<FFT_SIZE; i++)
 
110
        {
 
111
                if( InPut[i] > 16384.0 )        //flag overload if within 3dB of max
 
112
                        Ovrld = true;
 
113
                InBuf[i] = WindowTbl[i] * InPut[i];             //window the data
 
114
        }
 
115
//Calculate the FFT
 
116
        bitrv2(FFT_SIZE, WorkArea + 2, InBuf);
 
117
        cftfsub(FFT_SIZE, InBuf, SinCosTbl);
 
118
        rftfsub(FFT_SIZE, InBuf, WorkArea[1], SinCosTbl + WorkArea[0]);
 
119
        for( i=start; i<=stop; i++ )            //copy and scale into OutBuf[]
 
120
                OutBuf[i] = (int)(gain*pFFTAveBuf[i<<1]);
 
121
        return Ovrld;
 
122
}
 
123
 
 
124
void Cfft::rftfsub(int n, double *a, int nc, double *c)
 
125
{
 
126
    int j, k, kk, ks, m;
 
127
    double wkr, wki, xr, xi, yr, yi;
 
128
    
 
129
    ks = (nc << 2) / n;
 
130
    kk = 0;
 
131
    m = n >> 1;
 
132
        if(m_AverageOn)
 
133
        {
 
134
                for (k = 2; k <= m; k += 2 ) 
 
135
                {
 
136
                        j = n - k;
 
137
                        kk += ks;
 
138
                        wkr = 0.5 - c[nc - kk];
 
139
                        wki = c[kk];
 
140
                        xr = a[k] - a[j];
 
141
                        xi = a[k + 1] + a[j + 1];
 
142
                        yr = wkr * xr - wki * xi;
 
143
                        yi = wkr * xi + wki * xr;
 
144
                        a[k] -= yr;
 
145
                        xi = a[k]*a[k];
 
146
                        a[k+1] -= yi;
 
147
                        xi += ( a[k+1]*a[k+1]);
 
148
                        a[j] += yr;
 
149
                        xr = a[j]*a[j];
 
150
                        a[j+1] -= yi;
 
151
                        xr += (a[j+1]*a[j+1]);
 
152
                        pFFTAveBuf[k] = (1.0-1.0/m_AveSize)*pFFTAveBuf[k] +
 
153
                                                                        (1.0/m_AveSize)*(log10(xi+K_C) + K_B);
 
154
                        pFFTAveBuf[j] = (1.0-1.0/m_AveSize)*pFFTAveBuf[j] +
 
155
                                                                        (1.0/m_AveSize)*(log10(xr+K_C) + K_B);
 
156
                }
 
157
        }
 
158
        else
 
159
        {
 
160
 
 
161
                for (k = 2; k <= m; k += 2 ) 
 
162
                {
 
163
                        j = n - k;
 
164
                        kk += ks;
 
165
                        wkr = 0.5 - c[nc - kk];
 
166
                        wki = c[kk];
 
167
                        xr = a[k] - a[j];
 
168
                        xi = a[k + 1] + a[j + 1];
 
169
                        yr = wkr * xr - wki * xi;
 
170
                        yi = wkr * xi + wki * xr;
 
171
                        a[k] -= yr;
 
172
                        xi = a[k]*a[k];
 
173
                        a[k+1] -= yi;
 
174
                        xi += ( a[k+1]*a[k+1]);
 
175
                        a[j] += yr;
 
176
                        xr = a[j]*a[j];
 
177
                        a[j+1] -= yi;
 
178
                        xr += (a[j+1]*a[j+1]);
 
179
                        pFFTAveBuf[k] = log10(xi+K_C)+K_B;
 
180
                        pFFTAveBuf[j] = log10(xr+K_C)+K_B;
 
181
                }
 
182
        }
 
183
}
 
184
 
 
185
/* -------- initializing routines -------- */
 
186
void Cfft::makewt(int nw, int *ip, double *w)
 
187
{
 
188
    int nwh, j;
 
189
    double delta, x, y;
 
190
    
 
191
    ip[0] = nw;
 
192
    ip[1] = 1;
 
193
    if (nw > 2) {
 
194
        nwh = nw >> 1;
 
195
        delta = atan(1.0) / nwh;
 
196
        w[0] = 1;
 
197
        w[1] = 0;
 
198
        w[nwh] = cos(delta * nwh);
 
199
        w[nwh + 1] = w[nwh];
 
200
        for (j = 2; j < nwh; j += 2) {
 
201
            x = cos(delta * j);
 
202
            y = sin(delta * j);
 
203
            w[j] = x;
 
204
            w[j + 1] = y;
 
205
            w[nw - j] = y;
 
206
 
 
207
            w[nw - j + 1] = x;
 
208
        }
 
209
        bitrv2(nw, ip + 2, w);
 
210
    }
 
211
}
 
212
 
 
213
 
 
214
void Cfft::makect(int nc, int *ip, double *c)
 
215
{
 
216
    int nch, j;
 
217
    double delta;
 
218
    
 
219
    ip[1] = nc;
 
220
    if (nc > 1) {
 
221
        nch = nc >> 1;
 
222
        delta = atan(1.0) / nch;
 
223
        c[0] = cos(delta * nch);
 
224
        c[nch] = 0.5 * c[0];
 
225
        for (j = 1; j < nch; j++) {
 
226
            c[j] = 0.5 * cos(delta * j);
 
227
            c[nc - j] = 0.5 * sin(delta * j);
 
228
        }
 
229
    }
 
230
}
 
231
 
 
232
 
 
233
/* -------- child routines -------- */
 
234
void Cfft::bitrv2(int n, int *ip, double *a)
 
235
{
 
236
    int j, j1, k, k1, l, m, m2;
 
237
    double xr, xi;
 
238
    
 
239
    ip[0] = 0;
 
240
    l = n;
 
241
    m = 1;
 
242
    while ((m << 2) < l) {
 
243
        l >>= 1;
 
244
        for (j = 0; j < m; j++) {
 
245
            ip[m + j] = ip[j] + l;
 
246
        }
 
247
        m <<= 1;
 
248
    }
 
249
    if ((m << 2) > l) {
 
250
        for (k = 1; k < m; k++) {
 
251
            for (j = 0; j < k; j++) {
 
252
                j1 = (j << 1) + ip[k];
 
253
                k1 = (k << 1) + ip[j];
 
254
                xr = a[j1];
 
255
 
 
256
                xi = a[j1 + 1];
 
257
                a[j1] = a[k1];
 
258
                a[j1 + 1] = a[k1 + 1];
 
259
                a[k1] = xr;
 
260
                a[k1 + 1] = xi;
 
261
            }
 
262
        }
 
263
    } else {
 
264
        m2 = m << 1;
 
265
        for (k = 1; k < m; k++) {
 
266
            for (j = 0; j < k; j++) {
 
267
                j1 = (j << 1) + ip[k];
 
268
                k1 = (k << 1) + ip[j];
 
269
                xr = a[j1];
 
270
                xi = a[j1 + 1];
 
271
                a[j1] = a[k1];
 
272
                a[j1 + 1] = a[k1 + 1];
 
273
                a[k1] = xr;
 
274
                a[k1 + 1] = xi;
 
275
                j1 += m2;
 
276
                k1 += m2;
 
277
                xr = a[j1];
 
278
                xi = a[j1 + 1];
 
279
                a[j1] = a[k1];
 
280
                a[j1 + 1] = a[k1 + 1];
 
281
                a[k1] = xr;
 
282
                a[k1 + 1] = xi;
 
283
            }
 
284
        }
 
285
    }
 
286
}
 
287
 
 
288
void Cfft::cftfsub(int n, double *a, double *w)
 
289
{
 
290
    int j, j1, j2, j3, l;
 
291
    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
 
292
    
 
293
    l = 2;
 
294
    if (n > 8) {
 
295
        cft1st(n, a, w);
 
296
        l = 8;
 
297
        while ((l << 2) < n) {
 
298
            cftmdl(n, l, a, w);
 
299
            l <<= 2;
 
300
        }
 
301
    }
 
302
    if ((l << 2) == n) {
 
303
        for (j = 0; j < l; j += 2) {
 
304
            j1 = j + l;
 
305
            j2 = j1 + l;
 
306
            j3 = j2 + l;
 
307
            x0r = a[j] + a[j1];
 
308
            x0i = a[j + 1] + a[j1 + 1];
 
309
            x1r = a[j] - a[j1];
 
310
            x1i = a[j + 1] - a[j1 + 1];
 
311
            x2r = a[j2] + a[j3];
 
312
            x2i = a[j2 + 1] + a[j3 + 1];
 
313
            x3r = a[j2] - a[j3];
 
314
            x3i = a[j2 + 1] - a[j3 + 1];
 
315
            a[j] = x0r + x2r;
 
316
            a[j + 1] = x0i + x2i;
 
317
            a[j2] = x0r - x2r;
 
318
            a[j2 + 1] = x0i - x2i;
 
319
            a[j1] = x1r - x3i;
 
320
            a[j1 + 1] = x1i + x3r;
 
321
            a[j3] = x1r + x3i;
 
322
            a[j3 + 1] = x1i - x3r;
 
323
        }
 
324
    } else {
 
325
        for (j = 0; j < l; j += 2) {
 
326
            j1 = j + l;
 
327
            x0r = a[j] - a[j1];
 
328
            x0i = a[j + 1] - a[j1 + 1];
 
329
            a[j] += a[j1];
 
330
            a[j + 1] += a[j1 + 1];
 
331
            a[j1] = x0r;
 
332
            a[j1 + 1] = x0i;
 
333
        }
 
334
    }
 
335
}
 
336
 
 
337
 
 
338
 
 
339
void Cfft::cft1st(int n, double *a, double *w)
 
340
{
 
341
    int j, k1, k2;
 
342
    double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
 
343
    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
 
344
    
 
345
    x0r = a[0] + a[2];
 
346
    x0i = a[1] + a[3];
 
347
    x1r = a[0] - a[2];
 
348
    x1i = a[1] - a[3];
 
349
    x2r = a[4] + a[6];
 
350
    x2i = a[5] + a[7];
 
351
    x3r = a[4] - a[6];
 
352
    x3i = a[5] - a[7];
 
353
    a[0] = x0r + x2r;
 
354
    a[1] = x0i + x2i;
 
355
    a[4] = x0r - x2r;
 
356
    a[5] = x0i - x2i;
 
357
    a[2] = x1r - x3i;
 
358
    a[3] = x1i + x3r;
 
359
    a[6] = x1r + x3i;
 
360
    a[7] = x1i - x3r;
 
361
    wk1r = w[2];
 
362
    x0r = a[8] + a[10];
 
363
    x0i = a[9] + a[11];
 
364
    x1r = a[8] - a[10];
 
365
    x1i = a[9] - a[11];
 
366
    x2r = a[12] + a[14];
 
367
    x2i = a[13] + a[15];
 
368
    x3r = a[12] - a[14];
 
369
    x3i = a[13] - a[15];
 
370
    a[8] = x0r + x2r;
 
371
    a[9] = x0i + x2i;
 
372
    a[12] = x2i - x0i;
 
373
    a[13] = x0r - x2r;
 
374
    x0r = x1r - x3i;
 
375
    x0i = x1i + x3r;
 
376
    a[10] = wk1r * (x0r - x0i);
 
377
    a[11] = wk1r * (x0r + x0i);
 
378
    x0r = x3i + x1r;
 
379
    x0i = x3r - x1i;
 
380
    a[14] = wk1r * (x0i - x0r);
 
381
    a[15] = wk1r * (x0i + x0r);
 
382
    k1 = 0;
 
383
    for (j = 16; j < n; j += 16) {
 
384
        k1 += 2;
 
385
        k2 = k1 << 1;
 
386
        wk2r = w[k1];
 
387
        wk2i = w[k1 + 1];
 
388
        wk1r = w[k2];
 
389
        wk1i = w[k2 + 1];
 
390
        wk3r = wk1r - 2 * wk2i * wk1i;
 
391
        wk3i = 2 * wk2i * wk1r - wk1i;
 
392
        x0r = a[j] + a[j + 2];
 
393
        x0i = a[j + 1] + a[j + 3];
 
394
        x1r = a[j] - a[j + 2];
 
395
        x1i = a[j + 1] - a[j + 3];
 
396
        x2r = a[j + 4] + a[j + 6];
 
397
        x2i = a[j + 5] + a[j + 7];
 
398
        x3r = a[j + 4] - a[j + 6];
 
399
        x3i = a[j + 5] - a[j + 7];
 
400
        a[j] = x0r + x2r;
 
401
        a[j + 1] = x0i + x2i;
 
402
        x0r -= x2r;
 
403
        x0i -= x2i;
 
404
        a[j + 4] = wk2r * x0r - wk2i * x0i;
 
405
        a[j + 5] = wk2r * x0i + wk2i * x0r;
 
406
        x0r = x1r - x3i;
 
407
        x0i = x1i + x3r;
 
408
        a[j + 2] = wk1r * x0r - wk1i * x0i;
 
409
        a[j + 3] = wk1r * x0i + wk1i * x0r;
 
410
        x0r = x1r + x3i;
 
411
        x0i = x1i - x3r;
 
412
        a[j + 6] = wk3r * x0r - wk3i * x0i;
 
413
        a[j + 7] = wk3r * x0i + wk3i * x0r;
 
414
        wk1r = w[k2 + 2];
 
415
        wk1i = w[k2 + 3];
 
416
        wk3r = wk1r - 2 * wk2r * wk1i;
 
417
        wk3i = 2 * wk2r * wk1r - wk1i;
 
418
        x0r = a[j + 8] + a[j + 10];
 
419
        x0i = a[j + 9] + a[j + 11];
 
420
        x1r = a[j + 8] - a[j + 10];
 
421
        x1i = a[j + 9] - a[j + 11];
 
422
        x2r = a[j + 12] + a[j + 14];
 
423
        x2i = a[j + 13] + a[j + 15];
 
424
        x3r = a[j + 12] - a[j + 14];
 
425
        x3i = a[j + 13] - a[j + 15];
 
426
        a[j + 8] = x0r + x2r;
 
427
        a[j + 9] = x0i + x2i;
 
428
        x0r -= x2r;
 
429
        x0i -= x2i;
 
430
        a[j + 12] = -wk2i * x0r - wk2r * x0i;
 
431
        a[j + 13] = -wk2i * x0i + wk2r * x0r;
 
432
        x0r = x1r - x3i;
 
433
        x0i = x1i + x3r;
 
434
        a[j + 10] = wk1r * x0r - wk1i * x0i;
 
435
        a[j + 11] = wk1r * x0i + wk1i * x0r;
 
436
        x0r = x1r + x3i;
 
437
        x0i = x1i - x3r;
 
438
        a[j + 14] = wk3r * x0r - wk3i * x0i;
 
439
        a[j + 15] = wk3r * x0i + wk3i * x0r;
 
440
    }
 
441
}
 
442
 
 
443
 
 
444
void Cfft::cftmdl(int n, int l, double *a, double *w)
 
445
{
 
446
    int j, j1, j2, j3, k, k1, k2, m, m2;
 
447
    double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
 
448
    double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
 
449
    
 
450
    m = l << 2;
 
451
    for (j = 0; j < l; j += 2) {
 
452
        j1 = j + l;
 
453
        j2 = j1 + l;
 
454
        j3 = j2 + l;
 
455
        x0r = a[j] + a[j1];
 
456
        x0i = a[j + 1] + a[j1 + 1];
 
457
        x1r = a[j] - a[j1];
 
458
        x1i = a[j + 1] - a[j1 + 1];
 
459
        x2r = a[j2] + a[j3];
 
460
        x2i = a[j2 + 1] + a[j3 + 1];
 
461
        x3r = a[j2] - a[j3];
 
462
        x3i = a[j2 + 1] - a[j3 + 1];
 
463
        a[j] = x0r + x2r;
 
464
        a[j + 1] = x0i + x2i;
 
465
        a[j2] = x0r - x2r;
 
466
        a[j2 + 1] = x0i - x2i;
 
467
        a[j1] = x1r - x3i;
 
468
        a[j1 + 1] = x1i + x3r;
 
469
        a[j3] = x1r + x3i;
 
470
        a[j3 + 1] = x1i - x3r;
 
471
    }
 
472
    wk1r = w[2];
 
473
    for (j = m; j < l + m; j += 2) {
 
474
        j1 = j + l;
 
475
        j2 = j1 + l;
 
476
        j3 = j2 + l;
 
477
        x0r = a[j] + a[j1];
 
478
        x0i = a[j + 1] + a[j1 + 1];
 
479
        x1r = a[j] - a[j1];
 
480
        x1i = a[j + 1] - a[j1 + 1];
 
481
        x2r = a[j2] + a[j3];
 
482
        x2i = a[j2 + 1] + a[j3 + 1];
 
483
        x3r = a[j2] - a[j3];
 
484
        x3i = a[j2 + 1] - a[j3 + 1];
 
485
        a[j] = x0r + x2r;
 
486
        a[j + 1] = x0i + x2i;
 
487
        a[j2] = x2i - x0i;
 
488
        a[j2 + 1] = x0r - x2r;
 
489
        x0r = x1r - x3i;
 
490
        x0i = x1i + x3r;
 
491
        a[j1] = wk1r * (x0r - x0i);
 
492
        a[j1 + 1] = wk1r * (x0r + x0i);
 
493
        x0r = x3i + x1r;
 
494
        x0i = x3r - x1i;
 
495
        a[j3] = wk1r * (x0i - x0r);
 
496
        a[j3 + 1] = wk1r * (x0i + x0r);
 
497
    }
 
498
    k1 = 0;
 
499
    m2 = m << 1;
 
500
    for (k = m2; k < n; k += m2) {
 
501
        k1 += 2;
 
502
        k2 = k1 << 1;
 
503
        wk2r = w[k1];
 
504
        wk2i = w[k1 + 1];
 
505
        wk1r = w[k2];
 
506
        wk1i = w[k2 + 1];
 
507
        wk3r = wk1r - 2 * wk2i * wk1i;
 
508
        wk3i = 2 * wk2i * wk1r - wk1i;
 
509
        for (j = k; j < l + k; j += 2) {
 
510
            j1 = j + l;
 
511
            j2 = j1 + l;
 
512
            j3 = j2 + l;
 
513
            x0r = a[j] + a[j1];
 
514
            x0i = a[j + 1] + a[j1 + 1];
 
515
            x1r = a[j] - a[j1];
 
516
            x1i = a[j + 1] - a[j1 + 1];
 
517
            x2r = a[j2] + a[j3];
 
518
            x2i = a[j2 + 1] + a[j3 + 1];
 
519
            x3r = a[j2] - a[j3];
 
520
            x3i = a[j2 + 1] - a[j3 + 1];
 
521
            a[j] = x0r + x2r;
 
522
            a[j + 1] = x0i + x2i;
 
523
            x0r -= x2r;
 
524
            x0i -= x2i;
 
525
 
 
526
            a[j2] = wk2r * x0r - wk2i * x0i;
 
527
            a[j2 + 1] = wk2r * x0i + wk2i * x0r;
 
528
            x0r = x1r - x3i;
 
529
            x0i = x1i + x3r;
 
530
            a[j1] = wk1r * x0r - wk1i * x0i;
 
531
            a[j1 + 1] = wk1r * x0i + wk1i * x0r;
 
532
            x0r = x1r + x3i;
 
533
            x0i = x1i - x3r;
 
534
            a[j3] = wk3r * x0r - wk3i * x0i;
 
535
            a[j3 + 1] = wk3r * x0i + wk3i * x0r;
 
536
        }
 
537
        wk1r = w[k2 + 2];
 
538
        wk1i = w[k2 + 3];
 
539
        wk3r = wk1r - 2 * wk2r * wk1i;
 
540
        wk3i = 2 * wk2r * wk1r - wk1i;
 
541
        for (j = k + m; j < l + (k + m); j += 2) {
 
542
            j1 = j + l;
 
543
            j2 = j1 + l;
 
544
            j3 = j2 + l;
 
545
            x0r = a[j] + a[j1];
 
546
            x0i = a[j + 1] + a[j1 + 1];
 
547
            x1r = a[j] - a[j1];
 
548
            x1i = a[j + 1] - a[j1 + 1];
 
549
            x2r = a[j2] + a[j3];
 
550
            x2i = a[j2 + 1] + a[j3 + 1];
 
551
            x3r = a[j2] - a[j3];
 
552
            x3i = a[j2 + 1] - a[j3 + 1];
 
553
            a[j] = x0r + x2r;
 
554
            a[j + 1] = x0i + x2i;
 
555
            x0r -= x2r;
 
556
            x0i -= x2i;
 
557
            a[j2] = -wk2i * x0r - wk2r * x0i;
 
558
            a[j2 + 1] = -wk2i * x0i + wk2r * x0r;
 
559
            x0r = x1r - x3i;
 
560
            x0i = x1i + x3r;
 
561
            a[j1] = wk1r * x0r - wk1i * x0i;
 
562
            a[j1 + 1] = wk1r * x0i + wk1i * x0r;
 
563
            x0r = x1r + x3i;
 
564
            x0i = x1i - x3r;
 
565
            a[j3] = wk3r * x0r - wk3i * x0i;
 
566
            a[j3 + 1] = wk3r * x0i + wk3i * x0r;
 
567
        }
 
568
    }
 
569
}
 
570
 
 
571
 
 
572
 
 
573
 
 
574
 
 
575
 
 
576