?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
dlzma: added fix for very rare out-of-bounds write in LZMA encoder (it can happen...
[iv.d.git] / trestex / mbandeq_j.d
blob04fcceadceb8f5181c3c724964895eee902f9a90
1 // Shibatch Super Equalizer ver 0.03 for winamp
2 // written by Naoki Shibata shibatch@users.sourceforge.net
3 // LGPL
4 module mbandeq_j /*is aliced*/;
5 private:
6 import iv.alice;
7
8 //import fftsg;
9
10 //version = mbeq_debug_output;
11
12 version(mbeq_debug_output) import iv.cmdcon;
13
14
15 // ////////////////////////////////////////////////////////////////////////// //
16 // Real Discrete Fourier Transform wrapper
17 void rfft (int n, int isign, REAL* x) {
18 import std.math : sqrt;
19 __gshared int ipsize = 0, wsize=0;
20 __gshared int* ip = null;
21 __gshared REAL* w = null;
22 int newipsize, newwsize;
23
24 if (n == 0) {
25 import core.stdc.stdlib : free;
26 free(ip); ip = null; ipsize = 0;
27 free(w); w = null; wsize = 0;
28 return;
29 }
30
31 newipsize = 2+cast(int)sqrt(cast(REAL)(n/2));
32 if (newipsize > ipsize) {
33 import core.stdc.stdlib : realloc;
34 ipsize = newipsize;
35 ip = cast(int*)realloc(ip, int.sizeof*ipsize);
36 ip[0] = 0;
37 }
38
39 newwsize = n/2;
40 if (newwsize > wsize) {
41 import core.stdc.stdlib : realloc;
42 wsize = newwsize;
43 w = cast(REAL*)realloc(w, REAL.sizeof*wsize);
44 }
45
46 rdft(n, isign, x, ip, w);
47 }
48
49
50 // ////////////////////////////////////////////////////////////////////////// //
51 enum M = 15;
52
53 //#define RINT(x) ((x) >= 0 ? ((int)((x)+0.5)) : ((int)((x)-0.5)))
54 int RINT (REAL x) pure nothrow @safe @nogc { pragma(inline, true); return (x >= 0 ? cast(int)(x+0.5f) : cast(int)(x-0.5f)); }
55
56 enum DITHERLEN = 65536;
57
58 __gshared REAL[M+1] fact;
59 __gshared REAL aa = 96;
60 __gshared REAL iza;
61 __gshared REAL* lires, lires1, lires2;
62 __gshared REAL* rires, rires1, rires2;
63 __gshared REAL* irest;
64 __gshared REAL* fsamples;
65 __gshared int chg_ires, cur_ires;
66 __gshared int winlen, winlenbit, tabsize, nbufsamples;
67 __gshared short* inbuf;
68 __gshared REAL* outbuf;
69 enum dither = false;
70 static if (dither) {
71 __gshared REAL* ditherbuf;
72 __gshared uint ditherptr = 0;
73 __gshared REAL smphm1 = 0, smphm2 = 0;
74 }
75
76 enum NCH = 2;
77
78 enum NBANDS = 17;
79 static immutable REAL[NBANDS] bands = [
80 65.406392, 92.498606, 130.81278, 184.99721, 261.62557, 369.99442, 523.25113,
81 739.9884 , 1046.5023, 1479.9768, 2093.0045, 2959.9536, 4186.0091, 5919.9072,
82 8372.0181, 11839.814, 16744.036,
83 ];
84
85
86 REAL alpha (REAL a) nothrow @trusted @nogc {
87 import std.math : pow;
88 if (a <= 21) return 0;
89 if (a <= 50) return 0.5842*pow(a-21, 0.4)+0.07886*(a-21);
90 return 0.1102*(a-8.7);
91 }
92
93
94 REAL izero (REAL x) nothrow @trusted @nogc {
95 import std.math : pow;
96 REAL ret = 1;
97 foreach (immutable m; 1..M+1) {
98 REAL t = pow(x/2, m)/fact[m];
99 ret += t*t;
100 }
101 return ret;
102 }
103
104
105 /// wb: window length, bits
106 public void mbeqInit (int wb=14) {
107 import core.stdc.stdlib : malloc, calloc, free;
108
109 if (lires1 !is null) free(lires1);
110 if (lires2 !is null) free(lires2);
111 if (rires1 !is null) free(rires1);
112 if (rires2 !is null) free(rires2);
113 if (irest !is null) free(irest);
114 if (fsamples !is null) free(fsamples);
115 if (inbuf !is null) free(inbuf);
116 if (outbuf !is null) free(outbuf);
117 static if (dither) { if (ditherbuf !is null) free(ditherbuf); }
118
119 winlen = (1<<(wb-1))-1;
120 winlenbit = wb;
121 tabsize = 1<<wb;
122
123 //{ import core.stdc.stdio; printf("winlen=%u\n", winlen); }
124
125 lires1 = cast(REAL*)malloc(REAL.sizeof*tabsize);
126 lires2 = cast(REAL*)malloc(REAL.sizeof*tabsize);
127 rires1 = cast(REAL*)malloc(REAL.sizeof*tabsize);
128 rires2 = cast(REAL*)malloc(REAL.sizeof*tabsize);
129 irest = cast(REAL*)malloc(REAL.sizeof*tabsize);
130 fsamples = cast(REAL*)malloc(REAL.sizeof*tabsize);
131 inbuf = cast(short*)calloc(winlen*NCH, int.sizeof);
132 outbuf = cast(REAL*)calloc(tabsize*NCH, REAL.sizeof);
133 static if (dither) { ditherbuf = cast(REAL*)malloc(REAL.sizeof*DITHERLEN); }
134
135 lires = lires1;
136 rires = rires1;
137 cur_ires = 1;
138 chg_ires = 1;
139
140 static if (dither) {
141 foreach (immutable i; 0..DITHERLEN) {
142 import std.random;
143 //ditherbuf[i] = (REAL(rand())/RAND_MAX-0.5);
144 ditherbuf[i] = uniform!"[]"(0, 32760)/32760.0f-0.5f;
145 }
146 }
147
148 foreach (immutable i; 0..M+1) {
149 fact[i] = 1;
150 foreach (immutable j; 1..i+1) fact[i] *= j;
151 }
152
153 iza = izero(alpha(aa));
154 }
155
156
157 // -(N-1)/2 <= n <= (N-1)/2
158 REAL win (REAL n, int N) nothrow @trusted @nogc {
159 pragma(inline, true);
160 import std.math : sqrt;
161 return izero(alpha(aa)*sqrt(1-4*n*n/((N-1)*(N-1))))/iza;
162 }
163
164
165 REAL sinc (REAL x) nothrow @trusted @nogc {
166 pragma(inline, true);
167 import std.math : sin;
168 return (x == 0 ? 1 : sin(x)/x);
169 }
170
171
172 REAL hn_lpf (int n, REAL f, REAL fs) nothrow @trusted @nogc {
173 pragma(inline, true);
174 import std.math : PI;
175 immutable REAL t = 1/fs;
176 immutable REAL omega = 2*PI*f;
177 return 2*f*t*sinc(n*omega*t);
178 }
179
180
181 REAL hn_imp (int n) nothrow @trusted @nogc {
182 pragma(inline, true);
183 return (n == 0 ? 1.0 : 0.0);
184 }
185
186
187 REAL getLower (int i) nothrow @trusted @nogc {
188 pragma(inline, true);
189 if (i < 0) assert(0, "wtf?");
190 if (i > NBANDS+1) assert(0, "wtf?");
191 return (i == 0 ? 0 : bands[i-1]);
192 }
193
194
195 REAL getUpper (int i) nothrow @trusted @nogc {
196 pragma(inline, true);
197 if (i < 0) assert(0, "wtf?");
198 if (i > NBANDS+1) assert(0, "wtf?");
199 return (i == NBANDS ? 48000/2 : bands[i]);
200 }
201
202
203 REAL hn (int n, const(REAL)[] gains, REAL fs) nothrow @trusted @nogc {
204 REAL lhn = hn_lpf(n, bands[0], fs);
205 REAL ret = gains[0]*lhn;
206 foreach (immutable i, REAL gv; gains) {
207 if (i == 0) continue;
208 assert(i <= NBANDS);
209 //immutable REAL lower = getLower(i);
210 immutable REAL upper = getUpper(i);
211 if (upper >= fs/2) {
212 // the last one
213 ret += gv*(hn_imp(n)-lhn);
214 return ret;
215 }
216 REAL lhn2 = hn_lpf(n, upper, fs);
217 ret += gv*(lhn2-lhn);
218 lhn = lhn2;
219 }
220 assert(0, "wtf?!");
221 }
222
223
224 void processParam (REAL[] gains, const(REAL)[] bc, const(REAL)[] gainmods=null) {
225 import std.math : pow;
226 foreach (immutable i; 0..NBANDS+1) {
227 gains[i] = bc[i];
228 immutable REAL gm = (i < gainmods.length ? gainmods[i] : 1);
229 gains[i] *= pow(10, gm/20);
230 //if (i < gainmods.length) gains[i] *= pow(10, gainmods[i]/20);
231 }
232 }
233
234
235 ///
236 public void mbeqMakeTable (const(REAL)[] lbc, const(REAL)[] rbc, REAL fs, const(REAL)[] gainmods=null) {
237 int i, cires = cur_ires;
238 REAL* nires;
239 REAL[NBANDS+1] gains;
240
241 if (fs <= 0) return;
242
243 // left
244 processParam(gains[], lbc, gainmods);
245 version(mbeq_debug_output) foreach (immutable gi, immutable REAL gv; gains[]) conwriteln("L: ", getLower(gi, fs), "Hz to ", getUpper(gi, fs), "Hz, ", gv);
246 for (i = 0; i < winlen; ++i) irest[i] = hn(i-winlen/2, gains[], fs)*win(i-winlen/2, winlen);
247 for (; i < tabsize; ++i) irest[i] = 0;
248 rfft(tabsize, 1, irest);
249 nires = (cires == 1 ? lires2 : lires1);
250 for (i = 0; i < tabsize; ++i) nires[i] = irest[i];
251
252 // right
253 processParam(gains[], rbc, gainmods);
254 version(mbeq_debug_output) foreach (immutable gi, immutable REAL gv; gains[]) conwriteln("R: ", getLower(gi, fs), "Hz to ", getUpper(gi, fs), "Hz, ", gv);
255 for (i = 0; i < winlen; ++i) irest[i] = hn(i-winlen/2, gains[], fs)*win(i-winlen/2, winlen);
256 for (; i < tabsize; ++i) irest[i] = 0;
257 rfft(tabsize, 1, irest);
258 nires = (cires == 1 ? rires2 : rires1);
259 for (i = 0; i < tabsize; ++i) nires[i] = irest[i];
260
261 // done
262 chg_ires = (cires == 1 ? 2 : 1);
263 }
264
265
266 ///
267 public void mbeqQuit () {
268 import core.stdc.stdlib : free;
269
270 if (lires1 !is null) free(lires1);
271 if (lires2 !is null) free(lires2);
272 if (rires1 !is null) free(rires1);
273 if (rires2 !is null) free(rires2);
274 if (irest !is null) free(irest);
275 if (fsamples !is null) free(fsamples);
276 if (inbuf !is null) free(inbuf);
277 if (outbuf !is null) free(outbuf);
278
279 lires1 = null;
280 lires2 = null;
281 rires1 = null;
282 rires2 = null;
283 irest = null;
284 fsamples = null;
285 inbuf = null;
286 outbuf = null;
287
288 rfft(0, 0, null);
289 }
290
291
292 ///
293 public void mbeqClearbuf (int bps, int srate) {
294 nbufsamples = 0;
295 //foreach (immutable i; 0..tabsize*NCH) outbuf[i] = 0;
296 outbuf[0..tabsize*NCH] = 0;
297 }
298
299
300 ///
301 public int mbeqModifySamples (void* buf, int nsamples, int nch, int bps) {
302 int i, p, ch;
303 REAL* ires;
304 int amax = (1<<(bps-1))-1;
305 int amin = -(1<<(bps-1));
306 //static REAL smphm1 = 0, smphm2 = 0;
307
308 if (chg_ires) {
309 cur_ires = chg_ires;
310 lires = (cur_ires == 1 ? lires1 : lires2);
311 rires = (cur_ires == 1 ? rires1 : rires2);
312 chg_ires = 0;
313 }
314
315 p = 0;
316
317 while (nbufsamples+nsamples >= winlen) {
318 //version(mbeq_debug_output) conwriteln("nbufsamples+nsamples=", nbufsamples+nsamples, "; winlen=", winlen);
319 switch (bps) {
320 case 8:
321 for (i = 0; i < (winlen-nbufsamples)*nch; ++i) {
322 inbuf[nbufsamples*nch+i] = (cast(ubyte*)buf)[i+p*nch]-0x80;
323 REAL s = outbuf[nbufsamples*nch+i];
324 static if (dither) {
325 s -= smphm1;
326 REAL u = s;
327 s += ditherbuf[(ditherptr++)&(DITHERLEN-1)];
328 if (s < amin) s = amin;
329 if (amax < s) s = amax;
330 s = RINT(s);
331 smphm1 = s-u;
332 (cast(ubyte*)buf)[i+p*nch] = cast(ubyte)(s+0x80);
333 } else {
334 if (s < amin) s = amin;
335 if (amax < s) s = amax;
336 (cast(ubyte*)buf)[i+p*nch] = cast(ubyte)(RINT(s)+0x80);
337 }
338 }
339 for( i = winlen*nch; i < tabsize*nch; ++i) outbuf[i-winlen*nch] = outbuf[i];
340 break;
341 case 16:
342 for (i = 0; i < (winlen-nbufsamples)*nch; ++i) {
343 inbuf[nbufsamples*nch+i] = (cast(short*)buf)[i+p*nch];
344 REAL s = outbuf[nbufsamples*nch+i];
345 static if (dither) {
346 s -= smphm1;
347 REAL u = s;
348 s += ditherbuf[(ditherptr++)&(DITHERLEN-1)];
349 if (s < amin) s = amin;
350 if (amax < s) s = amax;
351 s = RINT(s);
352 smphm1 = s-u;
353 (cast(short*)buf)[i+p*nch] = cast(short)s;
354 } else {
355 if (s < amin) s = amin;
356 if (amax < s) s = amax;
357 (cast(short*)buf)[i+p*nch] = cast(short)RINT(s);
358 }
359 }
360 for (i = winlen*nch; i < tabsize*nch; ++i) outbuf[i-winlen*nch] = outbuf[i];
361 break;
362
363 case 24:
364 for (i = 0; i < (winlen-nbufsamples)*nch; ++i) {
365 (cast(int*)inbuf)[nbufsamples*nch+i] =
366 ((cast(ubyte*)buf)[(i+p*nch)*3])|
367 ((cast(ubyte*)buf)[(i+p*nch)*3+1]<<8)|
368 ((cast(byte*)buf)[(i+p*nch)*3+2]<<16);
369 REAL s = outbuf[nbufsamples*nch+i];
370 //static if (dither) s += ditherbuf[(ditherptr++)&(DITHERLEN-1)];
371 if (s < amin) s = amin;
372 if (amax < s) s = amax;
373 int s2 = RINT(s);
374 (cast(ubyte*)buf)[(i+p*nch)*3+0] = s2&255; s2 >>= 8;
375 (cast(ubyte*)buf)[(i+p*nch)*3+1] = s2&255; s2 >>= 8;
376 (cast(ubyte*)buf)[(i+p*nch)*3+2] = s2&255;
377 }
378 for (i = winlen*nch; i < tabsize*nch; ++i) outbuf[i-winlen*nch] = outbuf[i];
379 break;
380 default: assert(0);
381 }
382
383 p += winlen-nbufsamples;
384 //{ import core.stdc.stdio; printf("old nsamples: %d\n", nsamples); }
385 nsamples -= winlen-nbufsamples;
386 //{ import core.stdc.stdio; printf("new nsamples: %d\n", nsamples); }
387 nbufsamples = 0;
388
389 for (ch = 0; ch < nch; ++ch) {
390 ires = (ch == 0 ? lires : rires);
391
392 if (bps == 24) {
393 for (i = 0; i < winlen; ++i) fsamples[i] = (cast(int*)inbuf)[nch*i+ch];
394 } else {
395 for (i = 0; i < winlen; ++i) fsamples[i] = inbuf[nch*i+ch];
396 }
397
398 //for (i = winlen; i < tabsize; ++i) fsamples[i] = 0;
399 fsamples[winlen..tabsize] = 0;
400
401 rfft(tabsize, 1, fsamples);
402 fsamples[0] = ires[0]*fsamples[0];
403 fsamples[1] = ires[1]*fsamples[1];
404 for (i = 1; i < tabsize/2; ++i) {
405 REAL re = ires[i*2 ]*fsamples[i*2]-ires[i*2+1]*fsamples[i*2+1];
406 REAL im = ires[i*2+1]*fsamples[i*2]+ires[i*2 ]*fsamples[i*2+1];
407 fsamples[i*2 ] = re;
408 fsamples[i*2+1] = im;
409 }
410 rfft(tabsize, -1, fsamples);
411 /*if disabled:
412 {
413 for (i = winlen-1+winlen/2; i >= winlen/2; --i) fsamples[i] = fsamples[i-winlen/2]*tabsize/2;
414 for (; i >= 0; --i) fsamples[i] = 0;
415 }
416 */
417
418 for (i = 0; i < winlen; ++i) outbuf[i*nch+ch] += fsamples[i]/tabsize*2;
419 for (i = winlen; i < tabsize; ++i) outbuf[i*nch+ch] = fsamples[i]/tabsize*2;
420 }
421 }
422
423 switch (bps) {
424 case 8:
425 for (i = 0; i < nsamples*nch; ++i) {
426 inbuf[nbufsamples*nch+i] = (cast(ubyte*)buf)[i+p*nch]-0x80;
427 REAL s = outbuf[nbufsamples*nch+i];
428 static if (dither) {
429 s -= smphm1;
430 REAL u = s;
431 s += ditherbuf[(ditherptr++)&(DITHERLEN-1)];
432 if (s < amin) s = amin;
433 if (amax < s) s = amax;
434 s = RINT(s);
435 smphm1 = s-u;
436 (cast(ubyte*)buf)[i+p*nch] = cast(ubyte)(s+0x80);
437 } else {
438 if (s < amin) s = amin;
439 if (amax < s) s = amax;
440 (cast(ubyte*)buf)[i+p*nch] = cast(ubyte)(RINT(s)+0x80);
441 }
442 }
443 break;
444 case 16:
445 for (i = 0; i < nsamples*nch; ++i) {
446 //{ import core.stdc.stdio; printf("i=%u; nsamples*nch=%u; nbufsamples*nch+i=%u; i+p*nch=%u\n", i, nsamples*nch, nbufsamples*nch+i, i+p*nch); }
447 inbuf[nbufsamples*nch+i] = (cast(short*)buf)[i+p*nch];
448 REAL s = outbuf[nbufsamples*nch+i];
449 static if (dither) {
450 s -= smphm1;
451 REAL u = s;
452 s += ditherbuf[(ditherptr++)&(DITHERLEN-1)];
453 if (s < amin) s = amin;
454 if (amax < s) s = amax;
455 s = RINT(s);
456 smphm1 = s-u;
457 (cast(short*)buf)[i+p*nch] = cast(short)s;
458 } else {
459 if (s < amin) s = amin;
460 if (amax < s) s = amax;
461 (cast(short*)buf)[i+p*nch] = cast(short)RINT(s);
462 }
463 }
464 break;
465 case 24:
466 for (i = 0; i < nsamples*nch; ++i) {
467 (cast(int*)inbuf)[nbufsamples*nch+i] =
468 ((cast(ubyte*)buf)[(i+p*nch)*3])|
469 ((cast(ubyte*)buf)[(i+p*nch)*3+1]<<8)|
470 ((cast(byte*)buf)[(i+p*nch)*3+2]<<16);
471 REAL s = outbuf[nbufsamples*nch+i];
472 //static if (dither) s += ditherbuf[(ditherptr++)&(DITHERLEN-1)];
473 if (s < amin) s = amin;
474 if (amax < s) s = amax;
475 int s2 = RINT(s);
476 (cast(ubyte*)buf)[(i+p*nch)*3+0] = s2&255; s2 >>= 8;
477 (cast(ubyte*)buf)[(i+p*nch)*3+1] = s2&255; s2 >>= 8;
478 (cast(ubyte*)buf)[(i+p*nch)*3+2] = s2&255;
479 }
480 break;
481 default: assert(0);
482 }
483
484 p += nsamples;
485 nbufsamples += nsamples;
486
487 return p;
488 }
489
490
491 // ////////////////////////////////////////////////////////////////////////// //
492 private static immutable int[NBANDS+1] mbeqBandFreqs = [55, 77, 110, 156, 220, 311, 440, 622, 880, 1244, 1760, 2489, 3520, 4978, 7040, 9956, 14080, 19912]; ///
493 public __gshared int[NBANDS+2] mbeqLSliders; /// [0..96]; 0: preamp; 1..18: bands
494 public __gshared int[NBANDS+2] mbeqRSliders; /// [0..96]; 0: preamp; 1..18: bands
495 public __gshared REAL mbeqSampleRate = 48000; ///
496
497
498 public void mbeqSetBandsFromSliders () {
499 import std.math : pow;
500
501 REAL[NBANDS+1] lbands = 1.0;
502 REAL[NBANDS+1] rbands = 1.0;
503
504 immutable REAL lpreamp = (mbeqLSliders[0] == 96 ? 0 : pow(10, mbeqLSliders[0]/-20.0f));
505 immutable REAL rpreamp = (mbeqRSliders[0] == 96 ? 0 : pow(10, mbeqRSliders[0]/-20.0f));
506 version(mbeq_debug_output) conwriteln("lpreamp=", lpreamp, "; rpreamp=", rpreamp);
507 foreach (immutable i; 0..NBANDS+1) {
508 lbands[i] = (mbeqLSliders[i+1] == 96 ? 0 : lpreamp*pow(10, mbeqLSliders[i+1]/-20.0f));
509 rbands[i] = (mbeqRSliders[i+1] == 96 ? 0 : rpreamp*pow(10, mbeqRSliders[i+1]/-20.0f));
510 }
511 //mbeq_makeTable(lbands.ptr, rbands.ptr, paramroot, last_srate);
512 mbeqMakeTable(lbands[], rbands[], mbeqSampleRate);
513 version(mbeq_debug_output) conwriteln("lbands=", lbands);
514 version(mbeq_debug_output) conwriteln("rbands=", rbands);
515 }
516
517
518 public int mbeqBandCount () { pragma(inline, true); return NBANDS+2; } ///
519
520 /// last band is equal to samling rate
521 public int mbeqBandFreq (int idx) { pragma(inline, true); return (idx > 1 && idx < NBANDS+1 ? mbeqBandFreqs[idx-1] : (idx == NBANDS+1 ? 48000 : 0)); }
522
523
524 // ////////////////////////////////////////////////////////////////////////// //
525 private:
526 /*
527 * Copyright(C) 1996-2001 Takuya OOURA
528 * email: ooura@mmm.t.u-tokyo.ac.jp
529 * download: http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html
530 * You may use, copy, modify this code for any purpose and
531 * without fee. You may distribute this ORIGINAL package.
532 */
533 //module fftsg /*is aliced*/;
534 private nothrow @trusted @nogc:
535
536 /*public*/ alias REAL = float;
537
538 /*
539 Fast Fourier/Cosine/Sine Transform
540 dimension :one
541 data length :power of 2
542 decimation :frequency
543 radix :split-radix
544 data :inplace
545 table :use
546 functions
547 cdft: Complex Discrete Fourier Transform
548 rdft: Real Discrete Fourier Transform
549 ddct: Discrete Cosine Transform
550 ddst: Discrete Sine Transform
551 dfct: Cosine Transform of RDFT (Real Symmetric DFT)
552 dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
553 function prototypes
554 void cdft(int, int, REAL *, int *, REAL *);
555 void rdft(int, int, REAL *, int *, REAL *);
556 void ddct(int, int, REAL *, int *, REAL *);
557 void ddst(int, int, REAL *, int *, REAL *);
558 void dfct(int, REAL *, REAL *, int *, REAL *);
559 void dfst(int, REAL *, REAL *, int *, REAL *);
560
561
562 -------- Complex DFT (Discrete Fourier Transform) --------
563 [definition]
564 <case1>
565 X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
566 <case2>
567 X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
568 (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
569 [usage]
570 <case1>
571 ip[0] = 0; // first time only
572 cdft(2*n, 1, a, ip, w);
573 <case2>
574 ip[0] = 0; // first time only
575 cdft(2*n, -1, a, ip, w);
576 [parameters]
577 2*n :data length (int)
578 n >= 1, n = power of 2
579 a[0...2*n-1] :input/output data (REAL *)
580 input data
581 a[2*j] = Re(x[j]),
582 a[2*j+1] = Im(x[j]), 0<=j<n
583 output data
584 a[2*k] = Re(X[k]),
585 a[2*k+1] = Im(X[k]), 0<=k<n
586 ip[0...*] :work area for bit reversal (int *)
587 length of ip >= 2+sqrt(n)
588 strictly,
589 length of ip >=
590 2+(1<<(int)(log(n+0.5)/log(2))/2).
591 ip[0],ip[1] are pointers of the cos/sin table.
592 w[0...n/2-1] :cos/sin table (REAL *)
593 w[],ip[] are initialized if ip[0] == 0.
594 [remark]
595 Inverse of
596 cdft(2*n, -1, a, ip, w);
597 is
598 cdft(2*n, 1, a, ip, w);
599 for (j = 0; j <= 2 * n - 1; j++) {
600 a[j] *= 1.0 / n;
601 }
602 .
603
604
605 -------- Real DFT / Inverse of Real DFT --------
606 [definition]
607 <case1> RDFT
608 R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
609 I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
610 <case2> IRDFT (excluding scale)
611 a[k] = (R[0] + R[n/2]*cos(pi*k))/2 +
612 sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) +
613 sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
614 [usage]
615 <case1>
616 ip[0] = 0; // first time only
617 rdft(n, 1, a, ip, w);
618 <case2>
619 ip[0] = 0; // first time only
620 rdft(n, -1, a, ip, w);
621 [parameters]
622 n :data length (int)
623 n >= 2, n = power of 2
624 a[0...n-1] :input/output data (REAL *)
625 <case1>
626 output data
627 a[2*k] = R[k], 0<=k<n/2
628 a[2*k+1] = I[k], 0<k<n/2
629 a[1] = R[n/2]
630 <case2>
631 input data
632 a[2*j] = R[j], 0<=j<n/2
633 a[2*j+1] = I[j], 0<j<n/2
634 a[1] = R[n/2]
635 ip[0...*] :work area for bit reversal (int *)
636 length of ip >= 2+sqrt(n/2)
637 strictly,
638 length of ip >=
639 2+(1<<(int)(log(n/2+0.5)/log(2))/2).
640 ip[0],ip[1] are pointers of the cos/sin table.
641 w[0...n/2-1] :cos/sin table (REAL *)
642 w[],ip[] are initialized if ip[0] == 0.
643 [remark]
644 Inverse of
645 rdft(n, 1, a, ip, w);
646 is
647 rdft(n, -1, a, ip, w);
648 for (j = 0; j <= n - 1; j++) {
649 a[j] *= 2.0 / n;
650 }
651 .
652
653
654 -------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
655 [definition]
656 <case1> IDCT (excluding scale)
657 C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
658 <case2> DCT
659 C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
660 [usage]
661 <case1>
662 ip[0] = 0; // first time only
663 ddct(n, 1, a, ip, w);
664 <case2>
665 ip[0] = 0; // first time only
666 ddct(n, -1, a, ip, w);
667 [parameters]
668 n :data length (int)
669 n >= 2, n = power of 2
670 a[0...n-1] :input/output data (REAL *)
671 output data
672 a[k] = C[k], 0<=k<n
673 ip[0...*] :work area for bit reversal (int *)
674 length of ip >= 2+sqrt(n/2)
675 strictly,
676 length of ip >=
677 2+(1<<(int)(log(n/2+0.5)/log(2))/2).
678 ip[0],ip[1] are pointers of the cos/sin table.
679 w[0...n*5/4-1] :cos/sin table (REAL *)
680 w[],ip[] are initialized if ip[0] == 0.
681 [remark]
682 Inverse of
683 ddct(n, -1, a, ip, w);
684 is
685 a[0] *= 0.5;
686 ddct(n, 1, a, ip, w);
687 for (j = 0; j <= n - 1; j++) {
688 a[j] *= 2.0 / n;
689 }
690 .
691
692
693 -------- DST (Discrete Sine Transform) / Inverse of DST --------
694 [definition]
695 <case1> IDST (excluding scale)
696 S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
697 <case2> DST
698 S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
699 [usage]
700 <case1>
701 ip[0] = 0; // first time only
702 ddst(n, 1, a, ip, w);
703 <case2>
704 ip[0] = 0; // first time only
705 ddst(n, -1, a, ip, w);
706 [parameters]
707 n :data length (int)
708 n >= 2, n = power of 2
709 a[0...n-1] :input/output data (REAL *)
710 <case1>
711 input data
712 a[j] = A[j], 0<j<n
713 a[0] = A[n]
714 output data
715 a[k] = S[k], 0<=k<n
716 <case2>
717 output data
718 a[k] = S[k], 0<k<n
719 a[0] = S[n]
720 ip[0...*] :work area for bit reversal (int *)
721 length of ip >= 2+sqrt(n/2)
722 strictly,
723 length of ip >=
724 2+(1<<(int)(log(n/2+0.5)/log(2))/2).
725 ip[0],ip[1] are pointers of the cos/sin table.
726 w[0...n*5/4-1] :cos/sin table (REAL *)
727 w[],ip[] are initialized if ip[0] == 0.
728 [remark]
729 Inverse of
730 ddst(n, -1, a, ip, w);
731 is
732 a[0] *= 0.5;
733 ddst(n, 1, a, ip, w);
734 for (j = 0; j <= n - 1; j++) {
735 a[j] *= 2.0 / n;
736 }
737 .
738
739
740 -------- Cosine Transform of RDFT (Real Symmetric DFT) --------
741 [definition]
742 C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
743 [usage]
744 ip[0] = 0; // first time only
745 dfct(n, a, t, ip, w);
746 [parameters]
747 n :data length - 1 (int)
748 n >= 2, n = power of 2
749 a[0...n] :input/output data (REAL *)
750 output data
751 a[k] = C[k], 0<=k<=n
752 t[0...n/2] :work area (REAL *)
753 ip[0...*] :work area for bit reversal (int *)
754 length of ip >= 2+sqrt(n/4)
755 strictly,
756 length of ip >=
757 2+(1<<(int)(log(n/4+0.5)/log(2))/2).
758 ip[0],ip[1] are pointers of the cos/sin table.
759 w[0...n*5/8-1] :cos/sin table (REAL *)
760 w[],ip[] are initialized if ip[0] == 0.
761 [remark]
762 Inverse of
763 a[0] *= 0.5;
764 a[n] *= 0.5;
765 dfct(n, a, t, ip, w);
766 is
767 a[0] *= 0.5;
768 a[n] *= 0.5;
769 dfct(n, a, t, ip, w);
770 for (j = 0; j <= n; j++) {
771 a[j] *= 2.0 / n;
772 }
773 .
774
775
776 -------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
777 [definition]
778 S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
779 [usage]
780 ip[0] = 0; // first time only
781 dfst(n, a, t, ip, w);
782 [parameters]
783 n :data length + 1 (int)
784 n >= 2, n = power of 2
785 a[0...n-1] :input/output data (REAL *)
786 output data
787 a[k] = S[k], 0<k<n
788 (a[0] is used for work area)
789 t[0...n/2-1] :work area (REAL *)
790 ip[0...*] :work area for bit reversal (int *)
791 length of ip >= 2+sqrt(n/4)
792 strictly,
793 length of ip >=
794 2+(1<<(int)(log(n/4+0.5)/log(2))/2).
795 ip[0],ip[1] are pointers of the cos/sin table.
796 w[0...n*5/8-1] :cos/sin table (REAL *)
797 w[],ip[] are initialized if ip[0] == 0.
798 [remark]
799 Inverse of
800 dfst(n, a, t, ip, w);
801 is
802 dfst(n, a, t, ip, w);
803 for (j = 1; j <= n - 1; j++) {
804 a[j] *= 2.0 / n;
805 }
806 .
807
808
809 Appendix :
810 The cos/sin table is recalculated when the larger table required.
811 w[] and ip[] are compatible with all routines.
812 */
813
814
815 // Complex Discrete Fourier Transform
816 /*public*/ void cdft (int n, int isgn, REAL* a, int* ip, REAL* w) {
817 int nw = ip[0];
818 if (n > (nw << 2)) {
819 nw = n >> 2;
820 makewt(nw, ip, w);
821 }
822 if (isgn >= 0) {
823 cftfsub(n, a, ip, nw, w);
824 } else {
825 cftbsub(n, a, ip, nw, w);
826 }
827 }
828
829
830 // Real Discrete Fourier Transform
831 /*public*/ void rdft (int n, int isgn, REAL* a, int* ip, REAL* w) {
832 int nw = ip[0];
833 if (n > (nw << 2)) {
834 nw = n >> 2;
835 makewt(nw, ip, w);
836 }
837 int nc = ip[1];
838 if (n > (nc << 2)) {
839 nc = n >> 2;
840 makect(nc, ip, w + nw);
841 }
842 if (isgn >= 0) {
843 if (n > 4) {
844 cftfsub(n, a, ip, nw, w);
845 rftfsub(n, a, nc, w + nw);
846 } else if (n == 4) {
847 cftfsub(n, a, ip, nw, w);
848 }
849 REAL xi = a[0] - a[1];
850 a[0] += a[1];
851 a[1] = xi;
852 } else {
853 a[1] = 0.5 * (a[0] - a[1]);
854 a[0] -= a[1];
855 if (n > 4) {
856 rftbsub(n, a, nc, w + nw);
857 cftbsub(n, a, ip, nw, w);
858 } else if (n == 4) {
859 cftbsub(n, a, ip, nw, w);
860 }
861 }
862 }
863
864
865 // Discrete Cosine Transform
866 /*public*/ void ddct (int n, int isgn, REAL* a, int* ip, REAL* w) {
867 int nw = ip[0];
868 if (n > (nw << 2)) {
869 nw = n >> 2;
870 makewt(nw, ip, w);
871 }
872 int nc = ip[1];
873 if (n > nc) {
874 nc = n;
875 makect(nc, ip, w + nw);
876 }
877 if (isgn < 0) {
878 REAL xr = a[n - 1];
879 for (int j = n - 2; j >= 2; j -= 2) {
880 a[j + 1] = a[j] - a[j - 1];
881 a[j] += a[j - 1];
882 }
883 a[1] = a[0] - xr;
884 a[0] += xr;
885 if (n > 4) {
886 rftbsub(n, a, nc, w + nw);
887 cftbsub(n, a, ip, nw, w);
888 } else if (n == 4) {
889 cftbsub(n, a, ip, nw, w);
890 }
891 }
892 dctsub(n, a, nc, w + nw);
893 if (isgn >= 0) {
894 if (n > 4) {
895 cftfsub(n, a, ip, nw, w);
896 rftfsub(n, a, nc, w + nw);
897 } else if (n == 4) {
898 cftfsub(n, a, ip, nw, w);
899 }
900 REAL xr = a[0] - a[1];
901 a[0] += a[1];
902 for (int j = 2; j < n; j += 2) {
903 a[j - 1] = a[j] - a[j + 1];
904 a[j] += a[j + 1];
905 }
906 a[n - 1] = xr;
907 }
908 }
909
910
911 // Discrete Sine Transform
912 /*public*/ void ddst (int n, int isgn, REAL* a, int* ip, REAL* w) {
913 int nw = ip[0];
914 if (n > (nw << 2)) {
915 nw = n >> 2;
916 makewt(nw, ip, w);
917 }
918 int nc = ip[1];
919 if (n > nc) {
920 nc = n;
921 makect(nc, ip, w + nw);
922 }
923 if (isgn < 0) {
924 REAL xr = a[n - 1];
925 for (int j = n - 2; j >= 2; j -= 2) {
926 a[j + 1] = -a[j] - a[j - 1];
927 a[j] -= a[j - 1];
928 }
929 a[1] = a[0] + xr;
930 a[0] -= xr;
931 if (n > 4) {
932 rftbsub(n, a, nc, w + nw);
933 cftbsub(n, a, ip, nw, w);
934 } else if (n == 4) {
935 cftbsub(n, a, ip, nw, w);
936 }
937 }
938 dstsub(n, a, nc, w + nw);
939 if (isgn >= 0) {
940 if (n > 4) {
941 cftfsub(n, a, ip, nw, w);
942 rftfsub(n, a, nc, w + nw);
943 } else if (n == 4) {
944 cftfsub(n, a, ip, nw, w);
945 }
946 REAL xr = a[0] - a[1];
947 a[0] += a[1];
948 for (int j = 2; j < n; j += 2) {
949 a[j - 1] = -a[j] - a[j + 1];
950 a[j] -= a[j + 1];
951 }
952 a[n - 1] = -xr;
953 }
954 }
955
956
957 // Cosine Transform of RDFT (Real Symmetric DFT)
958 /*public*/ void dfct (int n, REAL* a, REAL* t, int* ip, REAL* w) {
959 int j, k, l, m, mh, nw, nc;
960 //REAL xr, xi, yr, yi;
961
962 nw = ip[0];
963 if (n > (nw << 3)) {
964 nw = n >> 3;
965 makewt(nw, ip, w);
966 }
967 nc = ip[1];
968 if (n > (nc << 1)) {
969 nc = n >> 1;
970 makect(nc, ip, w + nw);
971 }
972 m = n >> 1;
973 REAL yi = a[m];
974 REAL xi = a[0] + a[n];
975 a[0] -= a[n];
976 t[0] = xi - yi;
977 t[m] = xi + yi;
978 if (n > 2) {
979 mh = m >> 1;
980 for (j = 1; j < mh; j++) {
981 k = m - j;
982 REAL xr = a[j] - a[n - j];
983 xi = a[j] + a[n - j];
984 REAL yr = a[k] - a[n - k];
985 yi = a[k] + a[n - k];
986 a[j] = xr;
987 a[k] = yr;
988 t[j] = xi - yi;
989 t[k] = xi + yi;
990 }
991 t[mh] = a[mh] + a[n - mh];
992 a[mh] -= a[n - mh];
993 dctsub(m, a, nc, w + nw);
994 if (m > 4) {
995 cftfsub(m, a, ip, nw, w);
996 rftfsub(m, a, nc, w + nw);
997 } else if (m == 4) {
998 cftfsub(m, a, ip, nw, w);
999 }
1000 a[n - 1] = a[0] - a[1];
1001 a[1] = a[0] + a[1];
1002 for (j = m - 2; j >= 2; j -= 2) {
1003 a[2 * j + 1] = a[j] + a[j + 1];
1004 a[2 * j - 1] = a[j] - a[j + 1];
1005 }
1006 l = 2;
1007 m = mh;
1008 while (m >= 2) {
1009 dctsub(m, t, nc, w + nw);
1010 if (m > 4) {
1011 cftfsub(m, t, ip, nw, w);
1012 rftfsub(m, t, nc, w + nw);
1013 } else if (m == 4) {
1014 cftfsub(m, t, ip, nw, w);
1015 }
1016 a[n - l] = t[0] - t[1];
1017 a[l] = t[0] + t[1];
1018 k = 0;
1019 for (j = 2; j < m; j += 2) {
1020 k += l << 2;
1021 a[k - l] = t[j] - t[j + 1];
1022 a[k + l] = t[j] + t[j + 1];
1023 }
1024 l <<= 1;
1025 mh = m >> 1;
1026 for (j = 0; j < mh; j++) {
1027 k = m - j;
1028 t[j] = t[m + k] - t[m + j];
1029 t[k] = t[m + k] + t[m + j];
1030 }
1031 t[mh] = t[m + mh];
1032 m = mh;
1033 }
1034 a[l] = t[0];
1035 a[n] = t[2] - t[1];
1036 a[0] = t[2] + t[1];
1037 } else {
1038 a[1] = a[0];
1039 a[2] = t[0];
1040 a[0] = t[1];
1041 }
1042 }
1043
1044
1045 // Sine Transform of RDFT (Real Anti-symmetric DFT)
1046 /*public*/ void dfst (int n, REAL* a, REAL* t, int* ip, REAL* w) {
1047 int j, k, l, m, mh, nw, nc;
1048 //REAL xr, xi, yr, yi;
1049
1050 nw = ip[0];
1051 if (n > (nw << 3)) {
1052 nw = n >> 3;
1053 makewt(nw, ip, w);
1054 }
1055 nc = ip[1];
1056 if (n > (nc << 1)) {
1057 nc = n >> 1;
1058 makect(nc, ip, w + nw);
1059 }
1060 if (n > 2) {
1061 m = n >> 1;
1062 mh = m >> 1;
1063 for (j = 1; j < mh; j++) {
1064 k = m - j;
1065 REAL xr = a[j] + a[n - j];
1066 REAL xi = a[j] - a[n - j];
1067 REAL yr = a[k] + a[n - k];
1068 REAL yi = a[k] - a[n - k];
1069 a[j] = xr;
1070 a[k] = yr;
1071 t[j] = xi + yi;
1072 t[k] = xi - yi;
1073 }
1074 t[0] = a[mh] - a[n - mh];
1075 a[mh] += a[n - mh];
1076 a[0] = a[m];
1077 dstsub(m, a, nc, w + nw);
1078 if (m > 4) {
1079 cftfsub(m, a, ip, nw, w);
1080 rftfsub(m, a, nc, w + nw);
1081 } else if (m == 4) {
1082 cftfsub(m, a, ip, nw, w);
1083 }
1084 a[n - 1] = a[1] - a[0];
1085 a[1] = a[0] + a[1];
1086 for (j = m - 2; j >= 2; j -= 2) {
1087 a[2 * j + 1] = a[j] - a[j + 1];
1088 a[2 * j - 1] = -a[j] - a[j + 1];
1089 }
1090 l = 2;
1091 m = mh;
1092 while (m >= 2) {
1093 dstsub(m, t, nc, w + nw);
1094 if (m > 4) {
1095 cftfsub(m, t, ip, nw, w);
1096 rftfsub(m, t, nc, w + nw);
1097 } else if (m == 4) {
1098 cftfsub(m, t, ip, nw, w);
1099 }
1100 a[n - l] = t[1] - t[0];
1101 a[l] = t[0] + t[1];
1102 k = 0;
1103 for (j = 2; j < m; j += 2) {
1104 k += l << 2;
1105 a[k - l] = -t[j] - t[j + 1];
1106 a[k + l] = t[j] - t[j + 1];
1107 }
1108 l <<= 1;
1109 mh = m >> 1;
1110 for (j = 1; j < mh; j++) {
1111 k = m - j;
1112 t[j] = t[m + k] + t[m + j];
1113 t[k] = t[m + k] - t[m + j];
1114 }
1115 t[0] = t[m + mh];
1116 m = mh;
1117 }
1118 a[l] = t[0];
1119 }
1120 a[0] = 0;
1121 }
1122
1123
1124 // ////////////////////////////////////////////////////////////////////////// //
1125 private:
1126 /* -------- initializing routines -------- */
1127
1128 void makewt(int nw, int *ip, REAL *w)
1129 {
1130 import std.math : atan, cos, sin;
1131 int j, nwh, nw0, nw1;
1132 //REAL delta, wn4r, wk1r, wk1i, wk3r, wk3i;
1133
1134 ip[0] = nw;
1135 ip[1] = 1;
1136 if (nw > 2) {
1137 nwh = nw >> 1;
1138 REAL delta = atan(1.0) / nwh;
1139 REAL wn4r = cos(delta * nwh);
1140 w[0] = 1;
1141 w[1] = wn4r;
1142 if (nwh == 4) {
1143 w[2] = cos(delta * 2);
1144 w[3] = sin(delta * 2);
1145 } else if (nwh > 4) {
1146 makeipt(nw, ip);
1147 w[2] = 0.5 / cos(delta * 2);
1148 w[3] = 0.5 / cos(delta * 6);
1149 for (j = 4; j < nwh; j += 4) {
1150 w[j] = cos(delta * j);
1151 w[j + 1] = sin(delta * j);
1152 w[j + 2] = cos(3 * delta * j);
1153 w[j + 3] = -sin(3 * delta * j);
1154 }
1155 }
1156 nw0 = 0;
1157 while (nwh > 2) {
1158 nw1 = nw0 + nwh;
1159 nwh >>= 1;
1160 w[nw1] = 1;
1161 w[nw1 + 1] = wn4r;
1162 if (nwh == 4) {
1163 REAL wk1r = w[nw0 + 4];
1164 REAL wk1i = w[nw0 + 5];
1165 w[nw1 + 2] = wk1r;
1166 w[nw1 + 3] = wk1i;
1167 } else if (nwh > 4) {
1168 REAL wk1r = w[nw0 + 4];
1169 REAL wk9r = w[nw0 + 6];
1170 w[nw1 + 2] = 0.5 / wk1r;
1171 w[nw1 + 3] = 0.5 / wk9r;
1172 for (j = 4; j < nwh; j += 4) {
1173 wk1r = w[nw0 + 2 * j];
1174 REAL wk1i = w[nw0 + 2 * j + 1];
1175 REAL wk3r = w[nw0 + 2 * j + 2];
1176 REAL wk3i = w[nw0 + 2 * j + 3];
1177 w[nw1 + j] = wk1r;
1178 w[nw1 + j + 1] = wk1i;
1179 w[nw1 + j + 2] = wk3r;
1180 w[nw1 + j + 3] = wk3i;
1181 }
1182 }
1183 nw0 = nw1;
1184 }
1185 }
1186 }
1187
1188
1189 void makeipt(int nw, int *ip)
1190 {
1191 int j, l, m, m2, p, q;
1192
1193 ip[2] = 0;
1194 ip[3] = 16;
1195 m = 2;
1196 for (l = nw; l > 32; l >>= 2) {
1197 m2 = m << 1;
1198 q = m2 << 3;
1199 for (j = m; j < m2; j++) {
1200 p = ip[j] << 2;
1201 ip[m + j] = p;
1202 ip[m2 + j] = p + q;
1203 }
1204 m = m2;
1205 }
1206 }
1207
1208
1209 void makect(int nc, int *ip, REAL *c)
1210 {
1211 import std.math : atan, cos, sin;
1212 int j, nch;
1213 //REAL delta;
1214
1215 ip[1] = nc;
1216 if (nc > 1) {
1217 nch = nc >> 1;
1218 REAL delta = atan(1.0) / nch;
1219 c[0] = cos(delta * nch);
1220 c[nch] = 0.5 * c[0];
1221 for (j = 1; j < nch; j++) {
1222 c[j] = 0.5 * cos(delta * j);
1223 c[nc - j] = 0.5 * sin(delta * j);
1224 }
1225 }
1226 }
1227
1228
1229 /* -------- child routines -------- */
1230
1231 void cftfsub(int n, REAL *a, int *ip, int nw, REAL *w)
1232 {
1233
1234 if (n > 8) {
1235 if (n > 32) {
1236 cftf1st(n, a, &w[nw - (n >> 2)]);
1237 if (n > 512) {
1238 cftrec4(n, a, nw, w);
1239 } else if (n > 128) {
1240 cftleaf(n, 1, a, nw, w);
1241 } else {
1242 cftfx41(n, a, nw, w);
1243 }
1244 bitrv2(n, ip, a);
1245 } else if (n == 32) {
1246 cftf161(a, &w[nw - 8]);
1247 bitrv216(a);
1248 } else {
1249 cftf081(a, w);
1250 bitrv208(a);
1251 }
1252 } else if (n == 8) {
1253 cftf040(a);
1254 } else if (n == 4) {
1255 cftx020(a);
1256 }
1257 }
1258
1259
1260 void cftbsub(int n, REAL *a, int *ip, int nw, REAL *w)
1261 {
1262 if (n > 8) {
1263 if (n > 32) {
1264 cftb1st(n, a, &w[nw - (n >> 2)]);
1265 if (n > 512) {
1266 cftrec4(n, a, nw, w);
1267 } else if (n > 128) {
1268 cftleaf(n, 1, a, nw, w);
1269 } else {
1270 cftfx41(n, a, nw, w);
1271 }
1272 bitrv2conj(n, ip, a);
1273 } else if (n == 32) {
1274 cftf161(a, &w[nw - 8]);
1275 bitrv216neg(a);
1276 } else {
1277 cftf081(a, w);
1278 bitrv208neg(a);
1279 }
1280 } else if (n == 8) {
1281 cftb040(a);
1282 } else if (n == 4) {
1283 cftx020(a);
1284 }
1285 }
1286
1287
1288 void bitrv2(int n, int *ip, REAL *a)
1289 {
1290 int j, j1, k, k1, l, m, nh, nm;
1291 //REAL xr, xi, yr, yi;
1292
1293 m = 1;
1294 for (l = n >> 2; l > 8; l >>= 2) {
1295 m <<= 1;
1296 }
1297 nh = n >> 1;
1298 nm = 4 * m;
1299 if (l == 8) {
1300 for (k = 0; k < m; k++) {
1301 for (j = 0; j < k; j++) {
1302 j1 = 4 * j + 2 * ip[m + k];
1303 k1 = 4 * k + 2 * ip[m + j];
1304 REAL xr = a[j1];
1305 REAL xi = a[j1 + 1];
1306 REAL yr = a[k1];
1307 REAL yi = a[k1 + 1];
1308 a[j1] = yr;
1309 a[j1 + 1] = yi;
1310 a[k1] = xr;
1311 a[k1 + 1] = xi;
1312 j1 += nm;
1313 k1 += 2 * nm;
1314 xr = a[j1];
1315 xi = a[j1 + 1];
1316 yr = a[k1];
1317 yi = a[k1 + 1];
1318 a[j1] = yr;
1319 a[j1 + 1] = yi;
1320 a[k1] = xr;
1321 a[k1 + 1] = xi;
1322 j1 += nm;
1323 k1 -= nm;
1324 xr = a[j1];
1325 xi = a[j1 + 1];
1326 yr = a[k1];
1327 yi = a[k1 + 1];
1328 a[j1] = yr;
1329 a[j1 + 1] = yi;
1330 a[k1] = xr;
1331 a[k1 + 1] = xi;
1332 j1 += nm;
1333 k1 += 2 * nm;
1334 xr = a[j1];
1335 xi = a[j1 + 1];
1336 yr = a[k1];
1337 yi = a[k1 + 1];
1338 a[j1] = yr;
1339 a[j1 + 1] = yi;
1340 a[k1] = xr;
1341 a[k1 + 1] = xi;
1342 j1 += nh;
1343 k1 += 2;
1344 xr = a[j1];
1345 xi = a[j1 + 1];
1346 yr = a[k1];
1347 yi = a[k1 + 1];
1348 a[j1] = yr;
1349 a[j1 + 1] = yi;
1350 a[k1] = xr;
1351 a[k1 + 1] = xi;
1352 j1 -= nm;
1353 k1 -= 2 * nm;
1354 xr = a[j1];
1355 xi = a[j1 + 1];
1356 yr = a[k1];
1357 yi = a[k1 + 1];
1358 a[j1] = yr;
1359 a[j1 + 1] = yi;
1360 a[k1] = xr;
1361 a[k1 + 1] = xi;
1362 j1 -= nm;
1363 k1 += nm;
1364 xr = a[j1];
1365 xi = a[j1 + 1];
1366 yr = a[k1];
1367 yi = a[k1 + 1];
1368 a[j1] = yr;
1369 a[j1 + 1] = yi;
1370 a[k1] = xr;
1371 a[k1 + 1] = xi;
1372 j1 -= nm;
1373 k1 -= 2 * nm;
1374 xr = a[j1];
1375 xi = a[j1 + 1];
1376 yr = a[k1];
1377 yi = a[k1 + 1];
1378 a[j1] = yr;
1379 a[j1 + 1] = yi;
1380 a[k1] = xr;
1381 a[k1 + 1] = xi;
1382 j1 += 2;
1383 k1 += nh;
1384 xr = a[j1];
1385 xi = a[j1 + 1];
1386 yr = a[k1];
1387 yi = a[k1 + 1];
1388 a[j1] = yr;
1389 a[j1 + 1] = yi;
1390 a[k1] = xr;
1391 a[k1 + 1] = xi;
1392 j1 += nm;
1393 k1 += 2 * nm;
1394 xr = a[j1];
1395 xi = a[j1 + 1];
1396 yr = a[k1];
1397 yi = a[k1 + 1];
1398 a[j1] = yr;
1399 a[j1 + 1] = yi;
1400 a[k1] = xr;
1401 a[k1 + 1] = xi;
1402 j1 += nm;
1403 k1 -= nm;
1404 xr = a[j1];
1405 xi = a[j1 + 1];
1406 yr = a[k1];
1407 yi = a[k1 + 1];
1408 a[j1] = yr;
1409 a[j1 + 1] = yi;
1410 a[k1] = xr;
1411 a[k1 + 1] = xi;
1412 j1 += nm;
1413 k1 += 2 * nm;
1414 xr = a[j1];
1415 xi = a[j1 + 1];
1416 yr = a[k1];
1417 yi = a[k1 + 1];
1418 a[j1] = yr;
1419 a[j1 + 1] = yi;
1420 a[k1] = xr;
1421 a[k1 + 1] = xi;
1422 j1 -= nh;
1423 k1 -= 2;
1424 xr = a[j1];
1425 xi = a[j1 + 1];
1426 yr = a[k1];
1427 yi = a[k1 + 1];
1428 a[j1] = yr;
1429 a[j1 + 1] = yi;
1430 a[k1] = xr;
1431 a[k1 + 1] = xi;
1432 j1 -= nm;
1433 k1 -= 2 * nm;
1434 xr = a[j1];
1435 xi = a[j1 + 1];
1436 yr = a[k1];
1437 yi = a[k1 + 1];
1438 a[j1] = yr;
1439 a[j1 + 1] = yi;
1440 a[k1] = xr;
1441 a[k1 + 1] = xi;
1442 j1 -= nm;
1443 k1 += nm;
1444 xr = a[j1];
1445 xi = a[j1 + 1];
1446 yr = a[k1];
1447 yi = a[k1 + 1];
1448 a[j1] = yr;
1449 a[j1 + 1] = yi;
1450 a[k1] = xr;
1451 a[k1 + 1] = xi;
1452 j1 -= nm;
1453 k1 -= 2 * nm;
1454 xr = a[j1];
1455 xi = a[j1 + 1];
1456 yr = a[k1];
1457 yi = a[k1 + 1];
1458 a[j1] = yr;
1459 a[j1 + 1] = yi;
1460 a[k1] = xr;
1461 a[k1 + 1] = xi;
1462 }
1463 k1 = 4 * k + 2 * ip[m + k];
1464 j1 = k1 + 2;
1465 k1 += nh;
1466 REAL xr = a[j1];
1467 REAL xi = a[j1 + 1];
1468 REAL yr = a[k1];
1469 REAL yi = a[k1 + 1];
1470 a[j1] = yr;
1471 a[j1 + 1] = yi;
1472 a[k1] = xr;
1473 a[k1 + 1] = xi;
1474 j1 += nm;
1475 k1 += 2 * nm;
1476 xr = a[j1];
1477 xi = a[j1 + 1];
1478 yr = a[k1];
1479 yi = a[k1 + 1];
1480 a[j1] = yr;
1481 a[j1 + 1] = yi;
1482 a[k1] = xr;
1483 a[k1 + 1] = xi;
1484 j1 += nm;
1485 k1 -= nm;
1486 xr = a[j1];
1487 xi = a[j1 + 1];
1488 yr = a[k1];
1489 yi = a[k1 + 1];
1490 a[j1] = yr;
1491 a[j1 + 1] = yi;
1492 a[k1] = xr;
1493 a[k1 + 1] = xi;
1494 j1 -= 2;
1495 k1 -= nh;
1496 xr = a[j1];
1497 xi = a[j1 + 1];
1498 yr = a[k1];
1499 yi = a[k1 + 1];
1500 a[j1] = yr;
1501 a[j1 + 1] = yi;
1502 a[k1] = xr;
1503 a[k1 + 1] = xi;
1504 j1 += nh + 2;
1505 k1 += nh + 2;
1506 xr = a[j1];
1507 xi = a[j1 + 1];
1508 yr = a[k1];
1509 yi = a[k1 + 1];
1510 a[j1] = yr;
1511 a[j1 + 1] = yi;
1512 a[k1] = xr;
1513 a[k1 + 1] = xi;
1514 j1 -= nh - nm;
1515 k1 += 2 * nm - 2;
1516 xr = a[j1];
1517 xi = a[j1 + 1];
1518 yr = a[k1];
1519 yi = a[k1 + 1];
1520 a[j1] = yr;
1521 a[j1 + 1] = yi;
1522 a[k1] = xr;
1523 a[k1 + 1] = xi;
1524 }
1525 } else {
1526 for (k = 0; k < m; k++) {
1527 for (j = 0; j < k; j++) {
1528 j1 = 4 * j + ip[m + k];
1529 k1 = 4 * k + ip[m + j];
1530 REAL xr = a[j1];
1531 REAL xi = a[j1 + 1];
1532 REAL yr = a[k1];
1533 REAL yi = a[k1 + 1];
1534 a[j1] = yr;
1535 a[j1 + 1] = yi;
1536 a[k1] = xr;
1537 a[k1 + 1] = xi;
1538 j1 += nm;
1539 k1 += nm;
1540 xr = a[j1];
1541 xi = a[j1 + 1];
1542 yr = a[k1];
1543 yi = a[k1 + 1];
1544 a[j1] = yr;
1545 a[j1 + 1] = yi;
1546 a[k1] = xr;
1547 a[k1 + 1] = xi;
1548 j1 += nh;
1549 k1 += 2;
1550 xr = a[j1];
1551 xi = a[j1 + 1];
1552 yr = a[k1];
1553 yi = a[k1 + 1];
1554 a[j1] = yr;
1555 a[j1 + 1] = yi;
1556 a[k1] = xr;
1557 a[k1 + 1] = xi;
1558 j1 -= nm;
1559 k1 -= nm;
1560 xr = a[j1];
1561 xi = a[j1 + 1];
1562 yr = a[k1];
1563 yi = a[k1 + 1];
1564 a[j1] = yr;
1565 a[j1 + 1] = yi;
1566 a[k1] = xr;
1567 a[k1 + 1] = xi;
1568 j1 += 2;
1569 k1 += nh;
1570 xr = a[j1];
1571 xi = a[j1 + 1];
1572 yr = a[k1];
1573 yi = a[k1 + 1];
1574 a[j1] = yr;
1575 a[j1 + 1] = yi;
1576 a[k1] = xr;
1577 a[k1 + 1] = xi;
1578 j1 += nm;
1579 k1 += nm;
1580 xr = a[j1];
1581 xi = a[j1 + 1];
1582 yr = a[k1];
1583 yi = a[k1 + 1];
1584 a[j1] = yr;
1585 a[j1 + 1] = yi;
1586 a[k1] = xr;
1587 a[k1 + 1] = xi;
1588 j1 -= nh;
1589 k1 -= 2;
1590 xr = a[j1];
1591 xi = a[j1 + 1];
1592 yr = a[k1];
1593 yi = a[k1 + 1];
1594 a[j1] = yr;
1595 a[j1 + 1] = yi;
1596 a[k1] = xr;
1597 a[k1 + 1] = xi;
1598 j1 -= nm;
1599 k1 -= nm;
1600 xr = a[j1];
1601 xi = a[j1 + 1];
1602 yr = a[k1];
1603 yi = a[k1 + 1];
1604 a[j1] = yr;
1605 a[j1 + 1] = yi;
1606 a[k1] = xr;
1607 a[k1 + 1] = xi;
1608 }
1609 k1 = 4 * k + ip[m + k];
1610 j1 = k1 + 2;
1611 k1 += nh;
1612 REAL xr = a[j1];
1613 REAL xi = a[j1 + 1];
1614 REAL yr = a[k1];
1615 REAL yi = a[k1 + 1];
1616 a[j1] = yr;
1617 a[j1 + 1] = yi;
1618 a[k1] = xr;
1619 a[k1 + 1] = xi;
1620 j1 += nm;
1621 k1 += nm;
1622 xr = a[j1];
1623 xi = a[j1 + 1];
1624 yr = a[k1];
1625 yi = a[k1 + 1];
1626 a[j1] = yr;
1627 a[j1 + 1] = yi;
1628 a[k1] = xr;
1629 a[k1 + 1] = xi;
1630 }
1631 }
1632 }
1633
1634
1635 void bitrv2conj(int n, int *ip, REAL *a)
1636 {
1637 int j, j1, k, k1, l, m, nh, nm;
1638 //REAL xr, xi, yr, yi;
1639
1640 m = 1;
1641 for (l = n >> 2; l > 8; l >>= 2) {
1642 m <<= 1;
1643 }
1644 nh = n >> 1;
1645 nm = 4 * m;
1646 if (l == 8) {
1647 for (k = 0; k < m; k++) {
1648 for (j = 0; j < k; j++) {
1649 j1 = 4 * j + 2 * ip[m + k];
1650 k1 = 4 * k + 2 * ip[m + j];
1651 REAL xr = a[j1];
1652 REAL xi = -a[j1 + 1];
1653 REAL yr = a[k1];
1654 REAL yi = -a[k1 + 1];
1655 a[j1] = yr;
1656 a[j1 + 1] = yi;
1657 a[k1] = xr;
1658 a[k1 + 1] = xi;
1659 j1 += nm;
1660 k1 += 2 * nm;
1661 xr = a[j1];
1662 xi = -a[j1 + 1];
1663 yr = a[k1];
1664 yi = -a[k1 + 1];
1665 a[j1] = yr;
1666 a[j1 + 1] = yi;
1667 a[k1] = xr;
1668 a[k1 + 1] = xi;
1669 j1 += nm;
1670 k1 -= nm;
1671 xr = a[j1];
1672 xi = -a[j1 + 1];
1673 yr = a[k1];
1674 yi = -a[k1 + 1];
1675 a[j1] = yr;
1676 a[j1 + 1] = yi;
1677 a[k1] = xr;
1678 a[k1 + 1] = xi;
1679 j1 += nm;
1680 k1 += 2 * nm;
1681 xr = a[j1];
1682 xi = -a[j1 + 1];
1683 yr = a[k1];
1684 yi = -a[k1 + 1];
1685 a[j1] = yr;
1686 a[j1 + 1] = yi;
1687 a[k1] = xr;
1688 a[k1 + 1] = xi;
1689 j1 += nh;
1690 k1 += 2;
1691 xr = a[j1];
1692 xi = -a[j1 + 1];
1693 yr = a[k1];
1694 yi = -a[k1 + 1];
1695 a[j1] = yr;
1696 a[j1 + 1] = yi;
1697 a[k1] = xr;
1698 a[k1 + 1] = xi;
1699 j1 -= nm;
1700 k1 -= 2 * nm;
1701 xr = a[j1];
1702 xi = -a[j1 + 1];
1703 yr = a[k1];
1704 yi = -a[k1 + 1];
1705 a[j1] = yr;
1706 a[j1 + 1] = yi;
1707 a[k1] = xr;
1708 a[k1 + 1] = xi;
1709 j1 -= nm;
1710 k1 += nm;
1711 xr = a[j1];
1712 xi = -a[j1 + 1];
1713 yr = a[k1];
1714 yi = -a[k1 + 1];
1715 a[j1] = yr;
1716 a[j1 + 1] = yi;
1717 a[k1] = xr;
1718 a[k1 + 1] = xi;
1719 j1 -= nm;
1720 k1 -= 2 * nm;
1721 xr = a[j1];
1722 xi = -a[j1 + 1];
1723 yr = a[k1];
1724 yi = -a[k1 + 1];
1725 a[j1] = yr;
1726 a[j1 + 1] = yi;
1727 a[k1] = xr;
1728 a[k1 + 1] = xi;
1729 j1 += 2;
1730 k1 += nh;
1731 xr = a[j1];
1732 xi = -a[j1 + 1];
1733 yr = a[k1];
1734 yi = -a[k1 + 1];
1735 a[j1] = yr;
1736 a[j1 + 1] = yi;
1737 a[k1] = xr;
1738 a[k1 + 1] = xi;
1739 j1 += nm;
1740 k1 += 2 * nm;
1741 xr = a[j1];
1742 xi = -a[j1 + 1];
1743 yr = a[k1];
1744 yi = -a[k1 + 1];
1745 a[j1] = yr;
1746 a[j1 + 1] = yi;
1747 a[k1] = xr;
1748 a[k1 + 1] = xi;
1749 j1 += nm;
1750 k1 -= nm;
1751 xr = a[j1];
1752 xi = -a[j1 + 1];
1753 yr = a[k1];
1754 yi = -a[k1 + 1];
1755 a[j1] = yr;
1756 a[j1 + 1] = yi;
1757 a[k1] = xr;
1758 a[k1 + 1] = xi;
1759 j1 += nm;
1760 k1 += 2 * nm;
1761 xr = a[j1];
1762 xi = -a[j1 + 1];
1763 yr = a[k1];
1764 yi = -a[k1 + 1];
1765 a[j1] = yr;
1766 a[j1 + 1] = yi;
1767 a[k1] = xr;
1768 a[k1 + 1] = xi;
1769 j1 -= nh;
1770 k1 -= 2;
1771 xr = a[j1];
1772 xi = -a[j1 + 1];
1773 yr = a[k1];
1774 yi = -a[k1 + 1];
1775 a[j1] = yr;
1776 a[j1 + 1] = yi;
1777 a[k1] = xr;
1778 a[k1 + 1] = xi;
1779 j1 -= nm;
1780 k1 -= 2 * nm;
1781 xr = a[j1];
1782 xi = -a[j1 + 1];
1783 yr = a[k1];
1784 yi = -a[k1 + 1];
1785 a[j1] = yr;
1786 a[j1 + 1] = yi;
1787 a[k1] = xr;
1788 a[k1 + 1] = xi;
1789 j1 -= nm;
1790 k1 += nm;
1791 xr = a[j1];
1792 xi = -a[j1 + 1];
1793 yr = a[k1];
1794 yi = -a[k1 + 1];
1795 a[j1] = yr;
1796 a[j1 + 1] = yi;
1797 a[k1] = xr;
1798 a[k1 + 1] = xi;
1799 j1 -= nm;
1800 k1 -= 2 * nm;
1801 xr = a[j1];
1802 xi = -a[j1 + 1];
1803 yr = a[k1];
1804 yi = -a[k1 + 1];
1805 a[j1] = yr;
1806 a[j1 + 1] = yi;
1807 a[k1] = xr;
1808 a[k1 + 1] = xi;
1809 }
1810 k1 = 4 * k + 2 * ip[m + k];
1811 j1 = k1 + 2;
1812 k1 += nh;
1813 a[j1 - 1] = -a[j1 - 1];
1814 REAL xr = a[j1];
1815 REAL xi = -a[j1 + 1];
1816 REAL yr = a[k1];
1817 REAL yi = -a[k1 + 1];
1818 a[j1] = yr;
1819 a[j1 + 1] = yi;
1820 a[k1] = xr;
1821 a[k1 + 1] = xi;
1822 a[k1 + 3] = -a[k1 + 3];
1823 j1 += nm;
1824 k1 += 2 * nm;
1825 xr = a[j1];
1826 xi = -a[j1 + 1];
1827 yr = a[k1];
1828 yi = -a[k1 + 1];
1829 a[j1] = yr;
1830 a[j1 + 1] = yi;
1831 a[k1] = xr;
1832 a[k1 + 1] = xi;
1833 j1 += nm;
1834 k1 -= nm;
1835 xr = a[j1];
1836 xi = -a[j1 + 1];
1837 yr = a[k1];
1838 yi = -a[k1 + 1];
1839 a[j1] = yr;
1840 a[j1 + 1] = yi;
1841 a[k1] = xr;
1842 a[k1 + 1] = xi;
1843 j1 -= 2;
1844 k1 -= nh;
1845 xr = a[j1];
1846 xi = -a[j1 + 1];
1847 yr = a[k1];
1848 yi = -a[k1 + 1];
1849 a[j1] = yr;
1850 a[j1 + 1] = yi;
1851 a[k1] = xr;
1852 a[k1 + 1] = xi;
1853 j1 += nh + 2;
1854 k1 += nh + 2;
1855 xr = a[j1];
1856 xi = -a[j1 + 1];
1857 yr = a[k1];
1858 yi = -a[k1 + 1];
1859 a[j1] = yr;
1860 a[j1 + 1] = yi;
1861 a[k1] = xr;
1862 a[k1 + 1] = xi;
1863 j1 -= nh - nm;
1864 k1 += 2 * nm - 2;
1865 a[j1 - 1] = -a[j1 - 1];
1866 xr = a[j1];
1867 xi = -a[j1 + 1];
1868 yr = a[k1];
1869 yi = -a[k1 + 1];
1870 a[j1] = yr;
1871 a[j1 + 1] = yi;
1872 a[k1] = xr;
1873 a[k1 + 1] = xi;
1874 a[k1 + 3] = -a[k1 + 3];
1875 }
1876 } else {
1877 for (k = 0; k < m; k++) {
1878 for (j = 0; j < k; j++) {
1879 j1 = 4 * j + ip[m + k];
1880 k1 = 4 * k + ip[m + j];
1881 REAL xr = a[j1];
1882 REAL xi = -a[j1 + 1];
1883 REAL yr = a[k1];
1884 REAL yi = -a[k1 + 1];
1885 a[j1] = yr;
1886 a[j1 + 1] = yi;
1887 a[k1] = xr;
1888 a[k1 + 1] = xi;
1889 j1 += nm;
1890 k1 += nm;
1891 xr = a[j1];
1892 xi = -a[j1 + 1];
1893 yr = a[k1];
1894 yi = -a[k1 + 1];
1895 a[j1] = yr;
1896 a[j1 + 1] = yi;
1897 a[k1] = xr;
1898 a[k1 + 1] = xi;
1899 j1 += nh;
1900 k1 += 2;
1901 xr = a[j1];
1902 xi = -a[j1 + 1];
1903 yr = a[k1];
1904 yi = -a[k1 + 1];
1905 a[j1] = yr;
1906 a[j1 + 1] = yi;
1907 a[k1] = xr;
1908 a[k1 + 1] = xi;
1909 j1 -= nm;
1910 k1 -= nm;
1911 xr = a[j1];
1912 xi = -a[j1 + 1];
1913 yr = a[k1];
1914 yi = -a[k1 + 1];
1915 a[j1] = yr;
1916 a[j1 + 1] = yi;
1917 a[k1] = xr;
1918 a[k1 + 1] = xi;
1919 j1 += 2;
1920 k1 += nh;
1921 xr = a[j1];
1922 xi = -a[j1 + 1];
1923 yr = a[k1];
1924 yi = -a[k1 + 1];
1925 a[j1] = yr;
1926 a[j1 + 1] = yi;
1927 a[k1] = xr;
1928 a[k1 + 1] = xi;
1929 j1 += nm;
1930 k1 += nm;
1931 xr = a[j1];
1932 xi = -a[j1 + 1];
1933 yr = a[k1];
1934 yi = -a[k1 + 1];
1935 a[j1] = yr;
1936 a[j1 + 1] = yi;
1937 a[k1] = xr;
1938 a[k1 + 1] = xi;
1939 j1 -= nh;
1940 k1 -= 2;
1941 xr = a[j1];
1942 xi = -a[j1 + 1];
1943 yr = a[k1];
1944 yi = -a[k1 + 1];
1945 a[j1] = yr;
1946 a[j1 + 1] = yi;
1947 a[k1] = xr;
1948 a[k1 + 1] = xi;
1949 j1 -= nm;
1950 k1 -= nm;
1951 xr = a[j1];
1952 xi = -a[j1 + 1];
1953 yr = a[k1];
1954 yi = -a[k1 + 1];
1955 a[j1] = yr;
1956 a[j1 + 1] = yi;
1957 a[k1] = xr;
1958 a[k1 + 1] = xi;
1959 }
1960 k1 = 4 * k + ip[m + k];
1961 j1 = k1 + 2;
1962 k1 += nh;
1963 a[j1 - 1] = -a[j1 - 1];
1964 REAL xr = a[j1];
1965 REAL xi = -a[j1 + 1];
1966 REAL yr = a[k1];
1967 REAL yi = -a[k1 + 1];
1968 a[j1] = yr;
1969 a[j1 + 1] = yi;
1970 a[k1] = xr;
1971 a[k1 + 1] = xi;
1972 a[k1 + 3] = -a[k1 + 3];
1973 j1 += nm;
1974 k1 += nm;
1975 a[j1 - 1] = -a[j1 - 1];
1976 xr = a[j1];
1977 xi = -a[j1 + 1];
1978 yr = a[k1];
1979 yi = -a[k1 + 1];
1980 a[j1] = yr;
1981 a[j1 + 1] = yi;
1982 a[k1] = xr;
1983 a[k1 + 1] = xi;
1984 a[k1 + 3] = -a[k1 + 3];
1985 }
1986 }
1987 }
1988
1989
1990 void bitrv216(REAL *a)
1991 {
1992 /*
1993 REAL x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
1994 x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i,
1995 x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i;
1996 */
1997 immutable REAL x1r = a[2];
1998 immutable REAL x1i = a[3];
1999 immutable REAL x2r = a[4];
2000 immutable REAL x2i = a[5];
2001 immutable REAL x3r = a[6];
2002 immutable REAL x3i = a[7];
2003 immutable REAL x4r = a[8];
2004 immutable REAL x4i = a[9];
2005 immutable REAL x5r = a[10];
2006 immutable REAL x5i = a[11];
2007 immutable REAL x7r = a[14];
2008 immutable REAL x7i = a[15];
2009 immutable REAL x8r = a[16];
2010 immutable REAL x8i = a[17];
2011 immutable REAL x10r = a[20];
2012 immutable REAL x10i = a[21];
2013 immutable REAL x11r = a[22];
2014 immutable REAL x11i = a[23];
2015 immutable REAL x12r = a[24];
2016 immutable REAL x12i = a[25];
2017 immutable REAL x13r = a[26];
2018 immutable REAL x13i = a[27];
2019 immutable REAL x14r = a[28];
2020 immutable REAL x14i = a[29];
2021 a[2] = x8r;
2022 a[3] = x8i;
2023 a[4] = x4r;
2024 a[5] = x4i;
2025 a[6] = x12r;
2026 a[7] = x12i;
2027 a[8] = x2r;
2028 a[9] = x2i;
2029 a[10] = x10r;
2030 a[11] = x10i;
2031 a[14] = x14r;
2032 a[15] = x14i;
2033 a[16] = x1r;
2034 a[17] = x1i;
2035 a[20] = x5r;
2036 a[21] = x5i;
2037 a[22] = x13r;
2038 a[23] = x13i;
2039 a[24] = x3r;
2040 a[25] = x3i;
2041 a[26] = x11r;
2042 a[27] = x11i;
2043 a[28] = x7r;
2044 a[29] = x7i;
2045 }
2046
2047
2048 void bitrv216neg(REAL *a)
2049 {
2050 /*
2051 REAL x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
2052 x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i,
2053 x9r, x9i, x10r, x10i, x11r, x11i, x12r, x12i,
2054 x13r, x13i, x14r, x14i, x15r, x15i;
2055 */
2056 immutable REAL x1r = a[2];
2057 immutable REAL x1i = a[3];
2058 immutable REAL x2r = a[4];
2059 immutable REAL x2i = a[5];
2060 immutable REAL x3r = a[6];
2061 immutable REAL x3i = a[7];
2062 immutable REAL x4r = a[8];
2063 immutable REAL x4i = a[9];
2064 immutable REAL x5r = a[10];
2065 immutable REAL x5i = a[11];
2066 immutable REAL x6r = a[12];
2067 immutable REAL x6i = a[13];
2068 immutable REAL x7r = a[14];
2069 immutable REAL x7i = a[15];
2070 immutable REAL x8r = a[16];
2071 immutable REAL x8i = a[17];
2072 immutable REAL x9r = a[18];
2073 immutable REAL x9i = a[19];
2074 immutable REAL x10r = a[20];
2075 immutable REAL x10i = a[21];
2076 immutable REAL x11r = a[22];
2077 immutable REAL x11i = a[23];
2078 immutable REAL x12r = a[24];
2079 immutable REAL x12i = a[25];
2080 immutable REAL x13r = a[26];
2081 immutable REAL x13i = a[27];
2082 immutable REAL x14r = a[28];
2083 immutable REAL x14i = a[29];
2084 immutable REAL x15r = a[30];
2085 immutable REAL x15i = a[31];
2086 a[2] = x15r;
2087 a[3] = x15i;
2088 a[4] = x7r;
2089 a[5] = x7i;
2090 a[6] = x11r;
2091 a[7] = x11i;
2092 a[8] = x3r;
2093 a[9] = x3i;
2094 a[10] = x13r;
2095 a[11] = x13i;
2096 a[12] = x5r;
2097 a[13] = x5i;
2098 a[14] = x9r;
2099 a[15] = x9i;
2100 a[16] = x1r;
2101 a[17] = x1i;
2102 a[18] = x14r;
2103 a[19] = x14i;
2104 a[20] = x6r;
2105 a[21] = x6i;
2106 a[22] = x10r;
2107 a[23] = x10i;
2108 a[24] = x2r;
2109 a[25] = x2i;
2110 a[26] = x12r;
2111 a[27] = x12i;
2112 a[28] = x4r;
2113 a[29] = x4i;
2114 a[30] = x8r;
2115 a[31] = x8i;
2116 }
2117
2118
2119 void bitrv208(REAL *a)
2120 {
2121 //REAL x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i;
2122
2123 immutable REAL x1r = a[2];
2124 immutable REAL x1i = a[3];
2125 immutable REAL x3r = a[6];
2126 immutable REAL x3i = a[7];
2127 immutable REAL x4r = a[8];
2128 immutable REAL x4i = a[9];
2129 immutable REAL x6r = a[12];
2130 immutable REAL x6i = a[13];
2131 a[2] = x4r;
2132 a[3] = x4i;
2133 a[6] = x6r;
2134 a[7] = x6i;
2135 a[8] = x1r;
2136 a[9] = x1i;
2137 a[12] = x3r;
2138 a[13] = x3i;
2139 }
2140
2141
2142 void bitrv208neg(REAL *a)
2143 {
2144 //REAL x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, x5r, x5i, x6r, x6i, x7r, x7i;
2145
2146 immutable REAL x1r = a[2];
2147 immutable REAL x1i = a[3];
2148 immutable REAL x2r = a[4];
2149 immutable REAL x2i = a[5];
2150 immutable REAL x3r = a[6];
2151 immutable REAL x3i = a[7];
2152 immutable REAL x4r = a[8];
2153 immutable REAL x4i = a[9];
2154 immutable REAL x5r = a[10];
2155 immutable REAL x5i = a[11];
2156 immutable REAL x6r = a[12];
2157 immutable REAL x6i = a[13];
2158 immutable REAL x7r = a[14];
2159 immutable REAL x7i = a[15];
2160 a[2] = x7r;
2161 a[3] = x7i;
2162 a[4] = x3r;
2163 a[5] = x3i;
2164 a[6] = x5r;
2165 a[7] = x5i;
2166 a[8] = x1r;
2167 a[9] = x1i;
2168 a[10] = x6r;
2169 a[11] = x6i;
2170 a[12] = x2r;
2171 a[13] = x2i;
2172 a[14] = x4r;
2173 a[15] = x4i;
2174 }
2175
2176
2177 void cftf1st(int n, REAL *a, REAL *w)
2178 {
2179 int j, j0, j1, j2, j3, k, m, mh;
2180 //REAL wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i;
2181 //REAL x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
2182
2183 mh = n >> 3;
2184 m = 2 * mh;
2185 j1 = m;
2186 j2 = j1 + m;
2187 j3 = j2 + m;
2188 REAL x0r = a[0] + a[j2];
2189 REAL x0i = a[1] + a[j2 + 1];
2190 REAL x1r = a[0] - a[j2];
2191 REAL x1i = a[1] - a[j2 + 1];
2192 REAL x2r = a[j1] + a[j3];
2193 REAL x2i = a[j1 + 1] + a[j3 + 1];
2194 REAL x3r = a[j1] - a[j3];
2195 REAL x3i = a[j1 + 1] - a[j3 + 1];
2196 a[0] = x0r + x2r;
2197 a[1] = x0i + x2i;
2198 a[j1] = x0r - x2r;
2199 a[j1 + 1] = x0i - x2i;
2200 a[j2] = x1r - x3i;
2201 a[j2 + 1] = x1i + x3r;
2202 a[j3] = x1r + x3i;
2203 a[j3 + 1] = x1i - x3r;
2204 REAL wn4r = w[1];
2205 REAL csc1 = w[2];
2206 REAL csc3 = w[3];
2207 REAL wd1r = 1;
2208 REAL wd1i = 0;
2209 REAL wd3r = 1;
2210 REAL wd3i = 0;
2211 k = 0;
2212 for (j = 2; j < mh - 2; j += 4) {
2213 k += 4;
2214 REAL wk1r = csc1 * (wd1r + w[k]);
2215 REAL wk1i = csc1 * (wd1i + w[k + 1]);
2216 REAL wk3r = csc3 * (wd3r + w[k + 2]);
2217 REAL wk3i = csc3 * (wd3i + w[k + 3]);
2218 wd1r = w[k];
2219 wd1i = w[k + 1];
2220 wd3r = w[k + 2];
2221 wd3i = w[k + 3];
2222 j1 = j + m;
2223 j2 = j1 + m;
2224 j3 = j2 + m;
2225 x0r = a[j] + a[j2];
2226 x0i = a[j + 1] + a[j2 + 1];
2227 x1r = a[j] - a[j2];
2228 x1i = a[j + 1] - a[j2 + 1];
2229 REAL y0r = a[j + 2] + a[j2 + 2];
2230 REAL y0i = a[j + 3] + a[j2 + 3];
2231 REAL y1r = a[j + 2] - a[j2 + 2];
2232 REAL y1i = a[j + 3] - a[j2 + 3];
2233 x2r = a[j1] + a[j3];
2234 x2i = a[j1 + 1] + a[j3 + 1];
2235 x3r = a[j1] - a[j3];
2236 x3i = a[j1 + 1] - a[j3 + 1];
2237 REAL y2r = a[j1 + 2] + a[j3 + 2];
2238 REAL y2i = a[j1 + 3] + a[j3 + 3];
2239 REAL y3r = a[j1 + 2] - a[j3 + 2];
2240 REAL y3i = a[j1 + 3] - a[j3 + 3];
2241 a[j] = x0r + x2r;
2242 a[j + 1] = x0i + x2i;
2243 a[j + 2] = y0r + y2r;
2244 a[j + 3] = y0i + y2i;
2245 a[j1] = x0r - x2r;
2246 a[j1 + 1] = x0i - x2i;
2247 a[j1 + 2] = y0r - y2r;
2248 a[j1 + 3] = y0i - y2i;
2249 x0r = x1r - x3i;
2250 x0i = x1i + x3r;
2251 a[j2] = wk1r * x0r - wk1i * x0i;
2252 a[j2 + 1] = wk1r * x0i + wk1i * x0r;
2253 x0r = y1r - y3i;
2254 x0i = y1i + y3r;
2255 a[j2 + 2] = wd1r * x0r - wd1i * x0i;
2256 a[j2 + 3] = wd1r * x0i + wd1i * x0r;
2257 x0r = x1r + x3i;
2258 x0i = x1i - x3r;
2259 a[j3] = wk3r * x0r + wk3i * x0i;
2260 a[j3 + 1] = wk3r * x0i - wk3i * x0r;
2261 x0r = y1r + y3i;
2262 x0i = y1i - y3r;
2263 a[j3 + 2] = wd3r * x0r + wd3i * x0i;
2264 a[j3 + 3] = wd3r * x0i - wd3i * x0r;
2265 j0 = m - j;
2266 j1 = j0 + m;
2267 j2 = j1 + m;
2268 j3 = j2 + m;
2269 x0r = a[j0] + a[j2];
2270 x0i = a[j0 + 1] + a[j2 + 1];
2271 x1r = a[j0] - a[j2];
2272 x1i = a[j0 + 1] - a[j2 + 1];
2273 y0r = a[j0 - 2] + a[j2 - 2];
2274 y0i = a[j0 - 1] + a[j2 - 1];
2275 y1r = a[j0 - 2] - a[j2 - 2];
2276 y1i = a[j0 - 1] - a[j2 - 1];
2277 x2r = a[j1] + a[j3];
2278 x2i = a[j1 + 1] + a[j3 + 1];
2279 x3r = a[j1] - a[j3];
2280 x3i = a[j1 + 1] - a[j3 + 1];
2281 y2r = a[j1 - 2] + a[j3 - 2];
2282 y2i = a[j1 - 1] + a[j3 - 1];
2283 y3r = a[j1 - 2] - a[j3 - 2];
2284 y3i = a[j1 - 1] - a[j3 - 1];
2285 a[j0] = x0r + x2r;
2286 a[j0 + 1] = x0i + x2i;
2287 a[j0 - 2] = y0r + y2r;
2288 a[j0 - 1] = y0i + y2i;
2289 a[j1] = x0r - x2r;
2290 a[j1 + 1] = x0i - x2i;
2291 a[j1 - 2] = y0r - y2r;
2292 a[j1 - 1] = y0i - y2i;
2293 x0r = x1r - x3i;
2294 x0i = x1i + x3r;
2295 a[j2] = wk1i * x0r - wk1r * x0i;
2296 a[j2 + 1] = wk1i * x0i + wk1r * x0r;
2297 x0r = y1r - y3i;
2298 x0i = y1i + y3r;
2299 a[j2 - 2] = wd1i * x0r - wd1r * x0i;
2300 a[j2 - 1] = wd1i * x0i + wd1r * x0r;
2301 x0r = x1r + x3i;
2302 x0i = x1i - x3r;
2303 a[j3] = wk3i * x0r + wk3r * x0i;
2304 a[j3 + 1] = wk3i * x0i - wk3r * x0r;
2305 x0r = y1r + y3i;
2306 x0i = y1i - y3r;
2307 a[j3 - 2] = wd3i * x0r + wd3r * x0i;
2308 a[j3 - 1] = wd3i * x0i - wd3r * x0r;
2309 }
2310 REAL wk1r = csc1 * (wd1r + wn4r);
2311 REAL wk1i = csc1 * (wd1i + wn4r);
2312 REAL wk3r = csc3 * (wd3r - wn4r);
2313 REAL wk3i = csc3 * (wd3i - wn4r);
2314 j0 = mh;
2315 j1 = j0 + m;
2316 j2 = j1 + m;
2317 j3 = j2 + m;
2318 x0r = a[j0 - 2] + a[j2 - 2];
2319 x0i = a[j0 - 1] + a[j2 - 1];
2320 x1r = a[j0 - 2] - a[j2 - 2];
2321 x1i = a[j0 - 1] - a[j2 - 1];
2322 x2r = a[j1 - 2] + a[j3 - 2];
2323 x2i = a[j1 - 1] + a[j3 - 1];
2324 x3r = a[j1 - 2] - a[j3 - 2];
2325 x3i = a[j1 - 1] - a[j3 - 1];
2326 a[j0 - 2] = x0r + x2r;
2327 a[j0 - 1] = x0i + x2i;
2328 a[j1 - 2] = x0r - x2r;
2329 a[j1 - 1] = x0i - x2i;
2330 x0r = x1r - x3i;
2331 x0i = x1i + x3r;
2332 a[j2 - 2] = wk1r * x0r - wk1i * x0i;
2333 a[j2 - 1] = wk1r * x0i + wk1i * x0r;
2334 x0r = x1r + x3i;
2335 x0i = x1i - x3r;
2336 a[j3 - 2] = wk3r * x0r + wk3i * x0i;
2337 a[j3 - 1] = wk3r * x0i - wk3i * x0r;
2338 x0r = a[j0] + a[j2];
2339 x0i = a[j0 + 1] + a[j2 + 1];
2340 x1r = a[j0] - a[j2];
2341 x1i = a[j0 + 1] - a[j2 + 1];
2342 x2r = a[j1] + a[j3];
2343 x2i = a[j1 + 1] + a[j3 + 1];
2344 x3r = a[j1] - a[j3];
2345 x3i = a[j1 + 1] - a[j3 + 1];
2346 a[j0] = x0r + x2r;
2347 a[j0 + 1] = x0i + x2i;
2348 a[j1] = x0r - x2r;
2349 a[j1 + 1] = x0i - x2i;
2350 x0r = x1r - x3i;
2351 x0i = x1i + x3r;
2352 a[j2] = wn4r * (x0r - x0i);
2353 a[j2 + 1] = wn4r * (x0i + x0r);
2354 x0r = x1r + x3i;
2355 x0i = x1i - x3r;
2356 a[j3] = -wn4r * (x0r + x0i);
2357 a[j3 + 1] = -wn4r * (x0i - x0r);
2358 x0r = a[j0 + 2] + a[j2 + 2];
2359 x0i = a[j0 + 3] + a[j2 + 3];
2360 x1r = a[j0 + 2] - a[j2 + 2];
2361 x1i = a[j0 + 3] - a[j2 + 3];
2362 x2r = a[j1 + 2] + a[j3 + 2];
2363 x2i = a[j1 + 3] + a[j3 + 3];
2364 x3r = a[j1 + 2] - a[j3 + 2];
2365 x3i = a[j1 + 3] - a[j3 + 3];
2366 a[j0 + 2] = x0r + x2r;
2367 a[j0 + 3] = x0i + x2i;
2368 a[j1 + 2] = x0r - x2r;
2369 a[j1 + 3] = x0i - x2i;
2370 x0r = x1r - x3i;
2371 x0i = x1i + x3r;
2372 a[j2 + 2] = wk1i * x0r - wk1r * x0i;
2373 a[j2 + 3] = wk1i * x0i + wk1r * x0r;
2374 x0r = x1r + x3i;
2375 x0i = x1i - x3r;
2376 a[j3 + 2] = wk3i * x0r + wk3r * x0i;
2377 a[j3 + 3] = wk3i * x0i - wk3r * x0r;
2378 }
2379
2380
2381 void cftb1st(int n, REAL *a, REAL *w)
2382 {
2383 int j, j0, j1, j2, j3, k, m, mh;
2384 //REAL wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i;
2385 //REAL x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
2386
2387 mh = n >> 3;
2388 m = 2 * mh;
2389 j1 = m;
2390 j2 = j1 + m;
2391 j3 = j2 + m;
2392 REAL x0r = a[0] + a[j2];
2393 REAL x0i = -a[1] - a[j2 + 1];
2394 REAL x1r = a[0] - a[j2];
2395 REAL x1i = -a[1] + a[j2 + 1];
2396 REAL x2r = a[j1] + a[j3];
2397 REAL x2i = a[j1 + 1] + a[j3 + 1];
2398 REAL x3r = a[j1] - a[j3];
2399 REAL x3i = a[j1 + 1] - a[j3 + 1];
2400 a[0] = x0r + x2r;
2401 a[1] = x0i - x2i;
2402 a[j1] = x0r - x2r;
2403 a[j1 + 1] = x0i + x2i;
2404 a[j2] = x1r + x3i;
2405 a[j2 + 1] = x1i + x3r;
2406 a[j3] = x1r - x3i;
2407 a[j3 + 1] = x1i - x3r;
2408 REAL wn4r = w[1];
2409 REAL csc1 = w[2];
2410 REAL csc3 = w[3];
2411 REAL wd1r = 1;
2412 REAL wd1i = 0;
2413 REAL wd3r = 1;
2414 REAL wd3i = 0;
2415 k = 0;
2416 for (j = 2; j < mh - 2; j += 4) {
2417 k += 4;
2418 REAL wk1r = csc1 * (wd1r + w[k]);
2419 REAL wk1i = csc1 * (wd1i + w[k + 1]);
2420 REAL wk3r = csc3 * (wd3r + w[k + 2]);
2421 REAL wk3i = csc3 * (wd3i + w[k + 3]);
2422 wd1r = w[k];
2423 wd1i = w[k + 1];
2424 wd3r = w[k + 2];
2425 wd3i = w[k + 3];
2426 j1 = j + m;
2427 j2 = j1 + m;
2428 j3 = j2 + m;
2429 x0r = a[j] + a[j2];
2430 x0i = -a[j + 1] - a[j2 + 1];
2431 x1r = a[j] - a[j2];
2432 x1i = -a[j + 1] + a[j2 + 1];
2433 REAL y0r = a[j + 2] + a[j2 + 2];
2434 REAL y0i = -a[j + 3] - a[j2 + 3];
2435 REAL y1r = a[j + 2] - a[j2 + 2];
2436 REAL y1i = -a[j + 3] + a[j2 + 3];
2437 x2r = a[j1] + a[j3];
2438 x2i = a[j1 + 1] + a[j3 + 1];
2439 x3r = a[j1] - a[j3];
2440 x3i = a[j1 + 1] - a[j3 + 1];
2441 REAL y2r = a[j1 + 2] + a[j3 + 2];
2442 REAL y2i = a[j1 + 3] + a[j3 + 3];
2443 REAL y3r = a[j1 + 2] - a[j3 + 2];
2444 REAL y3i = a[j1 + 3] - a[j3 + 3];
2445 a[j] = x0r + x2r;
2446 a[j + 1] = x0i - x2i;
2447 a[j + 2] = y0r + y2r;
2448 a[j + 3] = y0i - y2i;
2449 a[j1] = x0r - x2r;
2450 a[j1 + 1] = x0i + x2i;
2451 a[j1 + 2] = y0r - y2r;
2452 a[j1 + 3] = y0i + y2i;
2453 x0r = x1r + x3i;
2454 x0i = x1i + x3r;
2455 a[j2] = wk1r * x0r - wk1i * x0i;
2456 a[j2 + 1] = wk1r * x0i + wk1i * x0r;
2457 x0r = y1r + y3i;
2458 x0i = y1i + y3r;
2459 a[j2 + 2] = wd1r * x0r - wd1i * x0i;
2460 a[j2 + 3] = wd1r * x0i + wd1i * x0r;
2461 x0r = x1r - x3i;
2462 x0i = x1i - x3r;
2463 a[j3] = wk3r * x0r + wk3i * x0i;
2464 a[j3 + 1] = wk3r * x0i - wk3i * x0r;
2465 x0r = y1r - y3i;
2466 x0i = y1i - y3r;
2467 a[j3 + 2] = wd3r * x0r + wd3i * x0i;
2468 a[j3 + 3] = wd3r * x0i - wd3i * x0r;
2469 j0 = m - j;
2470 j1 = j0 + m;
2471 j2 = j1 + m;
2472 j3 = j2 + m;
2473 x0r = a[j0] + a[j2];
2474 x0i = -a[j0 + 1] - a[j2 + 1];
2475 x1r = a[j0] - a[j2];
2476 x1i = -a[j0 + 1] + a[j2 + 1];
2477 y0r = a[j0 - 2] + a[j2 - 2];
2478 y0i = -a[j0 - 1] - a[j2 - 1];
2479 y1r = a[j0 - 2] - a[j2 - 2];
2480 y1i = -a[j0 - 1] + a[j2 - 1];
2481 x2r = a[j1] + a[j3];
2482 x2i = a[j1 + 1] + a[j3 + 1];
2483 x3r = a[j1] - a[j3];
2484 x3i = a[j1 + 1] - a[j3 + 1];
2485 y2r = a[j1 - 2] + a[j3 - 2];
2486 y2i = a[j1 - 1] + a[j3 - 1];
2487 y3r = a[j1 - 2] - a[j3 - 2];
2488 y3i = a[j1 - 1] - a[j3 - 1];
2489 a[j0] = x0r + x2r;
2490 a[j0 + 1] = x0i - x2i;
2491 a[j0 - 2] = y0r + y2r;
2492 a[j0 - 1] = y0i - y2i;
2493 a[j1] = x0r - x2r;
2494 a[j1 + 1] = x0i + x2i;
2495 a[j1 - 2] = y0r - y2r;
2496 a[j1 - 1] = y0i + y2i;
2497 x0r = x1r + x3i;
2498 x0i = x1i + x3r;
2499 a[j2] = wk1i * x0r - wk1r * x0i;
2500 a[j2 + 1] = wk1i * x0i + wk1r * x0r;
2501 x0r = y1r + y3i;
2502 x0i = y1i + y3r;
2503 a[j2 - 2] = wd1i * x0r - wd1r * x0i;
2504 a[j2 - 1] = wd1i * x0i + wd1r * x0r;
2505 x0r = x1r - x3i;
2506 x0i = x1i - x3r;
2507 a[j3] = wk3i * x0r + wk3r * x0i;
2508 a[j3 + 1] = wk3i * x0i - wk3r * x0r;
2509 x0r = y1r - y3i;
2510 x0i = y1i - y3r;
2511 a[j3 - 2] = wd3i * x0r + wd3r * x0i;
2512 a[j3 - 1] = wd3i * x0i - wd3r * x0r;
2513 }
2514 REAL wk1r = csc1 * (wd1r + wn4r);
2515 REAL wk1i = csc1 * (wd1i + wn4r);
2516 REAL wk3r = csc3 * (wd3r - wn4r);
2517 REAL wk3i = csc3 * (wd3i - wn4r);
2518 j0 = mh;
2519 j1 = j0 + m;
2520 j2 = j1 + m;
2521 j3 = j2 + m;
2522 x0r = a[j0 - 2] + a[j2 - 2];
2523 x0i = -a[j0 - 1] - a[j2 - 1];
2524 x1r = a[j0 - 2] - a[j2 - 2];
2525 x1i = -a[j0 - 1] + a[j2 - 1];
2526 x2r = a[j1 - 2] + a[j3 - 2];
2527 x2i = a[j1 - 1] + a[j3 - 1];
2528 x3r = a[j1 - 2] - a[j3 - 2];
2529 x3i = a[j1 - 1] - a[j3 - 1];
2530 a[j0 - 2] = x0r + x2r;
2531 a[j0 - 1] = x0i - x2i;
2532 a[j1 - 2] = x0r - x2r;
2533 a[j1 - 1] = x0i + x2i;
2534 x0r = x1r + x3i;
2535 x0i = x1i + x3r;
2536 a[j2 - 2] = wk1r * x0r - wk1i * x0i;
2537 a[j2 - 1] = wk1r * x0i + wk1i * x0r;
2538 x0r = x1r - x3i;
2539 x0i = x1i - x3r;
2540 a[j3 - 2] = wk3r * x0r + wk3i * x0i;
2541 a[j3 - 1] = wk3r * x0i - wk3i * x0r;
2542 x0r = a[j0] + a[j2];
2543 x0i = -a[j0 + 1] - a[j2 + 1];
2544 x1r = a[j0] - a[j2];
2545 x1i = -a[j0 + 1] + a[j2 + 1];
2546 x2r = a[j1] + a[j3];
2547 x2i = a[j1 + 1] + a[j3 + 1];
2548 x3r = a[j1] - a[j3];
2549 x3i = a[j1 + 1] - a[j3 + 1];
2550 a[j0] = x0r + x2r;
2551 a[j0 + 1] = x0i - x2i;
2552 a[j1] = x0r - x2r;
2553 a[j1 + 1] = x0i + x2i;
2554 x0r = x1r + x3i;
2555 x0i = x1i + x3r;
2556 a[j2] = wn4r * (x0r - x0i);
2557 a[j2 + 1] = wn4r * (x0i + x0r);
2558 x0r = x1r - x3i;
2559 x0i = x1i - x3r;
2560 a[j3] = -wn4r * (x0r + x0i);
2561 a[j3 + 1] = -wn4r * (x0i - x0r);
2562 x0r = a[j0 + 2] + a[j2 + 2];
2563 x0i = -a[j0 + 3] - a[j2 + 3];
2564 x1r = a[j0 + 2] - a[j2 + 2];
2565 x1i = -a[j0 + 3] + a[j2 + 3];
2566 x2r = a[j1 + 2] + a[j3 + 2];
2567 x2i = a[j1 + 3] + a[j3 + 3];
2568 x3r = a[j1 + 2] - a[j3 + 2];
2569 x3i = a[j1 + 3] - a[j3 + 3];
2570 a[j0 + 2] = x0r + x2r;
2571 a[j0 + 3] = x0i - x2i;
2572 a[j1 + 2] = x0r - x2r;
2573 a[j1 + 3] = x0i + x2i;
2574 x0r = x1r + x3i;
2575 x0i = x1i + x3r;
2576 a[j2 + 2] = wk1i * x0r - wk1r * x0i;
2577 a[j2 + 3] = wk1i * x0i + wk1r * x0r;
2578 x0r = x1r - x3i;
2579 x0i = x1i - x3r;
2580 a[j3 + 2] = wk3i * x0r + wk3r * x0i;
2581 a[j3 + 3] = wk3i * x0i - wk3r * x0r;
2582 }
2583
2584
2585 void cftrec4(int n, REAL *a, int nw, REAL *w)
2586 {
2587 int isplt, j, k, m;
2588
2589 m = n;
2590 while (m > 512) {
2591 m >>= 2;
2592 cftmdl1(m, &a[n - m], &w[nw - (m >> 1)]);
2593 }
2594 cftleaf(m, 1, &a[n - m], nw, w);
2595 k = 0;
2596 for (j = n - m; j > 0; j -= m) {
2597 k++;
2598 isplt = cfttree(m, j, k, a, nw, w);
2599 cftleaf(m, isplt, &a[j - m], nw, w);
2600 }
2601 }
2602
2603
2604 int cfttree(int n, int j, int k, REAL *a, int nw, REAL *w)
2605 {
2606 int i, isplt, m;
2607
2608 if ((k & 3) != 0) {
2609 isplt = k & 1;
2610 if (isplt != 0) {
2611 cftmdl1(n, &a[j - n], &w[nw - (n >> 1)]);
2612 } else {
2613 cftmdl2(n, &a[j - n], &w[nw - n]);
2614 }
2615 } else {
2616 m = n;
2617 for (i = k; (i & 3) == 0; i >>= 2) {
2618 m <<= 2;
2619 }
2620 isplt = i & 1;
2621 if (isplt != 0) {
2622 while (m > 128) {
2623 cftmdl1(m, &a[j - m], &w[nw - (m >> 1)]);
2624 m >>= 2;
2625 }
2626 } else {
2627 while (m > 128) {
2628 cftmdl2(m, &a[j - m], &w[nw - m]);
2629 m >>= 2;
2630 }
2631 }
2632 }
2633 return isplt;
2634 }
2635
2636
2637 void cftleaf(int n, int isplt, REAL *a, int nw, REAL *w)
2638 {
2639 if (n == 512) {
2640 cftmdl1(128, a, &w[nw - 64]);
2641 cftf161(a, &w[nw - 8]);
2642 cftf162(&a[32], &w[nw - 32]);
2643 cftf161(&a[64], &w[nw - 8]);
2644 cftf161(&a[96], &w[nw - 8]);
2645 cftmdl2(128, &a[128], &w[nw - 128]);
2646 cftf161(&a[128], &w[nw - 8]);
2647 cftf162(&a[160], &w[nw - 32]);
2648 cftf161(&a[192], &w[nw - 8]);
2649 cftf162(&a[224], &w[nw - 32]);
2650 cftmdl1(128, &a[256], &w[nw - 64]);
2651 cftf161(&a[256], &w[nw - 8]);
2652 cftf162(&a[288], &w[nw - 32]);
2653 cftf161(&a[320], &w[nw - 8]);
2654 cftf161(&a[352], &w[nw - 8]);
2655 if (isplt != 0) {
2656 cftmdl1(128, &a[384], &w[nw - 64]);
2657 cftf161(&a[480], &w[nw - 8]);
2658 } else {
2659 cftmdl2(128, &a[384], &w[nw - 128]);
2660 cftf162(&a[480], &w[nw - 32]);
2661 }
2662 cftf161(&a[384], &w[nw - 8]);
2663 cftf162(&a[416], &w[nw - 32]);
2664 cftf161(&a[448], &w[nw - 8]);
2665 } else {
2666 cftmdl1(64, a, &w[nw - 32]);
2667 cftf081(a, &w[nw - 8]);
2668 cftf082(&a[16], &w[nw - 8]);
2669 cftf081(&a[32], &w[nw - 8]);
2670 cftf081(&a[48], &w[nw - 8]);
2671 cftmdl2(64, &a[64], &w[nw - 64]);
2672 cftf081(&a[64], &w[nw - 8]);
2673 cftf082(&a[80], &w[nw - 8]);
2674 cftf081(&a[96], &w[nw - 8]);
2675 cftf082(&a[112], &w[nw - 8]);
2676 cftmdl1(64, &a[128], &w[nw - 32]);
2677 cftf081(&a[128], &w[nw - 8]);
2678 cftf082(&a[144], &w[nw - 8]);
2679 cftf081(&a[160], &w[nw - 8]);
2680 cftf081(&a[176], &w[nw - 8]);
2681 if (isplt != 0) {
2682 cftmdl1(64, &a[192], &w[nw - 32]);
2683 cftf081(&a[240], &w[nw - 8]);
2684 } else {
2685 cftmdl2(64, &a[192], &w[nw - 64]);
2686 cftf082(&a[240], &w[nw - 8]);
2687 }
2688 cftf081(&a[192], &w[nw - 8]);
2689 cftf082(&a[208], &w[nw - 8]);
2690 cftf081(&a[224], &w[nw - 8]);
2691 }
2692 }
2693
2694
2695 void cftmdl1(int n, REAL *a, REAL *w)
2696 {
2697 int j, j0, j1, j2, j3, k, m, mh;
2698 //REAL wn4r, wk1r, wk1i, wk3r, wk3i;
2699 //REAL x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
2700
2701 mh = n >> 3;
2702 m = 2 * mh;
2703 j1 = m;
2704 j2 = j1 + m;
2705 j3 = j2 + m;
2706 REAL x0r = a[0] + a[j2];
2707 REAL x0i = a[1] + a[j2 + 1];
2708 REAL x1r = a[0] - a[j2];
2709 REAL x1i = a[1] - a[j2 + 1];
2710 REAL x2r = a[j1] + a[j3];
2711 REAL x2i = a[j1 + 1] + a[j3 + 1];
2712 REAL x3r = a[j1] - a[j3];
2713 REAL x3i = a[j1 + 1] - a[j3 + 1];
2714 a[0] = x0r + x2r;
2715 a[1] = x0i + x2i;
2716 a[j1] = x0r - x2r;
2717 a[j1 + 1] = x0i - x2i;
2718 a[j2] = x1r - x3i;
2719 a[j2 + 1] = x1i + x3r;
2720 a[j3] = x1r + x3i;
2721 a[j3 + 1] = x1i - x3r;
2722 REAL wn4r = w[1];
2723 k = 0;
2724 for (j = 2; j < mh; j += 2) {
2725 k += 4;
2726 REAL wk1r = w[k];
2727 REAL wk1i = w[k + 1];
2728 REAL wk3r = w[k + 2];
2729 REAL wk3i = w[k + 3];
2730 j1 = j + m;
2731 j2 = j1 + m;
2732 j3 = j2 + m;
2733 x0r = a[j] + a[j2];
2734 x0i = a[j + 1] + a[j2 + 1];
2735 x1r = a[j] - a[j2];
2736 x1i = a[j + 1] - a[j2 + 1];
2737 x2r = a[j1] + a[j3];
2738 x2i = a[j1 + 1] + a[j3 + 1];
2739 x3r = a[j1] - a[j3];
2740 x3i = a[j1 + 1] - a[j3 + 1];
2741 a[j] = x0r + x2r;
2742 a[j + 1] = x0i + x2i;
2743 a[j1] = x0r - x2r;
2744 a[j1 + 1] = x0i - x2i;
2745 x0r = x1r - x3i;
2746 x0i = x1i + x3r;
2747 a[j2] = wk1r * x0r - wk1i * x0i;
2748 a[j2 + 1] = wk1r * x0i + wk1i * x0r;
2749 x0r = x1r + x3i;
2750 x0i = x1i - x3r;
2751 a[j3] = wk3r * x0r + wk3i * x0i;
2752 a[j3 + 1] = wk3r * x0i - wk3i * x0r;
2753 j0 = m - j;
2754 j1 = j0 + m;
2755 j2 = j1 + m;
2756 j3 = j2 + m;
2757 x0r = a[j0] + a[j2];
2758 x0i = a[j0 + 1] + a[j2 + 1];
2759 x1r = a[j0] - a[j2];
2760 x1i = a[j0 + 1] - a[j2 + 1];
2761 x2r = a[j1] + a[j3];
2762 x2i = a[j1 + 1] + a[j3 + 1];
2763 x3r = a[j1] - a[j3];
2764 x3i = a[j1 + 1] - a[j3 + 1];
2765 a[j0] = x0r + x2r;
2766 a[j0 + 1] = x0i + x2i;
2767 a[j1] = x0r - x2r;
2768 a[j1 + 1] = x0i - x2i;
2769 x0r = x1r - x3i;
2770 x0i = x1i + x3r;
2771 a[j2] = wk1i * x0r - wk1r * x0i;
2772 a[j2 + 1] = wk1i * x0i + wk1r * x0r;
2773 x0r = x1r + x3i;
2774 x0i = x1i - x3r;
2775 a[j3] = wk3i * x0r + wk3r * x0i;
2776 a[j3 + 1] = wk3i * x0i - wk3r * x0r;
2777 }
2778 j0 = mh;
2779 j1 = j0 + m;
2780 j2 = j1 + m;
2781 j3 = j2 + m;
2782 x0r = a[j0] + a[j2];
2783 x0i = a[j0 + 1] + a[j2 + 1];
2784 x1r = a[j0] - a[j2];
2785 x1i = a[j0 + 1] - a[j2 + 1];
2786 x2r = a[j1] + a[j3];
2787 x2i = a[j1 + 1] + a[j3 + 1];
2788 x3r = a[j1] - a[j3];
2789 x3i = a[j1 + 1] - a[j3 + 1];
2790 a[j0] = x0r + x2r;
2791 a[j0 + 1] = x0i + x2i;
2792 a[j1] = x0r - x2r;
2793 a[j1 + 1] = x0i - x2i;
2794 x0r = x1r - x3i;
2795 x0i = x1i + x3r;
2796 a[j2] = wn4r * (x0r - x0i);
2797 a[j2 + 1] = wn4r * (x0i + x0r);
2798 x0r = x1r + x3i;
2799 x0i = x1i - x3r;
2800 a[j3] = -wn4r * (x0r + x0i);
2801 a[j3 + 1] = -wn4r * (x0i - x0r);
2802 }
2803
2804
2805 void cftmdl2(int n, REAL *a, REAL *w)
2806 {
2807 int j, j0, j1, j2, j3, k, kr, m, mh;
2808 //REAL wn4r, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i;
2809 //REAL x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y2r, y2i;
2810
2811 mh = n >> 3;
2812 m = 2 * mh;
2813 REAL wn4r = w[1];
2814 j1 = m;
2815 j2 = j1 + m;
2816 j3 = j2 + m;
2817 REAL x0r = a[0] - a[j2 + 1];
2818 REAL x0i = a[1] + a[j2];
2819 REAL x1r = a[0] + a[j2 + 1];
2820 REAL x1i = a[1] - a[j2];
2821 REAL x2r = a[j1] - a[j3 + 1];
2822 REAL x2i = a[j1 + 1] + a[j3];
2823 REAL x3r = a[j1] + a[j3 + 1];
2824 REAL x3i = a[j1 + 1] - a[j3];
2825 REAL y0r = wn4r * (x2r - x2i);
2826 REAL y0i = wn4r * (x2i + x2r);
2827 a[0] = x0r + y0r;
2828 a[1] = x0i + y0i;
2829 a[j1] = x0r - y0r;
2830 a[j1 + 1] = x0i - y0i;
2831 y0r = wn4r * (x3r - x3i);
2832 y0i = wn4r * (x3i + x3r);
2833 a[j2] = x1r - y0i;
2834 a[j2 + 1] = x1i + y0r;
2835 a[j3] = x1r + y0i;
2836 a[j3 + 1] = x1i - y0r;
2837 k = 0;
2838 kr = 2 * m;
2839 for (j = 2; j < mh; j += 2) {
2840 k += 4;
2841 REAL wk1r = w[k];
2842 REAL wk1i = w[k + 1];
2843 REAL wk3r = w[k + 2];
2844 REAL wk3i = w[k + 3];
2845 kr -= 4;
2846 REAL wd1i = w[kr];
2847 REAL wd1r = w[kr + 1];
2848 REAL wd3i = w[kr + 2];
2849 REAL wd3r = w[kr + 3];
2850 j1 = j + m;
2851 j2 = j1 + m;
2852 j3 = j2 + m;
2853 x0r = a[j] - a[j2 + 1];
2854 x0i = a[j + 1] + a[j2];
2855 x1r = a[j] + a[j2 + 1];
2856 x1i = a[j + 1] - a[j2];
2857 x2r = a[j1] - a[j3 + 1];
2858 x2i = a[j1 + 1] + a[j3];
2859 x3r = a[j1] + a[j3 + 1];
2860 x3i = a[j1 + 1] - a[j3];
2861 y0r = wk1r * x0r - wk1i * x0i;
2862 y0i = wk1r * x0i + wk1i * x0r;
2863 REAL y2r = wd1r * x2r - wd1i * x2i;
2864 REAL y2i = wd1r * x2i + wd1i * x2r;
2865 a[j] = y0r + y2r;
2866 a[j + 1] = y0i + y2i;
2867 a[j1] = y0r - y2r;
2868 a[j1 + 1] = y0i - y2i;
2869 y0r = wk3r * x1r + wk3i * x1i;
2870 y0i = wk3r * x1i - wk3i * x1r;
2871 y2r = wd3r * x3r + wd3i * x3i;
2872 y2i = wd3r * x3i - wd3i * x3r;
2873 a[j2] = y0r + y2r;
2874 a[j2 + 1] = y0i + y2i;
2875 a[j3] = y0r - y2r;
2876 a[j3 + 1] = y0i - y2i;
2877 j0 = m - j;
2878 j1 = j0 + m;
2879 j2 = j1 + m;
2880 j3 = j2 + m;
2881 x0r = a[j0] - a[j2 + 1];
2882 x0i = a[j0 + 1] + a[j2];
2883 x1r = a[j0] + a[j2 + 1];
2884 x1i = a[j0 + 1] - a[j2];
2885 x2r = a[j1] - a[j3 + 1];
2886 x2i = a[j1 + 1] + a[j3];
2887 x3r = a[j1] + a[j3 + 1];
2888 x3i = a[j1 + 1] - a[j3];
2889 y0r = wd1i * x0r - wd1r * x0i;
2890 y0i = wd1i * x0i + wd1r * x0r;
2891 y2r = wk1i * x2r - wk1r * x2i;
2892 y2i = wk1i * x2i + wk1r * x2r;
2893 a[j0] = y0r + y2r;
2894 a[j0 + 1] = y0i + y2i;
2895 a[j1] = y0r - y2r;
2896 a[j1 + 1] = y0i - y2i;
2897 y0r = wd3i * x1r + wd3r * x1i;
2898 y0i = wd3i * x1i - wd3r * x1r;
2899 y2r = wk3i * x3r + wk3r * x3i;
2900 y2i = wk3i * x3i - wk3r * x3r;
2901 a[j2] = y0r + y2r;
2902 a[j2 + 1] = y0i + y2i;
2903 a[j3] = y0r - y2r;
2904 a[j3 + 1] = y0i - y2i;
2905 }
2906 REAL wk1r = w[m];
2907 REAL wk1i = w[m + 1];
2908 j0 = mh;
2909 j1 = j0 + m;
2910 j2 = j1 + m;
2911 j3 = j2 + m;
2912 x0r = a[j0] - a[j2 + 1];
2913 x0i = a[j0 + 1] + a[j2];
2914 x1r = a[j0] + a[j2 + 1];
2915 x1i = a[j0 + 1] - a[j2];
2916 x2r = a[j1] - a[j3 + 1];
2917 x2i = a[j1 + 1] + a[j3];
2918 x3r = a[j1] + a[j3 + 1];
2919 x3i = a[j1 + 1] - a[j3];
2920 y0r = wk1r * x0r - wk1i * x0i;
2921 y0i = wk1r * x0i + wk1i * x0r;
2922 REAL y2r = wk1i * x2r - wk1r * x2i;
2923 REAL y2i = wk1i * x2i + wk1r * x2r;
2924 a[j0] = y0r + y2r;
2925 a[j0 + 1] = y0i + y2i;
2926 a[j1] = y0r - y2r;
2927 a[j1 + 1] = y0i - y2i;
2928 y0r = wk1i * x1r - wk1r * x1i;
2929 y0i = wk1i * x1i + wk1r * x1r;
2930 y2r = wk1r * x3r - wk1i * x3i;
2931 y2i = wk1r * x3i + wk1i * x3r;
2932 a[j2] = y0r - y2r;
2933 a[j2 + 1] = y0i - y2i;
2934 a[j3] = y0r + y2r;
2935 a[j3 + 1] = y0i + y2i;
2936 }
2937
2938
2939 void cftfx41(int n, REAL *a, int nw, REAL *w)
2940 {
2941 if (n == 128) {
2942 cftf161(a, &w[nw - 8]);
2943 cftf162(&a[32], &w[nw - 32]);
2944 cftf161(&a[64], &w[nw - 8]);
2945 cftf161(&a[96], &w[nw - 8]);
2946 } else {
2947 cftf081(a, &w[nw - 8]);
2948 cftf082(&a[16], &w[nw - 8]);
2949 cftf081(&a[32], &w[nw - 8]);
2950 cftf081(&a[48], &w[nw - 8]);
2951 }
2952 }
2953
2954
2955 void cftf161(REAL *a, REAL *w)
2956 {
2957 /*
2958 REAL wn4r, wk1r, wk1i,
2959 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
2960 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
2961 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i,
2962 y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i,
2963 y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
2964 */
2965 REAL wn4r = w[1];
2966 REAL wk1r = w[2];
2967 REAL wk1i = w[3];
2968 REAL x0r = a[0] + a[16];
2969 REAL x0i = a[1] + a[17];
2970 REAL x1r = a[0] - a[16];
2971 REAL x1i = a[1] - a[17];
2972 REAL x2r = a[8] + a[24];
2973 REAL x2i = a[9] + a[25];
2974 REAL x3r = a[8] - a[24];
2975 REAL x3i = a[9] - a[25];
2976 REAL y0r = x0r + x2r;
2977 REAL y0i = x0i + x2i;
2978 REAL y4r = x0r - x2r;
2979 REAL y4i = x0i - x2i;
2980 REAL y8r = x1r - x3i;
2981 REAL y8i = x1i + x3r;
2982 REAL y12r = x1r + x3i;
2983 REAL y12i = x1i - x3r;
2984 x0r = a[2] + a[18];
2985 x0i = a[3] + a[19];
2986 x1r = a[2] - a[18];
2987 x1i = a[3] - a[19];
2988 x2r = a[10] + a[26];
2989 x2i = a[11] + a[27];
2990 x3r = a[10] - a[26];
2991 x3i = a[11] - a[27];
2992 REAL y1r = x0r + x2r;
2993 REAL y1i = x0i + x2i;
2994 REAL y5r = x0r - x2r;
2995 REAL y5i = x0i - x2i;
2996 x0r = x1r - x3i;
2997 x0i = x1i + x3r;
2998 REAL y9r = wk1r * x0r - wk1i * x0i;
2999 REAL y9i = wk1r * x0i + wk1i * x0r;
3000 x0r = x1r + x3i;
3001 x0i = x1i - x3r;
3002 REAL y13r = wk1i * x0r - wk1r * x0i;
3003 REAL y13i = wk1i * x0i + wk1r * x0r;
3004 x0r = a[4] + a[20];
3005 x0i = a[5] + a[21];
3006 x1r = a[4] - a[20];
3007 x1i = a[5] - a[21];
3008 x2r = a[12] + a[28];
3009 x2i = a[13] + a[29];
3010 x3r = a[12] - a[28];
3011 x3i = a[13] - a[29];
3012 REAL y2r = x0r + x2r;
3013 REAL y2i = x0i + x2i;
3014 REAL y6r = x0r - x2r;
3015 REAL y6i = x0i - x2i;
3016 x0r = x1r - x3i;
3017 x0i = x1i + x3r;
3018 REAL y10r = wn4r * (x0r - x0i);
3019 REAL y10i = wn4r * (x0i + x0r);
3020 x0r = x1r + x3i;
3021 x0i = x1i - x3r;
3022 REAL y14r = wn4r * (x0r + x0i);
3023 REAL y14i = wn4r * (x0i - x0r);
3024 x0r = a[6] + a[22];
3025 x0i = a[7] + a[23];
3026 x1r = a[6] - a[22];
3027 x1i = a[7] - a[23];
3028 x2r = a[14] + a[30];
3029 x2i = a[15] + a[31];
3030 x3r = a[14] - a[30];
3031 x3i = a[15] - a[31];
3032 REAL y3r = x0r + x2r;
3033 REAL y3i = x0i + x2i;
3034 REAL y7r = x0r - x2r;
3035 REAL y7i = x0i - x2i;
3036 x0r = x1r - x3i;
3037 x0i = x1i + x3r;
3038 REAL y11r = wk1i * x0r - wk1r * x0i;
3039 REAL y11i = wk1i * x0i + wk1r * x0r;
3040 x0r = x1r + x3i;
3041 x0i = x1i - x3r;
3042 REAL y15r = wk1r * x0r - wk1i * x0i;
3043 REAL y15i = wk1r * x0i + wk1i * x0r;
3044 x0r = y12r - y14r;
3045 x0i = y12i - y14i;
3046 x1r = y12r + y14r;
3047 x1i = y12i + y14i;
3048 x2r = y13r - y15r;
3049 x2i = y13i - y15i;
3050 x3r = y13r + y15r;
3051 x3i = y13i + y15i;
3052 a[24] = x0r + x2r;
3053 a[25] = x0i + x2i;
3054 a[26] = x0r - x2r;
3055 a[27] = x0i - x2i;
3056 a[28] = x1r - x3i;
3057 a[29] = x1i + x3r;
3058 a[30] = x1r + x3i;
3059 a[31] = x1i - x3r;
3060 x0r = y8r + y10r;
3061 x0i = y8i + y10i;
3062 x1r = y8r - y10r;
3063 x1i = y8i - y10i;
3064 x2r = y9r + y11r;
3065 x2i = y9i + y11i;
3066 x3r = y9r - y11r;
3067 x3i = y9i - y11i;
3068 a[16] = x0r + x2r;
3069 a[17] = x0i + x2i;
3070 a[18] = x0r - x2r;
3071 a[19] = x0i - x2i;
3072 a[20] = x1r - x3i;
3073 a[21] = x1i + x3r;
3074 a[22] = x1r + x3i;
3075 a[23] = x1i - x3r;
3076 x0r = y5r - y7i;
3077 x0i = y5i + y7r;
3078 x2r = wn4r * (x0r - x0i);
3079 x2i = wn4r * (x0i + x0r);
3080 x0r = y5r + y7i;
3081 x0i = y5i - y7r;
3082 x3r = wn4r * (x0r - x0i);
3083 x3i = wn4r * (x0i + x0r);
3084 x0r = y4r - y6i;
3085 x0i = y4i + y6r;
3086 x1r = y4r + y6i;
3087 x1i = y4i - y6r;
3088 a[8] = x0r + x2r;
3089 a[9] = x0i + x2i;
3090 a[10] = x0r - x2r;
3091 a[11] = x0i - x2i;
3092 a[12] = x1r - x3i;
3093 a[13] = x1i + x3r;
3094 a[14] = x1r + x3i;
3095 a[15] = x1i - x3r;
3096 x0r = y0r + y2r;
3097 x0i = y0i + y2i;
3098 x1r = y0r - y2r;
3099 x1i = y0i - y2i;
3100 x2r = y1r + y3r;
3101 x2i = y1i + y3i;
3102 x3r = y1r - y3r;
3103 x3i = y1i - y3i;
3104 a[0] = x0r + x2r;
3105 a[1] = x0i + x2i;
3106 a[2] = x0r - x2r;
3107 a[3] = x0i - x2i;
3108 a[4] = x1r - x3i;
3109 a[5] = x1i + x3r;
3110 a[6] = x1r + x3i;
3111 a[7] = x1i - x3r;
3112 }
3113
3114
3115 void cftf162(REAL *a, REAL *w)
3116 {
3117 /*
3118 REAL wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i,
3119 x0r, x0i, x1r, x1i, x2r, x2i,
3120 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
3121 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i,
3122 y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i,
3123 y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
3124 */
3125 REAL wn4r = w[1];
3126 REAL wk1r = w[4];
3127 REAL wk1i = w[5];
3128 REAL wk3r = w[6];
3129 REAL wk3i = -w[7];
3130 REAL wk2r = w[8];
3131 REAL wk2i = w[9];
3132 REAL x1r = a[0] - a[17];
3133 REAL x1i = a[1] + a[16];
3134 REAL x0r = a[8] - a[25];
3135 REAL x0i = a[9] + a[24];
3136 REAL x2r = wn4r * (x0r - x0i);
3137 REAL x2i = wn4r * (x0i + x0r);
3138 REAL y0r = x1r + x2r;
3139 REAL y0i = x1i + x2i;
3140 REAL y4r = x1r - x2r;
3141 REAL y4i = x1i - x2i;
3142 x1r = a[0] + a[17];
3143 x1i = a[1] - a[16];
3144 x0r = a[8] + a[25];
3145 x0i = a[9] - a[24];
3146 x2r = wn4r * (x0r - x0i);
3147 x2i = wn4r * (x0i + x0r);
3148 REAL y8r = x1r - x2i;
3149 REAL y8i = x1i + x2r;
3150 REAL y12r = x1r + x2i;
3151 REAL y12i = x1i - x2r;
3152 x0r = a[2] - a[19];
3153 x0i = a[3] + a[18];
3154 x1r = wk1r * x0r - wk1i * x0i;
3155 x1i = wk1r * x0i + wk1i * x0r;
3156 x0r = a[10] - a[27];
3157 x0i = a[11] + a[26];
3158 x2r = wk3i * x0r - wk3r * x0i;
3159 x2i = wk3i * x0i + wk3r * x0r;
3160 REAL y1r = x1r + x2r;
3161 REAL y1i = x1i + x2i;
3162 REAL y5r = x1r - x2r;
3163 REAL y5i = x1i - x2i;
3164 x0r = a[2] + a[19];
3165 x0i = a[3] - a[18];
3166 x1r = wk3r * x0r - wk3i * x0i;
3167 x1i = wk3r * x0i + wk3i * x0r;
3168 x0r = a[10] + a[27];
3169 x0i = a[11] - a[26];
3170 x2r = wk1r * x0r + wk1i * x0i;
3171 x2i = wk1r * x0i - wk1i * x0r;
3172 REAL y9r = x1r - x2r;
3173 REAL y9i = x1i - x2i;
3174 REAL y13r = x1r + x2r;
3175 REAL y13i = x1i + x2i;
3176 x0r = a[4] - a[21];
3177 x0i = a[5] + a[20];
3178 x1r = wk2r * x0r - wk2i * x0i;
3179 x1i = wk2r * x0i + wk2i * x0r;
3180 x0r = a[12] - a[29];
3181 x0i = a[13] + a[28];
3182 x2r = wk2i * x0r - wk2r * x0i;
3183 x2i = wk2i * x0i + wk2r * x0r;
3184 REAL y2r = x1r + x2r;
3185 REAL y2i = x1i + x2i;
3186 REAL y6r = x1r - x2r;
3187 REAL y6i = x1i - x2i;
3188 x0r = a[4] + a[21];
3189 x0i = a[5] - a[20];
3190 x1r = wk2i * x0r - wk2r * x0i;
3191 x1i = wk2i * x0i + wk2r * x0r;
3192 x0r = a[12] + a[29];
3193 x0i = a[13] - a[28];
3194 x2r = wk2r * x0r - wk2i * x0i;
3195 x2i = wk2r * x0i + wk2i * x0r;
3196 REAL y10r = x1r - x2r;
3197 REAL y10i = x1i - x2i;
3198 REAL y14r = x1r + x2r;
3199 REAL y14i = x1i + x2i;
3200 x0r = a[6] - a[23];
3201 x0i = a[7] + a[22];
3202 x1r = wk3r * x0r - wk3i * x0i;
3203 x1i = wk3r * x0i + wk3i * x0r;
3204 x0r = a[14] - a[31];
3205 x0i = a[15] + a[30];
3206 x2r = wk1i * x0r - wk1r * x0i;
3207 x2i = wk1i * x0i + wk1r * x0r;
3208 REAL y3r = x1r + x2r;
3209 REAL y3i = x1i + x2i;
3210 REAL y7r = x1r - x2r;
3211 REAL y7i = x1i - x2i;
3212 x0r = a[6] + a[23];
3213 x0i = a[7] - a[22];
3214 x1r = wk1i * x0r + wk1r * x0i;
3215 x1i = wk1i * x0i - wk1r * x0r;
3216 x0r = a[14] + a[31];
3217 x0i = a[15] - a[30];
3218 x2r = wk3i * x0r - wk3r * x0i;
3219 x2i = wk3i * x0i + wk3r * x0r;
3220 REAL y11r = x1r + x2r;
3221 REAL y11i = x1i + x2i;
3222 REAL y15r = x1r - x2r;
3223 REAL y15i = x1i - x2i;
3224 x1r = y0r + y2r;
3225 x1i = y0i + y2i;
3226 x2r = y1r + y3r;
3227 x2i = y1i + y3i;
3228 a[0] = x1r + x2r;
3229 a[1] = x1i + x2i;
3230 a[2] = x1r - x2r;
3231 a[3] = x1i - x2i;
3232 x1r = y0r - y2r;
3233 x1i = y0i - y2i;
3234 x2r = y1r - y3r;
3235 x2i = y1i - y3i;
3236 a[4] = x1r - x2i;
3237 a[5] = x1i + x2r;
3238 a[6] = x1r + x2i;
3239 a[7] = x1i - x2r;
3240 x1r = y4r - y6i;
3241 x1i = y4i + y6r;
3242 x0r = y5r - y7i;
3243 x0i = y5i + y7r;
3244 x2r = wn4r * (x0r - x0i);
3245 x2i = wn4r * (x0i + x0r);
3246 a[8] = x1r + x2r;
3247 a[9] = x1i + x2i;
3248 a[10] = x1r - x2r;
3249 a[11] = x1i - x2i;
3250 x1r = y4r + y6i;
3251 x1i = y4i - y6r;
3252 x0r = y5r + y7i;
3253 x0i = y5i - y7r;
3254 x2r = wn4r * (x0r - x0i);
3255 x2i = wn4r * (x0i + x0r);
3256 a[12] = x1r - x2i;
3257 a[13] = x1i + x2r;
3258 a[14] = x1r + x2i;
3259 a[15] = x1i - x2r;
3260 x1r = y8r + y10r;
3261 x1i = y8i + y10i;
3262 x2r = y9r - y11r;
3263 x2i = y9i - y11i;
3264 a[16] = x1r + x2r;
3265 a[17] = x1i + x2i;
3266 a[18] = x1r - x2r;
3267 a[19] = x1i - x2i;
3268 x1r = y8r - y10r;
3269 x1i = y8i - y10i;
3270 x2r = y9r + y11r;
3271 x2i = y9i + y11i;
3272 a[20] = x1r - x2i;
3273 a[21] = x1i + x2r;
3274 a[22] = x1r + x2i;
3275 a[23] = x1i - x2r;
3276 x1r = y12r - y14i;
3277 x1i = y12i + y14r;
3278 x0r = y13r + y15i;
3279 x0i = y13i - y15r;
3280 x2r = wn4r * (x0r - x0i);
3281 x2i = wn4r * (x0i + x0r);
3282 a[24] = x1r + x2r;
3283 a[25] = x1i + x2i;
3284 a[26] = x1r - x2r;
3285 a[27] = x1i - x2i;
3286 x1r = y12r + y14i;
3287 x1i = y12i - y14r;
3288 x0r = y13r - y15i;
3289 x0i = y13i + y15r;
3290 x2r = wn4r * (x0r - x0i);
3291 x2i = wn4r * (x0i + x0r);
3292 a[28] = x1r - x2i;
3293 a[29] = x1i + x2r;
3294 a[30] = x1r + x2i;
3295 a[31] = x1i - x2r;
3296 }
3297
3298
3299 void cftf081(REAL *a, REAL *w)
3300 {
3301 /*
3302 REAL wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
3303 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
3304 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
3305 */
3306 REAL wn4r = w[1];
3307 REAL x0r = a[0] + a[8];
3308 REAL x0i = a[1] + a[9];
3309 REAL x1r = a[0] - a[8];
3310 REAL x1i = a[1] - a[9];
3311 REAL x2r = a[4] + a[12];
3312 REAL x2i = a[5] + a[13];
3313 REAL x3r = a[4] - a[12];
3314 REAL x3i = a[5] - a[13];
3315 REAL y0r = x0r + x2r;
3316 REAL y0i = x0i + x2i;
3317 REAL y2r = x0r - x2r;
3318 REAL y2i = x0i - x2i;
3319 REAL y1r = x1r - x3i;
3320 REAL y1i = x1i + x3r;
3321 REAL y3r = x1r + x3i;
3322 REAL y3i = x1i - x3r;
3323 x0r = a[2] + a[10];
3324 x0i = a[3] + a[11];
3325 x1r = a[2] - a[10];
3326 x1i = a[3] - a[11];
3327 x2r = a[6] + a[14];
3328 x2i = a[7] + a[15];
3329 x3r = a[6] - a[14];
3330 x3i = a[7] - a[15];
3331 REAL y4r = x0r + x2r;
3332 REAL y4i = x0i + x2i;
3333 REAL y6r = x0r - x2r;
3334 REAL y6i = x0i - x2i;
3335 x0r = x1r - x3i;
3336 x0i = x1i + x3r;
3337 x2r = x1r + x3i;
3338 x2i = x1i - x3r;
3339 REAL y5r = wn4r * (x0r - x0i);
3340 REAL y5i = wn4r * (x0r + x0i);
3341 REAL y7r = wn4r * (x2r - x2i);
3342 REAL y7i = wn4r * (x2r + x2i);
3343 a[8] = y1r + y5r;
3344 a[9] = y1i + y5i;
3345 a[10] = y1r - y5r;
3346 a[11] = y1i - y5i;
3347 a[12] = y3r - y7i;
3348 a[13] = y3i + y7r;
3349 a[14] = y3r + y7i;
3350 a[15] = y3i - y7r;
3351 a[0] = y0r + y4r;
3352 a[1] = y0i + y4i;
3353 a[2] = y0r - y4r;
3354 a[3] = y0i - y4i;
3355 a[4] = y2r - y6i;
3356 a[5] = y2i + y6r;
3357 a[6] = y2r + y6i;
3358 a[7] = y2i - y6r;
3359 }
3360
3361
3362 void cftf082(REAL *a, REAL *w)
3363 {
3364 /*
3365 REAL wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i,
3366 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
3367 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
3368 */
3369 REAL wn4r = w[1];
3370 REAL wk1r = w[2];
3371 REAL wk1i = w[3];
3372 REAL y0r = a[0] - a[9];
3373 REAL y0i = a[1] + a[8];
3374 REAL y1r = a[0] + a[9];
3375 REAL y1i = a[1] - a[8];
3376 REAL x0r = a[4] - a[13];
3377 REAL x0i = a[5] + a[12];
3378 REAL y2r = wn4r * (x0r - x0i);
3379 REAL y2i = wn4r * (x0i + x0r);
3380 x0r = a[4] + a[13];
3381 x0i = a[5] - a[12];
3382 REAL y3r = wn4r * (x0r - x0i);
3383 REAL y3i = wn4r * (x0i + x0r);
3384 x0r = a[2] - a[11];
3385 x0i = a[3] + a[10];
3386 REAL y4r = wk1r * x0r - wk1i * x0i;
3387 REAL y4i = wk1r * x0i + wk1i * x0r;
3388 x0r = a[2] + a[11];
3389 x0i = a[3] - a[10];
3390 REAL y5r = wk1i * x0r - wk1r * x0i;
3391 REAL y5i = wk1i * x0i + wk1r * x0r;
3392 x0r = a[6] - a[15];
3393 x0i = a[7] + a[14];
3394 REAL y6r = wk1i * x0r - wk1r * x0i;
3395 REAL y6i = wk1i * x0i + wk1r * x0r;
3396 x0r = a[6] + a[15];
3397 x0i = a[7] - a[14];
3398 REAL y7r = wk1r * x0r - wk1i * x0i;
3399 REAL y7i = wk1r * x0i + wk1i * x0r;
3400 x0r = y0r + y2r;
3401 x0i = y0i + y2i;
3402 REAL x1r = y4r + y6r;
3403 REAL x1i = y4i + y6i;
3404 a[0] = x0r + x1r;
3405 a[1] = x0i + x1i;
3406 a[2] = x0r - x1r;
3407 a[3] = x0i - x1i;
3408 x0r = y0r - y2r;
3409 x0i = y0i - y2i;
3410 x1r = y4r - y6r;
3411 x1i = y4i - y6i;
3412 a[4] = x0r - x1i;
3413 a[5] = x0i + x1r;
3414 a[6] = x0r + x1i;
3415 a[7] = x0i - x1r;
3416 x0r = y1r - y3i;
3417 x0i = y1i + y3r;
3418 x1r = y5r - y7r;
3419 x1i = y5i - y7i;
3420 a[8] = x0r + x1r;
3421 a[9] = x0i + x1i;
3422 a[10] = x0r - x1r;
3423 a[11] = x0i - x1i;
3424 x0r = y1r + y3i;
3425 x0i = y1i - y3r;
3426 x1r = y5r + y7r;
3427 x1i = y5i + y7i;
3428 a[12] = x0r - x1i;
3429 a[13] = x0i + x1r;
3430 a[14] = x0r + x1i;
3431 a[15] = x0i - x1r;
3432 }
3433
3434
3435 void cftf040(REAL *a)
3436 {
3437 //REAL x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
3438
3439 immutable REAL x0r = a[0] + a[4];
3440 immutable REAL x0i = a[1] + a[5];
3441 immutable REAL x1r = a[0] - a[4];
3442 immutable REAL x1i = a[1] - a[5];
3443 immutable REAL x2r = a[2] + a[6];
3444 immutable REAL x2i = a[3] + a[7];
3445 immutable REAL x3r = a[2] - a[6];
3446 immutable REAL x3i = a[3] - a[7];
3447 a[0] = x0r + x2r;
3448 a[1] = x0i + x2i;
3449 a[2] = x1r - x3i;
3450 a[3] = x1i + x3r;
3451 a[4] = x0r - x2r;
3452 a[5] = x0i - x2i;
3453 a[6] = x1r + x3i;
3454 a[7] = x1i - x3r;
3455 }
3456
3457
3458 void cftb040(REAL *a)
3459 {
3460 //REAL x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
3461
3462 immutable REAL x0r = a[0] + a[4];
3463 immutable REAL x0i = a[1] + a[5];
3464 immutable REAL x1r = a[0] - a[4];
3465 immutable REAL x1i = a[1] - a[5];
3466 immutable REAL x2r = a[2] + a[6];
3467 immutable REAL x2i = a[3] + a[7];
3468 immutable REAL x3r = a[2] - a[6];
3469 immutable REAL x3i = a[3] - a[7];
3470 a[0] = x0r + x2r;
3471 a[1] = x0i + x2i;
3472 a[2] = x1r + x3i;
3473 a[3] = x1i - x3r;
3474 a[4] = x0r - x2r;
3475 a[5] = x0i - x2i;
3476 a[6] = x1r - x3i;
3477 a[7] = x1i + x3r;
3478 }
3479
3480
3481 void cftx020(REAL *a)
3482 {
3483 //REAL x0r, x0i;
3484
3485 immutable REAL x0r = a[0] - a[2];
3486 immutable REAL x0i = a[1] - a[3];
3487 a[0] += a[2];
3488 a[1] += a[3];
3489 a[2] = x0r;
3490 a[3] = x0i;
3491 }
3492
3493
3494 void rftfsub(int n, REAL *a, int nc, REAL *c)
3495 {
3496 int j, k, kk, ks, m;
3497 //REAL wkr, wki, xr, xi, yr, yi;
3498
3499 m = n >> 1;
3500 ks = 2 * nc / m;
3501 kk = 0;
3502 for (j = 2; j < m; j += 2) {
3503 k = n - j;
3504 kk += ks;
3505 REAL wkr = 0.5 - c[nc - kk];
3506 REAL wki = c[kk];
3507 REAL xr = a[j] - a[k];
3508 REAL xi = a[j + 1] + a[k + 1];
3509 REAL yr = wkr * xr - wki * xi;
3510 REAL yi = wkr * xi + wki * xr;
3511 a[j] -= yr;
3512 a[j + 1] -= yi;
3513 a[k] += yr;
3514 a[k + 1] -= yi;
3515 }
3516 }
3517
3518
3519 void rftbsub(int n, REAL *a, int nc, REAL *c)
3520 {
3521 int j, k, kk, ks, m;
3522 //REAL wkr, wki, xr, xi, yr, yi;
3523
3524 m = n >> 1;
3525 ks = 2 * nc / m;
3526 kk = 0;
3527 for (j = 2; j < m; j += 2) {
3528 k = n - j;
3529 kk += ks;
3530 REAL wkr = 0.5 - c[nc - kk];
3531 REAL wki = c[kk];
3532 REAL xr = a[j] - a[k];
3533 REAL xi = a[j + 1] + a[k + 1];
3534 REAL yr = wkr * xr + wki * xi;
3535 REAL yi = wkr * xi - wki * xr;
3536 a[j] -= yr;
3537 a[j + 1] -= yi;
3538 a[k] += yr;
3539 a[k + 1] -= yi;
3540 }
3541 }
3542
3543
3544 void dctsub(int n, REAL *a, int nc, REAL *c)
3545 {
3546 int j, k, kk, ks, m;
3547 //REAL wkr, wki, xr;
3548
3549 m = n >> 1;
3550 ks = nc / n;
3551 kk = 0;
3552 for (j = 1; j < m; j++) {
3553 k = n - j;
3554 kk += ks;
3555 REAL wkr = c[kk] - c[nc - kk];
3556 REAL wki = c[kk] + c[nc - kk];
3557 REAL xr = wki * a[j] - wkr * a[k];
3558 a[j] = wkr * a[j] + wki * a[k];
3559 a[k] = xr;
3560 }
3561 a[m] *= c[0];
3562 }
3563
3564
3565 void dstsub(int n, REAL *a, int nc, REAL *c)
3566 {
3567 int j, k, kk, ks, m;
3568 //REAL wkr, wki, xr;
3569
3570 m = n >> 1;
3571 ks = nc / n;
3572 kk = 0;
3573 for (j = 1; j < m; j++) {
3574 k = n - j;
3575 kk += ks;
3576 REAL wkr = c[kk] - c[nc - kk];
3577 REAL wki = c[kk] + c[nc - kk];
3578 REAL xr = wki * a[k] - wkr * a[j];
3579 a[k] = wkr * a[k] + wki * a[j];
3580 a[j] = xr;
3581 }
3582 a[m] *= c[0];
3583 }
various I.V. D modules