WDL/fft.c

   1 /*
   2   WDL - fft.cpp
   3   Copyright (C) 2006 and later Cockos Incorporated
   4   Copyright 1999 D. J. Bernstein
   5
   6   This software is provided 'as-is', without any express or implied
   7   warranty.  In no event will the authors be held liable for any damages
   8   arising from the use of this software.
   9
  10   Permission is granted to anyone to use this software for any purpose,
  11   including commercial applications, and to alter it and redistribute it
  12   freely, subject to the following restrictions:
  13
  14   1. The origin of this software must not be misrepresented; you must not
  15      claim that you wrote the original software. If you use this software
  16      in a product, an acknowledgment in the product documentation would be
  17      appreciated but is not required.
  18   2. Altered source versions must be plainly marked as such, and must not be
  19      misrepresented as being the original software.
  20   3. This notice may not be removed or altered from any source distribution.
  21
  22
  23
  24   This file implements the WDL FFT library. These routines are based on the
  25   DJBFFT library, which are   Copyright 1999 D. J. Bernstein, djb@pobox.com
  26
  27   The DJB FFT web page is:  http://cr.yp.to/djbfft.html
  28
  29 */
  30
  31
  32 // this is based on djbfft
  33
  34 #include <math.h>
  35 #include "fft.h"
  36
  37
  38 #define FFT_MAXBITLEN 15
  39
  40 #ifdef _MSC_VER
  41 #define inline __inline
  42 #endif
  43
  44 #define PI 3.1415926535897932384626433832795
  45
  46 static WDL_FFT_COMPLEX d16[3];
  47 static WDL_FFT_COMPLEX d32[7];
  48 static WDL_FFT_COMPLEX d64[15];
  49 static WDL_FFT_COMPLEX d128[31];
  50 static WDL_FFT_COMPLEX d256[63];
  51 static WDL_FFT_COMPLEX d512[127];
  52 static WDL_FFT_COMPLEX d1024[127];
  53 static WDL_FFT_COMPLEX d2048[255];
  54 static WDL_FFT_COMPLEX d4096[511];
  55 static WDL_FFT_COMPLEX d8192[1023];
  56 static WDL_FFT_COMPLEX d16384[2047];
  57 static WDL_FFT_COMPLEX d32768[4095];
  58
  59
  60 #define sqrthalf (d16[1].re)
  61
  62 #define VOL *(volatile WDL_FFT_REAL *)&
  63
  64 #define TRANSFORM(a0,a1,a2,a3,wre,wim) { \
  65   t6 = a2.re; \
  66   t1 = a0.re - t6; \
  67   t6 += a0.re; \
  68   a0.re = t6; \
  69   t3 = a3.im; \
  70   t4 = a1.im - t3; \
  71   t8 = t1 - t4; \
  72   t1 += t4; \
  73   t3 += a1.im; \
  74   a1.im = t3; \
  75   t5 = wre; \
  76   t7 = t8 * t5; \
  77   t4 = t1 * t5; \
  78   t8 *= wim; \
  79   t2 = a3.re; \
  80   t3 = a1.re - t2; \
  81   t2 += a1.re; \
  82   a1.re = t2; \
  83   t1 *= wim; \
  84   t6 = a2.im; \
  85   t2 = a0.im - t6; \
  86   t6 += a0.im; \
  87   a0.im = t6; \
  88   t6 = t2 + t3; \
  89   t2 -= t3; \
  90   t3 = t6 * wim; \
  91   t7 -= t3; \
  92   a2.re = t7; \
  93   t6 *= t5; \
  94   t6 += t8; \
  95   a2.im = t6; \
  96   t5 *= t2; \
  97   t5 -= t1; \
  98   a3.im = t5; \
  99   t2 *= wim; \
 100   t4 += t2; \
 101   a3.re = t4; \
 102   }
 103
 104 #define TRANSFORMHALF(a0,a1,a2,a3) { \
 105   t1 = a2.re; \
 106   t5 = a0.re - t1; \
 107   t1 += a0.re; \
 108   a0.re = t1; \
 109   t4 = a3.im; \
 110   t8 = a1.im - t4; \
 111   t1 = t5 - t8; \
 112   t5 += t8; \
 113   t4 += a1.im; \
 114   a1.im = t4; \
 115   t3 = a3.re; \
 116   t7 = a1.re - t3; \
 117   t3 += a1.re; \
 118   a1.re = t3; \
 119   t8 = a2.im; \
 120   t6 = a0.im - t8; \
 121   t2 = t6 + t7; \
 122   t6 -= t7; \
 123   t8 += a0.im; \
 124   a0.im = t8; \
 125   t4 = t6 + t5; \
 126   t3 = sqrthalf; \
 127   t4 *= t3; \
 128   a3.re = t4; \
 129   t6 -= t5; \
 130   t6 *= t3; \
 131   a3.im = t6; \
 132   t7 = t1 - t2; \
 133   t7 *= t3; \
 134   a2.re = t7; \
 135   t2 += t1; \
 136   t2 *= t3; \
 137   a2.im = t2; \
 138   }
 139
 140 #define TRANSFORMZERO(a0,a1,a2,a3) { \
 141   t5 = a2.re; \
 142   t1 = a0.re - t5; \
 143   t5 += a0.re; \
 144   a0.re = t5; \
 145   t8 = a3.im; \
 146   t4 = a1.im - t8; \
 147   t7 = a3.re; \
 148   t6 = t1 - t4; \
 149   a2.re = t6; \
 150   t1 += t4; \
 151   a3.re = t1; \
 152   t8 += a1.im; \
 153   a1.im = t8; \
 154   t3 = a1.re - t7; \
 155   t7 += a1.re; \
 156   a1.re = t7; \
 157   t6 = a2.im; \
 158   t2 = a0.im - t6; \
 159   t7 = t2 + t3; \
 160   a2.im = t7; \
 161   t2 -= t3; \
 162   a3.im = t2; \
 163   t6 += a0.im; \
 164   a0.im = t6; \
 165   }
 166
 167 #define UNTRANSFORM(a0,a1,a2,a3,wre,wim) { \
 168   t6 = VOL wre; \
 169   t1 = VOL a2.re; \
 170   t1 *= t6; \
 171   t8 = VOL wim; \
 172   t3 = VOL a2.im; \
 173   t3 *= t8; \
 174   t2 = VOL a2.im; \
 175   t4 = VOL a2.re; \
 176   t5 = VOL a3.re; \
 177   t5 *= t6; \
 178   t7 = VOL a3.im; \
 179   t1 += t3; \
 180   t7 *= t8; \
 181   t5 -= t7; \
 182   t3 = t5 + t1; \
 183   t5 -= t1; \
 184   t2 *= t6; \
 185   t6 *= a3.im; \
 186   t4 *= t8; \
 187   t2 -= t4; \
 188   t8 *= a3.re; \
 189   t6 += t8; \
 190   t1 = a0.re - t3; \
 191   t3 += a0.re; \
 192   a0.re = t3; \
 193   t7 = a1.im - t5; \
 194   t5 += a1.im; \
 195   a1.im = t5; \
 196   t4 = t2 - t6; \
 197   t6 += t2; \
 198   t8 = a1.re - t4; \
 199   t4 += a1.re; \
 200   a1.re = t4; \
 201   t2 = a0.im - t6; \
 202   t6 += a0.im; \
 203   a0.im = t6; \
 204   a2.re = t1; \
 205   a3.im = t7; \
 206   a3.re = t8; \
 207   a2.im = t2; \
 208   }
 209
 210
 211 #define UNTRANSFORMHALF(a0,a1,a2,a3) { \
 212   t6 = sqrthalf; \
 213   t1 = a2.re; \
 214   t2 = a2.im - t1; \
 215   t2 *= t6; \
 216   t1 += a2.im; \
 217   t1 *= t6; \
 218   t4 = a3.im; \
 219   t3 = a3.re - t4; \
 220   t3 *= t6; \
 221   t4 += a3.re; \
 222   t4 *= t6; \
 223   t8 = t3 - t1; \
 224   t7 = t2 - t4; \
 225   t1 += t3; \
 226   t2 += t4; \
 227   t4 = a1.im - t8; \
 228   a3.im = t4; \
 229   t8 += a1.im; \
 230   a1.im = t8; \
 231   t3 = a1.re - t7; \
 232   a3.re = t3; \
 233   t7 += a1.re; \
 234   a1.re = t7; \
 235   t5 = a0.re - t1; \
 236   a2.re = t5; \
 237   t1 += a0.re; \
 238   a0.re = t1; \
 239   t6 = a0.im - t2; \
 240   a2.im = t6; \
 241   t2 += a0.im; \
 242   a0.im = t2; \
 243   }
 244
 245 #define UNTRANSFORMZERO(a0,a1,a2,a3) { \
 246   t2 = a3.im; \
 247   t3 = a2.im - t2; \
 248   t2 += a2.im; \
 249   t1 = a2.re; \
 250   t4 = a3.re - t1; \
 251   t1 += a3.re; \
 252   t5 = a0.re - t1; \
 253   a2.re = t5; \
 254   t6 = a0.im - t2; \
 255   a2.im = t6; \
 256   t7 = a1.re - t3; \
 257   a3.re = t7; \
 258   t8 = a1.im - t4; \
 259   a3.im = t8; \
 260   t1 += a0.re; \
 261   a0.re = t1; \
 262   t2 += a0.im; \
 263   a0.im = t2; \
 264   t3 += a1.re; \
 265   a1.re = t3; \
 266   t4 += a1.im; \
 267   a1.im = t4; \
 268   }
 269
 270 #define R(a0,a1,b0,b1,wre,wim) { \
 271   t1 = a0 - a1; \
 272   t2 = b0 - b1; \
 273   t5 = t1 * wim; \
 274   t6 = t2 * wim; \
 275   t3 = VOL a0; \
 276   t1 *= wre; \
 277   t3 += a1; \
 278   t2 *= wre; \
 279   t1 -= t6; \
 280   t4 = VOL b0; \
 281   t2 += t5; \
 282   t4 += b1; \
 283   a0 = t3; \
 284   b1 = t2; \
 285   a1 = t4; \
 286   b0 = t1; \
 287   }
 288
 289 #define RHALF(a0,a1,b0,b1) { \
 290   t1 = a0 - a1; \
 291   t2 = b0 - b1; \
 292   t3 = a0 + a1; \
 293   t5 = t1 - t2; \
 294   t1 += t2; \
 295   t4 = VOL b0; \
 296   t5 *= sqrthalf; \
 297   t4 += b1; \
 298   t1 *= sqrthalf; \
 299   a0 = t3; \
 300   b1 = t1; \
 301   a1 = t4; \
 302   b0 = t5; \
 303   }
 304
 305 #define RZERO(a0,a1,b0,b1) { \
 306   t1 = a0 - a1; \
 307   t2 = b0 - b1; \
 308   t3 = a0 + a1; \
 309   t4 = b0 + b1; \
 310   b0 = t1; \
 311   a0 = t3; \
 312   b1 = t2; \
 313   a1 = t4; \
 314   }
 315
 316 #define V(a0,a1,b0,b1,wre,wim) { \
 317   t5 = b0 * wre; \
 318   t1 = b1 * wim; \
 319   t6 = b1 * wre; \
 320   t5 += t1; \
 321   t3 = b0 * wim; \
 322   t2 = a0 - t5; \
 323   t6 -= t3; \
 324   t5 += a0; \
 325   t4 = a1 - t6; \
 326   t6 += a1; \
 327   a1 = t2; \
 328   a0 = t5; \
 329   b1 = t4; \
 330   b0 = t6; \
 331   }
 332
 333 #define VHALF(a0,a1,b0,b1) { \
 334   t5 = b0 + b1; \
 335   t6 = b1 - b0; \
 336   t5 *= sqrthalf; \
 337   t2 = VOL a0; \
 338   t6 *= sqrthalf; \
 339   t2 -= t5; \
 340   t5 += a0; \
 341   t4 = a1 - t6; \
 342   t6 += a1; \
 343   a1 = t2; \
 344   a0 = t5; \
 345   b0 = t6; \
 346   b1 = t4; \
 347   }
 348
 349 #define VZERO(a0,a1,b0,b1) { \
 350   t1 = a0 + b0; \
 351   t2 = a0 - b0; \
 352   t3 = a1 + b1; \
 353   t4 = a1 - b1; \
 354   a0 = t1; \
 355   b0 = t3; \
 356   a1 = t2; \
 357   b1 = t4; \
 358   }
 359
 360 static void c2(register WDL_FFT_COMPLEX *a)
 361 {
 362   register WDL_FFT_REAL t1;
 363
 364   t1 = a[1].re;
 365   a[1].re = a[0].re - t1;
 366   a[0].re += t1;
 367
 368   t1 = a[1].im;
 369   a[1].im = a[0].im - t1;
 370   a[0].im += t1;
 371 }
 372
 373 static inline void c4(register WDL_FFT_COMPLEX *a)
 374 {
 375   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8;
 376
 377   t5 = a[2].re;
 378   t1 = a[0].re - t5;
 379   t7 = a[3].re;
 380   t5 += a[0].re;
 381   t3 = a[1].re - t7;
 382   t7 += a[1].re;
 383   t8 = t5 + t7;
 384   a[0].re = t8;
 385   t5 -= t7;
 386   a[1].re = t5;
 387   t6 = a[2].im;
 388   t2 = a[0].im - t6;
 389   t6 += a[0].im;
 390   t5 = a[3].im;
 391   a[2].im = t2 + t3;
 392   t2 -= t3;
 393   a[3].im = t2;
 394   t4 = a[1].im - t5;
 395   a[3].re = t1 + t4;
 396   t1 -= t4;
 397   a[2].re = t1;
 398   t5 += a[1].im;
 399   a[0].im = t6 + t5;
 400   t6 -= t5;
 401   a[1].im = t6;
 402 }
 403
 404 static void c8(register WDL_FFT_COMPLEX *a)
 405 {
 406   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8;
 407
 408   t7 = a[4].im;
 409   t4 = a[0].im - t7;
 410   t7 += a[0].im;
 411   a[0].im = t7;
 412
 413   t8 = a[6].re;
 414   t5 = a[2].re - t8;
 415   t8 += a[2].re;
 416   a[2].re = t8;
 417
 418   t7 = a[6].im;
 419   a[6].im = t4 - t5;
 420   t4 += t5;
 421   a[4].im = t4;
 422
 423   t6 = a[2].im - t7;
 424   t7 += a[2].im;
 425   a[2].im = t7;
 426
 427   t8 = a[4].re;
 428   t3 = a[0].re - t8;
 429   t8 += a[0].re;
 430   a[0].re = t8;
 431
 432   a[4].re = t3 - t6;
 433   t3 += t6;
 434   a[6].re = t3;
 435
 436   t7 = a[5].re;
 437   t3 = a[1].re - t7;
 438   t7 += a[1].re;
 439   a[1].re = t7;
 440
 441   t8 = a[7].im;
 442   t6 = a[3].im - t8;
 443   t8 += a[3].im;
 444   a[3].im = t8;
 445   t1 = t3 - t6;
 446   t3 += t6;
 447
 448   t7 = a[5].im;
 449   t4 = a[1].im - t7;
 450   t7 += a[1].im;
 451   a[1].im = t7;
 452
 453   t8 = a[7].re;
 454   t5 = a[3].re - t8;
 455   t8 += a[3].re;
 456   a[3].re = t8;
 457
 458   t2 = t4 - t5;
 459   t4 += t5;
 460
 461   t6 = t1 - t4;
 462   t8 = sqrthalf;
 463   t6 *= t8;
 464   a[5].re = a[4].re - t6;
 465   t1 += t4;
 466   t1 *= t8;
 467   a[5].im = a[4].im - t1;
 468   t6 += a[4].re;
 469   a[4].re = t6;
 470   t1 += a[4].im;
 471   a[4].im = t1;
 472
 473   t5 = t2 - t3;
 474   t5 *= t8;
 475   a[7].im = a[6].im - t5;
 476   t2 += t3;
 477   t2 *= t8;
 478   a[7].re = a[6].re - t2;
 479   t2 += a[6].re;
 480   a[6].re = t2;
 481   t5 += a[6].im;
 482   a[6].im = t5;
 483
 484   c4(a);
 485 }
 486
 487 static void c16(register WDL_FFT_COMPLEX *a)
 488 {
 489   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8;
 490
 491   TRANSFORMZERO(a[0],a[4],a[8],a[12]);
 492   TRANSFORM(a[1],a[5],a[9],a[13],d16[0].re,d16[0].im);
 493   TRANSFORMHALF(a[2],a[6],a[10],a[14]);
 494   TRANSFORM(a[3],a[7],a[11],a[15],d16[0].im,d16[0].re);
 495   c4(a + 8);
 496   c4(a + 12);
 497
 498   c8(a);
 499 }
 500
 501 /* a[0...8n-1], w[0...2n-2]; n >= 2 */
 502 static void cpass(register WDL_FFT_COMPLEX *a,register const WDL_FFT_COMPLEX *w,register unsigned int n)
 503 {
 504   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8;
 505   register WDL_FFT_COMPLEX *a1;
 506   register WDL_FFT_COMPLEX *a2;
 507   register WDL_FFT_COMPLEX *a3;
 508
 509   a2 = a + 4 * n;
 510   a1 = a + 2 * n;
 511   a3 = a2 + 2 * n;
 512   --n;
 513
 514   TRANSFORMZERO(a[0],a1[0],a2[0],a3[0]);
 515   TRANSFORM(a[1],a1[1],a2[1],a3[1],w[0].re,w[0].im);
 516
 517   for (;;) {
 518     TRANSFORM(a[2],a1[2],a2[2],a3[2],w[1].re,w[1].im);
 519     TRANSFORM(a[3],a1[3],a2[3],a3[3],w[2].re,w[2].im);
 520     if (!--n) break;
 521     a += 2;
 522     a1 += 2;
 523     a2 += 2;
 524     a3 += 2;
 525     w += 2;
 526   }
 527 }
 528
 529 static void c32(register WDL_FFT_COMPLEX *a)
 530 {
 531   cpass(a,d32,4);
 532   c8(a + 16);
 533   c8(a + 24);
 534   c16(a);
 535 }
 536
 537 static void c64(register WDL_FFT_COMPLEX *a)
 538 {
 539   cpass(a,d64,8);
 540   c16(a + 32);
 541   c16(a + 48);
 542   c32(a);
 543 }
 544
 545 static void c128(register WDL_FFT_COMPLEX *a)
 546 {
 547   cpass(a,d128,16);
 548   c32(a + 64);
 549   c32(a + 96);
 550   c64(a);
 551 }
 552
 553 static void c256(register WDL_FFT_COMPLEX *a)
 554 {
 555   cpass(a,d256,32);
 556   c64(a + 128);
 557   c64(a + 192);
 558   c128(a);
 559 }
 560
 561 static void c512(register WDL_FFT_COMPLEX *a)
 562 {
 563   cpass(a,d512,64);
 564   c128(a + 384);
 565   c128(a + 256);
 566   c256(a);
 567 }
 568
 569 /* a[0...8n-1], w[0...n-2]; n even, n >= 4 */
 570 static void cpassbig(register WDL_FFT_COMPLEX *a,register const WDL_FFT_COMPLEX *w,register unsigned int n)
 571 {
 572   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8;
 573   register WDL_FFT_COMPLEX *a1;
 574   register WDL_FFT_COMPLEX *a2;
 575   register WDL_FFT_COMPLEX *a3;
 576   register unsigned int k;
 577
 578   a2 = a + 4 * n;
 579   a1 = a + 2 * n;
 580   a3 = a2 + 2 * n;
 581   k = n - 2;
 582
 583   TRANSFORMZERO(a[0],a1[0],a2[0],a3[0]);
 584   TRANSFORM(a[1],a1[1],a2[1],a3[1],w[0].re,w[0].im);
 585   a += 2;
 586   a1 += 2;
 587   a2 += 2;
 588   a3 += 2;
 589
 590   do {
 591     TRANSFORM(a[0],a1[0],a2[0],a3[0],w[1].re,w[1].im);
 592     TRANSFORM(a[1],a1[1],a2[1],a3[1],w[2].re,w[2].im);
 593     a += 2;
 594     a1 += 2;
 595     a2 += 2;
 596     a3 += 2;
 597     w += 2;
 598   } while (k -= 2);
 599
 600   TRANSFORMHALF(a[0],a1[0],a2[0],a3[0]);
 601   TRANSFORM(a[1],a1[1],a2[1],a3[1],w[0].im,w[0].re);
 602   a += 2;
 603   a1 += 2;
 604   a2 += 2;
 605   a3 += 2;
 606
 607   k = n - 2;
 608   do {
 609     TRANSFORM(a[0],a1[0],a2[0],a3[0],w[-1].im,w[-1].re);
 610     TRANSFORM(a[1],a1[1],a2[1],a3[1],w[-2].im,w[-2].re);
 611     a += 2;
 612     a1 += 2;
 613     a2 += 2;
 614     a3 += 2;
 615     w -= 2;
 616   } while (k -= 2);
 617 }
 618
 619
 620 static void c1024(register WDL_FFT_COMPLEX *a)
 621 {
 622   cpassbig(a,d1024,128);
 623   c256(a + 768);
 624   c256(a + 512);
 625   c512(a);
 626 }
 627
 628 static void c2048(register WDL_FFT_COMPLEX *a)
 629 {
 630   cpassbig(a,d2048,256);
 631   c512(a + 1536);
 632   c512(a + 1024);
 633   c1024(a);
 634 }
 635
 636 static void c4096(register WDL_FFT_COMPLEX *a)
 637 {
 638   cpassbig(a,d4096,512);
 639   c1024(a + 3072);
 640   c1024(a + 2048);
 641   c2048(a);
 642 }
 643
 644 static void c8192(register WDL_FFT_COMPLEX *a)
 645 {
 646   cpassbig(a,d8192,1024);
 647   c2048(a + 6144);
 648   c2048(a + 4096);
 649   c4096(a);
 650 }
 651
 652 static void c16384(register WDL_FFT_COMPLEX *a)
 653 {
 654   cpassbig(a,d16384,2048);
 655   c4096(a + 8192 + 4096);
 656   c4096(a + 8192);
 657   c8192(a);
 658 }
 659
 660 static void c32768(register WDL_FFT_COMPLEX *a)
 661 {
 662   cpassbig(a,d32768,4096);
 663   c8192(a + 16384 + 8192);
 664   c8192(a + 16384);
 665   c16384(a);
 666 }
 667
 668 #if 0
 669 static void mulr4(WDL_FFT_REAL *a,WDL_FFT_REAL *b)
 670 {
 671   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6;
 672
 673   t1 = a[2] * b[2];
 674   t2 = a[3] * b[3];
 675   t3 = a[3] * b[2];
 676   t4 = a[2] * b[3];
 677   t5 = a[0] * b[0];
 678   t6 = a[1] * b[1];
 679   t1 -= t2;
 680   t3 += t4;
 681   a[0] = t5;
 682   a[1] = t6;
 683   a[2] = t1;
 684   a[3] = t3;
 685 }
 686
 687 #endif
 688 /* n even, n > 0 */
 689 void WDL_fft_complexmul(WDL_FFT_COMPLEX *a,WDL_FFT_COMPLEX *b,int n)
 690 {
 691   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8;
 692   if (n<2 || (n&1)) return;
 693
 694   do {
 695     t1 = a[0].re * b[0].re;
 696     t2 = a[0].im * b[0].im;
 697     t3 = a[0].im * b[0].re;
 698     t4 = a[0].re * b[0].im;
 699     t5 = a[1].re * b[1].re;
 700     t6 = a[1].im * b[1].im;
 701     t7 = a[1].im * b[1].re;
 702     t8 = a[1].re * b[1].im;
 703     t1 -= t2;
 704     t3 += t4;
 705     t5 -= t6;
 706     t7 += t8;
 707     a[0].re = t1;
 708     a[1].re = t5;
 709     a[0].im = t3;
 710     a[1].im = t7;
 711     a += 2;
 712     b += 2;
 713   } while (n -= 2);
 714 }
 715
 716 void WDL_fft_complexmul2(WDL_FFT_COMPLEX *c, WDL_FFT_COMPLEX *a, WDL_FFT_COMPLEX *b, int n)
 717 {
 718   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8;
 719   if (n<2 || (n&1)) return;
 720
 721   do {
 722     t1 = a[0].re * b[0].re;
 723     t2 = a[0].im * b[0].im;
 724     t3 = a[0].im * b[0].re;
 725     t4 = a[0].re * b[0].im;
 726     t5 = a[1].re * b[1].re;
 727     t6 = a[1].im * b[1].im;
 728     t7 = a[1].im * b[1].re;
 729     t8 = a[1].re * b[1].im;
 730     t1 -= t2;
 731     t3 += t4;
 732     t5 -= t6;
 733     t7 += t8;
 734     c[0].re = t1;
 735     c[1].re = t5;
 736     c[0].im = t3;
 737     c[1].im = t7;
 738     a += 2;
 739     b += 2;
 740     c += 2;
 741   } while (n -= 2);
 742 }
 743 void WDL_fft_complexmul3(WDL_FFT_COMPLEX *c, WDL_FFT_COMPLEX *a, WDL_FFT_COMPLEX *b, int n)
 744 {
 745   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8;
 746   if (n<2 || (n&1)) return;
 747
 748   do {
 749     t1 = a[0].re * b[0].re;
 750     t2 = a[0].im * b[0].im;
 751     t3 = a[0].im * b[0].re;
 752     t4 = a[0].re * b[0].im;
 753     t5 = a[1].re * b[1].re;
 754     t6 = a[1].im * b[1].im;
 755     t7 = a[1].im * b[1].re;
 756     t8 = a[1].re * b[1].im;
 757     t1 -= t2;
 758     t3 += t4;
 759     t5 -= t6;
 760     t7 += t8;
 761     c[0].re += t1;
 762     c[1].re += t5;
 763     c[0].im += t3;
 764     c[1].im += t7;
 765     a += 2;
 766     b += 2;
 767     c += 2;
 768   } while (n -= 2);
 769 }
 770
 771
 772 static inline void u4(register WDL_FFT_COMPLEX *a)
 773 {
 774   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8;
 775
 776   t1 = VOL a[1].re;
 777   t3 = a[0].re - t1;
 778   t6 = VOL a[2].re;
 779   t1 += a[0].re;
 780   t8 = a[3].re - t6;
 781   t6 += a[3].re;
 782   a[0].re = t1 + t6;
 783   t1 -= t6;
 784   a[2].re = t1;
 785
 786   t2 = VOL a[1].im;
 787   t4 = a[0].im - t2;
 788   t2 += a[0].im;
 789   t5 = VOL a[3].im;
 790   a[1].im = t4 + t8;
 791   t4 -= t8;
 792   a[3].im = t4;
 793
 794   t7 = a[2].im - t5;
 795   t5 += a[2].im;
 796   a[1].re = t3 + t7;
 797   t3 -= t7;
 798   a[3].re = t3;
 799   a[0].im = t2 + t5;
 800   t2 -= t5;
 801   a[2].im = t2;
 802 }
 803
 804 static void u8(register WDL_FFT_COMPLEX *a)
 805 {
 806   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8;
 807
 808   u4(a);
 809
 810   t1 = a[5].re;
 811   a[5].re = a[4].re - t1;
 812   t1 += a[4].re;
 813
 814   t3 = a[7].re;
 815   a[7].re = a[6].re - t3;
 816   t3 += a[6].re;
 817
 818   t8 = t3 - t1;
 819   t1 += t3;
 820
 821   t6 = a[2].im - t8;
 822   t8 += a[2].im;
 823   a[2].im = t8;
 824
 825   t5 = a[0].re - t1;
 826   a[4].re = t5;
 827   t1 += a[0].re;
 828   a[0].re = t1;
 829
 830   t2 = a[5].im;
 831   a[5].im = a[4].im - t2;
 832   t2 += a[4].im;
 833
 834   t4 = a[7].im;
 835   a[7].im = a[6].im - t4;
 836   t4 += a[6].im;
 837
 838   a[6].im = t6;
 839
 840   t7 = t2 - t4;
 841   t2 += t4;
 842
 843   t3 = a[2].re - t7;
 844   a[6].re = t3;
 845   t7 += a[2].re;
 846   a[2].re = t7;
 847
 848   t6 = a[0].im - t2;
 849   a[4].im = t6;
 850   t2 += a[0].im;
 851   a[0].im = t2;
 852
 853   t6 = sqrthalf;
 854
 855   t1 = a[5].re;
 856   t2 = a[5].im - t1;
 857   t2 *= t6;
 858   t1 += a[5].im;
 859   t1 *= t6;
 860   t4 = a[7].im;
 861   t3 = a[7].re - t4;
 862   t3 *= t6;
 863   t4 += a[7].re;
 864   t4 *= t6;
 865
 866   t8 = t3 - t1;
 867   t1 += t3;
 868   t7 = t2 - t4;
 869   t2 += t4;
 870
 871   t4 = a[3].im - t8;
 872   a[7].im = t4;
 873   t5 = a[1].re - t1;
 874   a[5].re = t5;
 875   t3 = a[3].re - t7;
 876   a[7].re = t3;
 877   t6 = a[1].im - t2;
 878   a[5].im = t6;
 879
 880   t8 += a[3].im;
 881   a[3].im = t8;
 882   t1 += a[1].re;
 883   a[1].re = t1;
 884   t7 += a[3].re;
 885   a[3].re = t7;
 886   t2 += a[1].im;
 887   a[1].im = t2;
 888 }
 889
 890 static void u16(register WDL_FFT_COMPLEX *a)
 891 {
 892   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8;
 893
 894   u8(a);
 895   u4(a + 8);
 896   u4(a + 12);
 897
 898   UNTRANSFORMZERO(a[0],a[4],a[8],a[12]);
 899   UNTRANSFORMHALF(a[2],a[6],a[10],a[14]);
 900   UNTRANSFORM(a[1],a[5],a[9],a[13],d16[0].re,d16[0].im);
 901   UNTRANSFORM(a[3],a[7],a[11],a[15],d16[0].im,d16[0].re);
 902 }
 903
 904 /* a[0...8n-1], w[0...2n-2] */
 905 static void upass(register WDL_FFT_COMPLEX *a,register const WDL_FFT_COMPLEX *w,register unsigned int n)
 906 {
 907   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8;
 908   register WDL_FFT_COMPLEX *a1;
 909   register WDL_FFT_COMPLEX *a2;
 910   register WDL_FFT_COMPLEX *a3;
 911
 912   a2 = a + 4 * n;
 913   a1 = a + 2 * n;
 914   a3 = a2 + 2 * n;
 915   n -= 1;
 916
 917   UNTRANSFORMZERO(a[0],a1[0],a2[0],a3[0]);
 918   UNTRANSFORM(a[1],a1[1],a2[1],a3[1],w[0].re,w[0].im);
 919
 920   for (;;) {
 921     UNTRANSFORM(a[2],a1[2],a2[2],a3[2],w[1].re,w[1].im);
 922     UNTRANSFORM(a[3],a1[3],a2[3],a3[3],w[2].re,w[2].im);
 923     if (!--n) break;
 924     a += 2;
 925     a1 += 2;
 926     a2 += 2;
 927     a3 += 2;
 928     w += 2;
 929   }
 930 }
 931
 932 static void u32(register WDL_FFT_COMPLEX *a)
 933 {
 934   u16(a);
 935   u8(a + 16);
 936   u8(a + 24);
 937   upass(a,d32,4);
 938 }
 939
 940 static void u64(register WDL_FFT_COMPLEX *a)
 941 {
 942   u32(a);
 943   u16(a + 32);
 944   u16(a + 48);
 945   upass(a,d64,8);
 946 }
 947
 948 static void u128(register WDL_FFT_COMPLEX *a)
 949 {
 950   u64(a);
 951   u32(a + 64);
 952   u32(a + 96);
 953   upass(a,d128,16);
 954 }
 955
 956 static void u256(register WDL_FFT_COMPLEX *a)
 957 {
 958   u128(a);
 959   u64(a + 128);
 960   u64(a + 192);
 961   upass(a,d256,32);
 962 }
 963
 964 static void u512(register WDL_FFT_COMPLEX *a)
 965 {
 966   u256(a);
 967   u128(a + 256);
 968   u128(a + 384);
 969   upass(a,d512,64);
 970 }
 971
 972
 973 /* a[0...8n-1], w[0...n-2]; n even, n >= 4 */
 974 static void upassbig(register WDL_FFT_COMPLEX *a,register const WDL_FFT_COMPLEX *w,register unsigned int n)
 975 {
 976   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8;
 977   register WDL_FFT_COMPLEX *a1;
 978   register WDL_FFT_COMPLEX *a2;
 979   register WDL_FFT_COMPLEX *a3;
 980   register unsigned int k;
 981
 982   a2 = a + 4 * n;
 983   a1 = a + 2 * n;
 984   a3 = a2 + 2 * n;
 985   k = n - 2;
 986
 987   UNTRANSFORMZERO(a[0],a1[0],a2[0],a3[0]);
 988   UNTRANSFORM(a[1],a1[1],a2[1],a3[1],w[0].re,w[0].im);
 989   a += 2;
 990   a1 += 2;
 991   a2 += 2;
 992   a3 += 2;
 993
 994   do {
 995     UNTRANSFORM(a[0],a1[0],a2[0],a3[0],w[1].re,w[1].im);
 996     UNTRANSFORM(a[1],a1[1],a2[1],a3[1],w[2].re,w[2].im);
 997     a += 2;
 998     a1 += 2;
 999     a2 += 2;
1000     a3 += 2;
1001     w += 2;
1002   } while (k -= 2);
1003
1004   UNTRANSFORMHALF(a[0],a1[0],a2[0],a3[0]);
1005   UNTRANSFORM(a[1],a1[1],a2[1],a3[1],w[0].im,w[0].re);
1006   a += 2;
1007   a1 += 2;
1008   a2 += 2;
1009   a3 += 2;
1010
1011   k = n - 2;
1012   do {
1013     UNTRANSFORM(a[0],a1[0],a2[0],a3[0],w[-1].im,w[-1].re);
1014     UNTRANSFORM(a[1],a1[1],a2[1],a3[1],w[-2].im,w[-2].re);
1015     a += 2;
1016     a1 += 2;
1017     a2 += 2;
1018     a3 += 2;
1019     w -= 2;
1020   } while (k -= 2);
1021 }
1022
1023
1024
1025 static void u1024(register WDL_FFT_COMPLEX *a)
1026 {
1027   u512(a);
1028   u256(a + 512);
1029   u256(a + 768);
1030   upassbig(a,d1024,128);
1031 }
1032
1033 static void u2048(register WDL_FFT_COMPLEX *a)
1034 {
1035   u1024(a);
1036   u512(a + 1024);
1037   u512(a + 1536);
1038   upassbig(a,d2048,256);
1039 }
1040
1041
1042 static void u4096(register WDL_FFT_COMPLEX *a)
1043 {
1044   u2048(a);
1045   u1024(a + 2048);
1046   u1024(a + 3072);
1047   upassbig(a,d4096,512);
1048 }
1049
1050 static void u8192(register WDL_FFT_COMPLEX *a)
1051 {
1052   u4096(a);
1053   u2048(a + 4096);
1054   u2048(a + 6144);
1055   upassbig(a,d8192,1024);
1056 }
1057
1058 static void u16384(register WDL_FFT_COMPLEX *a)
1059 {
1060   u8192(a);
1061   u4096(a + 8192);
1062   u4096(a + 8192 + 4096);
1063   upassbig(a,d16384,2048);
1064 }
1065
1066 static void u32768(register WDL_FFT_COMPLEX *a)
1067 {
1068   u16384(a);
1069   u8192(a + 16384);
1070   u8192(a + 16384  + 8192 );
1071   upassbig(a,d32768,4096);
1072 }
1073
1074
1075 static void __fft_gen(WDL_FFT_COMPLEX *buf, int sz, int isfull)
1076 {
1077   int x;
1078   double div=PI*0.25/(sz+1.0);
1079
1080   if (isfull) div*=2.0;
1081
1082   for (x = 0; x < sz; x ++)
1083   {
1084     buf[x].re = (WDL_FFT_REAL) cos((x+1)*div);
1085     buf[x].im = (WDL_FFT_REAL) sin((x+1)*div);
1086   }
1087 }
1088
1089 #ifndef WDL_FFT_NO_PERMUTE
1090
1091 static unsigned int fftfreq_c(unsigned int i,unsigned int n)
1092 {
1093   unsigned int m;
1094
1095   if (n <= 2) return i;
1096
1097   m = n >> 1;
1098   if (i < m) return fftfreq_c(i,m) << 1;
1099
1100   i -= m;
1101   m >>= 1;
1102   if (i < m) return (fftfreq_c(i,m) << 2) + 1;
1103   i -= m;
1104   return ((fftfreq_c(i,m) << 2) - 1) & (n - 1);
1105 }
1106
1107 static int _idxperm[2<<FFT_MAXBITLEN];
1108
1109 static void idx_perm_calc(int offs, int n)
1110 {
1111         int i, j;
1112         _idxperm[offs] = 0;
1113         for (i = 1; i < n; ++i) {
1114                 j = fftfreq_c(i, n);
1115                 _idxperm[offs+n-j] = i;
1116         }
1117 }
1118
1119 int WDL_fft_permute(int fftsize, int idx)
1120 {
1121   return _idxperm[fftsize+idx-2];
1122 }
1123 int *WDL_fft_permute_tab(int fftsize)
1124 {
1125   return _idxperm + fftsize - 2;
1126 }
1127
1128
1129 #endif
1130
1131 void WDL_fft_init()
1132 {
1133   static int ffttabinit;
1134   if (!ffttabinit)
1135   {
1136     int i, offs;
1137         ffttabinit=1;
1138
1139 #define fft_gen(x,y) __fft_gen(x,sizeof(x)/sizeof(x[0]),y)
1140     fft_gen(d16,1);
1141     fft_gen(d32,1);
1142     fft_gen(d64,1);
1143     fft_gen(d128,1);
1144     fft_gen(d256,1);
1145     fft_gen(d512,1);
1146     fft_gen(d1024,0);
1147     fft_gen(d2048,0);
1148     fft_gen(d4096,0);
1149     fft_gen(d8192,0);
1150     fft_gen(d16384,0);
1151     fft_gen(d32768,0);
1152 #undef fft_gen
1153
1154 #ifndef WDL_FFT_NO_PERMUTE
1155           offs = 0;
1156           for (i = 2; i <= 32768; i *= 2)
1157     {
1158                   idx_perm_calc(offs, i);
1159                   offs += i;
1160           }
1161 #endif
1162
1163   }
1164 }
1165
1166 void WDL_fft(WDL_FFT_COMPLEX *buf, int len, int isInverse)
1167 {
1168   switch (len)
1169   {
1170     case 2: c2(buf); break;
1171 #define TMP(x) case x: if (!isInverse) c##x(buf); else u##x(buf); break;
1172     TMP(4)
1173     TMP(8)
1174     TMP(16)
1175     TMP(32)
1176     TMP(64)
1177     TMP(128)
1178     TMP(256)
1179     TMP(512)
1180     TMP(1024)
1181     TMP(2048)
1182     TMP(4096)
1183     TMP(8192)
1184     TMP(16384)
1185     TMP(32768)
1186 #undef TMP
1187   }
1188 }
1189
1190
1191 #if 0
1192 /* n multiple of 4, n >= 8 */
1193 void WDL_fft_realmul(WDL_FFT_REAL *a,WDL_FFT_REAL *b,int n)
1194 {
1195   if (n<8 || (n&3)) return;
1196   mulr4(a,b);
1197   WDL_fft_complexmul((WDL_FFT_COMPLEX *)(a + 4),(WDL_FFT_COMPLEX *)(b + 4),(n - 4) / 2);
1198 }
1199
1200
1201
1202
1203 //////////// begin WDL_FFT_REAL modes
1204 static void r2(register WDL_FFT_REAL *a)
1205 {
1206   register WDL_FFT_REAL t1, t2;
1207
1208   t1 = a[0] + a[1];
1209   t2 = a[0] - a[1];
1210   a[0] = t1;
1211   a[1] = t2;
1212 }
1213
1214 static void r4(register WDL_FFT_REAL *a)
1215 {
1216   register WDL_FFT_REAL t1, t2, t3, t4, t6;
1217
1218   t3 = a[0] + a[1];
1219   t4 = a[2] + a[3];
1220   t1 = a[0] - a[1];
1221   t2 = a[2] - a[3];
1222   t6 = t3 - t4;
1223   t3 += t4;
1224   a[2] = t1;
1225   a[3] = t2;
1226   a[0] = t3;
1227   a[1] = t6;
1228 }
1229
1230 static void r8(register WDL_FFT_REAL *a)
1231 {
1232   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8;
1233
1234   t2 = a[0] + a[1];
1235   t8 = a[4] + a[5];
1236   t3 = a[2] - a[3];
1237   t6 = t2 - t8;
1238   t2 += t8;
1239   t1 = a[2] + a[3];
1240   t7 = a[6] + a[7];
1241   a[2] = t6;
1242   t5 = t1 - t7;
1243   t1 += t7;
1244   t4 = a[0] - a[1];
1245   a[3] = t5;
1246   t8 = t2 - t1;
1247   t2 += t1;
1248   t7 = a[6] - a[7];
1249   a[1] = t8;
1250   t6 = t3 - t7;
1251   t3 += t7;
1252   a[0] = t2;
1253   t6 *= sqrthalf;
1254   t8 = a[4] - a[5];
1255   t3 *= sqrthalf;
1256   t1 = t4 - t6;
1257   t4 += t6;
1258   t2 = t8 - t3;
1259   t8 += t3;
1260   a[6] = t1;
1261   a[4] = t4;
1262   a[7] = t2;
1263   a[5] = t8;
1264 }
1265
1266 /* a[0...8n-1], w[0...2n-1]; n even, n >= 4 */
1267 static void rpass(register WDL_FFT_REAL *a,register const WDL_FFT_COMPLEX *w,register unsigned int n)
1268 {
1269   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6;
1270   register WDL_FFT_REAL *b;
1271   register unsigned int k;
1272
1273   b = a + 4 * n;
1274   k = n - 2;
1275
1276   RZERO(a[0],a[1],b[0],b[1]);
1277   R(a[2],a[3],b[2],b[3],w[0].re,w[0].im);
1278   R(a[4],a[5],b[4],b[5],w[1].re,w[1].im);
1279   R(a[6],a[7],b[6],b[7],w[2].re,w[2].im);
1280
1281   for (;;) {
1282     R(a[8],a[9],b[8],b[9],w[3].re,w[3].im);
1283     R(a[10],a[11],b[10],b[11],w[4].re,w[4].im);
1284     R(a[12],a[13],b[12],b[13],w[5].re,w[5].im);
1285     R(a[14],a[15],b[14],b[15],w[6].re,w[6].im);
1286     if (!(k -= 2)) break;
1287     a += 8;
1288     b += 8;
1289     w += 4;
1290   }
1291 }
1292
1293 static void r16(register WDL_FFT_REAL *a)
1294 {
1295   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6;
1296
1297   RZERO(a[0],a[1],a[8],a[9]);
1298   R(a[2],a[3],a[10],a[11],d16[0].re,d16[0].im);
1299   R(a[4],a[5],a[12],a[13],d16[1].re,d16[1].im);
1300   R(a[6],a[7],a[14],a[15],d16[2].re,d16[2].im);
1301   r8(a);
1302   c4((WDL_FFT_COMPLEX *)(a + 8));
1303 }
1304
1305 static void r32(register WDL_FFT_REAL *a)
1306 {
1307   rpass(a,d32,4);
1308   r16(a);
1309   c8((WDL_FFT_COMPLEX *)(a + 16));
1310 }
1311
1312 static void r64(register WDL_FFT_REAL *a)
1313 {
1314   rpass(a,d64,8);
1315   r32(a);
1316   c16((WDL_FFT_COMPLEX *)(a + 32));
1317 }
1318
1319 static void r128(register WDL_FFT_REAL *a)
1320 {
1321   rpass(a,d128,16);
1322   r64(a);
1323   c32((WDL_FFT_COMPLEX *)(a + 64));
1324 }
1325
1326 static void r256(register WDL_FFT_REAL *a)
1327 {
1328   rpass(a,d256,32);
1329   r128(a);
1330   c64((WDL_FFT_COMPLEX *)(a + 128));
1331 }
1332
1333 static void r512(register WDL_FFT_REAL *a)
1334 {
1335   rpass(a,d512,64);
1336   r256(a);
1337   c128((WDL_FFT_COMPLEX *)(a + 256));
1338 }
1339
1340
1341 /* a[0...8n-1], w[0...n-1]; n even, n >= 8 */
1342 static void rpassbig(register WDL_FFT_REAL *a,register const WDL_FFT_COMPLEX *w,register unsigned int n)
1343 {
1344   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6;
1345   register WDL_FFT_REAL *b;
1346   register unsigned int k;
1347
1348   b = a + 4 * n;
1349
1350   RZERO(a[0],a[1],b[0],b[1]);
1351   R(a[2],a[3],b[2],b[3],w[0].re,w[0].im);
1352
1353   k = n - 2;
1354   do {
1355     R(a[4],a[5],b[4],b[5],w[1].re,w[1].im);
1356     R(a[6],a[7],b[6],b[7],w[2].re,w[2].im);
1357     a += 4;
1358     b += 4;
1359     w += 2;
1360   } while (k -= 2);
1361
1362   RHALF(a[4],a[5],b[4],b[5]);
1363   R(a[6],a[7],b[6],b[7],w[0].im,w[0].re);
1364
1365   k = n - 2;
1366   do {
1367     R(a[8],a[9],b[8],b[9],w[-1].im,w[-1].re);
1368     R(a[10],a[11],b[10],b[11],w[-2].im,w[-2].re);
1369     a += 4;
1370     b += 4;
1371     w -= 2;
1372   } while (k -= 2);
1373 }
1374
1375
1376 static void r1024(register WDL_FFT_REAL *a)
1377 {
1378   rpassbig(a,d1024,128);
1379   r512(a);
1380   c256((WDL_FFT_COMPLEX *)(a + 512));
1381 }
1382
1383 static void r2048(register WDL_FFT_REAL *a)
1384 {
1385   rpassbig(a,d2048,256);
1386   r1024(a);
1387   c512((WDL_FFT_COMPLEX *)(a + 1024));
1388 }
1389
1390
1391 static void r4096(register WDL_FFT_REAL *a)
1392 {
1393   rpassbig(a,d4096,512);
1394   r2048(a);
1395   c1024((WDL_FFT_COMPLEX *)(a + 2048));
1396 }
1397
1398 static void r8192(register WDL_FFT_REAL *a)
1399 {
1400   rpassbig(a,d8192,1024);
1401   r4096(a);
1402   c2048((WDL_FFT_COMPLEX *)(a + 4096));
1403 }
1404
1405 static void r16384(register WDL_FFT_REAL *a)
1406 {
1407   rpassbig(a,d16384,2048);
1408   r8192(a);
1409   c4096((WDL_FFT_COMPLEX *)(a + 8192));
1410 }
1411
1412 static void r32768(register WDL_FFT_REAL *a)
1413 {
1414   rpassbig(a,d32768,4096);
1415   r16384(a);
1416   c8192((WDL_FFT_COMPLEX *)(a + 16384));
1417 }
1418
1419
1420
1421 static void v4(register WDL_FFT_REAL *a)
1422 {
1423   register WDL_FFT_REAL t1, t3, t5, t6;
1424
1425   t5 = a[0] + a[1];
1426   t6 = a[0] - a[1];
1427   t1 = t5 + a[2];
1428   t5 -= a[2];
1429   t3 = t6 + a[3];
1430   t6 -= a[3];
1431   a[0] = t1;
1432   a[1] = t5;
1433   a[2] = t3;
1434   a[3] = t6;
1435 }
1436
1437 static void v8(register WDL_FFT_REAL *a)
1438 {
1439   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8;
1440
1441   t5 = a[0] + a[1];
1442   t2 = a[4] + a[6];
1443   t8 = t5 + a[2];
1444   t5 -= a[2];
1445   t1 = a[0] - a[1];
1446   t7 = t8 + t2;
1447   t8 -= t2;
1448   t3 = a[4] - a[6];
1449   a[0] = t7;
1450   t6 = a[5] + a[7];
1451   a[1] = t8;
1452   t7 = t5 + t6;
1453   t5 -= t6;
1454   t4 = a[5] - a[7];
1455   a[4] = t7;
1456   t6 = t4 - t3;
1457   t3 += t4;
1458   a[5] = t5;
1459   t3 *= sqrthalf;
1460   t2 = t1 + a[3];
1461   t1 -= a[3];
1462   t6 *= sqrthalf;
1463   t7 = t2 - t3;
1464   t3 += t2;
1465   t8 = t1 - t6;
1466   t6 += t1;
1467   a[3] = t7;
1468   a[7] = t8;
1469   a[2] = t3;
1470   a[6] = t6;
1471 }
1472
1473 /* a[0...8n-1], w[0...2n-1]; n even, n >= 4 */
1474 static void vpass(register WDL_FFT_REAL *a,register const WDL_FFT_COMPLEX *w,register unsigned int n)
1475 {
1476   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6;
1477   register WDL_FFT_REAL *b;
1478   register unsigned int k;
1479
1480   b = a + 4 * n;
1481   k = n - 2;
1482
1483   VZERO(a[0],a[1],b[0],b[1]);
1484   V(a[2],a[3],b[2],b[3],w[0].re,w[0].im);
1485   V(a[4],a[5],b[4],b[5],w[1].re,w[1].im);
1486   V(a[6],a[7],b[6],b[7],w[2].re,w[2].im);
1487
1488   for (;;) {
1489     V(a[8],a[9],b[8],b[9],w[3].re,w[3].im);
1490     V(a[10],a[11],b[10],b[11],w[4].re,w[4].im);
1491     V(a[12],a[13],b[12],b[13],w[5].re,w[5].im);
1492     V(a[14],a[15],b[14],b[15],w[6].re,w[6].im);
1493     if (!(k -= 2)) break;
1494     a += 8;
1495     b += 8;
1496     w += 4;
1497   }
1498 }
1499
1500 static void v16(register WDL_FFT_REAL *a)
1501 {
1502   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6;
1503
1504   u4((WDL_FFT_COMPLEX *)(a + 8));
1505   v8(a);
1506   VZERO(a[0],a[1],a[8],a[9]);
1507   V(a[2],a[3],a[10],a[11],d16[0].re,d16[0].im);
1508   V(a[4],a[5],a[12],a[13],d16[1].re,d16[1].im);
1509   V(a[6],a[7],a[14],a[15],d16[2].re,d16[2].im);
1510 }
1511
1512 static void v32(register WDL_FFT_REAL *a)
1513 {
1514   u8((WDL_FFT_COMPLEX *)(a + 16));
1515   v16(a);
1516   vpass(a,d32,4);
1517 }
1518
1519 static void v64(register WDL_FFT_REAL *a)
1520 {
1521   u16((WDL_FFT_COMPLEX *)(a + 32));
1522   v32(a);
1523   vpass(a,d64,8);
1524 }
1525
1526 static void v128(register WDL_FFT_REAL *a)
1527 {
1528   u32((WDL_FFT_COMPLEX *)(a + 64));
1529   v64(a);
1530   vpass(a,d128,16);
1531 }
1532
1533 static void v256(register WDL_FFT_REAL *a)
1534 {
1535   u64((WDL_FFT_COMPLEX *)(a + 128));
1536   v128(a);
1537   vpass(a,d256,32);
1538 }
1539
1540 static void v512(register WDL_FFT_REAL *a)
1541 {
1542   u128((WDL_FFT_COMPLEX *)(a + 256));
1543   v256(a);
1544   vpass(a,d512,64);
1545 }
1546
1547
1548
1549 /* a[0...8n-1], w[0...n-1]; n even, n >= 8 */
1550 static void vpassbig(register WDL_FFT_REAL *a,register const WDL_FFT_COMPLEX *w,register unsigned int n)
1551 {
1552   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6;
1553   register WDL_FFT_REAL *b;
1554   register unsigned int k;
1555
1556   b = a + 4 * n;
1557
1558   VZERO(a[0],a[1],b[0],b[1]);
1559   V(a[2],a[3],b[2],b[3],w[0].re,w[0].im);
1560
1561   k = n - 2;
1562   do {
1563     V(a[4],a[5],b[4],b[5],w[1].re,w[1].im);
1564     V(a[6],a[7],b[6],b[7],w[2].re,w[2].im);
1565     a += 4;
1566     b += 4;
1567     w += 2;
1568   } while (k -= 2);
1569
1570   VHALF(a[4],a[5],b[4],b[5]);
1571   V(a[6],a[7],b[6],b[7],w[0].im,w[0].re);
1572
1573   k = n - 2;
1574   do {
1575     V(a[8],a[9],b[8],b[9],w[-1].im,w[-1].re);
1576     V(a[10],a[11],b[10],b[11],w[-2].im,w[-2].re);
1577     a += 4;
1578     b += 4;
1579     w -= 2;
1580   } while (k -= 2);
1581 }
1582
1583
1584 static void v1024(register WDL_FFT_REAL *a)
1585 {
1586   u256((WDL_FFT_COMPLEX *)(a + 512));
1587   v512(a);
1588   vpassbig(a,d1024,128);
1589 }
1590
1591 static void v2048(register WDL_FFT_REAL *a)
1592 {
1593   u512((WDL_FFT_COMPLEX *)(a + 1024));
1594   v1024(a);
1595   vpassbig(a,d2048,256);
1596 }
1597
1598
1599 static void v4096(register WDL_FFT_REAL *a)
1600 {
1601   u1024((WDL_FFT_COMPLEX *)(a + 2048));
1602   v2048(a);
1603   vpassbig(a,d4096,512);
1604 }
1605
1606 static void v8192(register WDL_FFT_REAL *a)
1607 {
1608   u2048((WDL_FFT_COMPLEX *)(a + 4096));
1609   v4096(a);
1610   vpassbig(a,d8192,1024);
1611 }
1612
1613 static void v16384(register WDL_FFT_REAL *a)
1614 {
1615   u4096((WDL_FFT_COMPLEX *)(a + 8192));
1616   v8192(a);
1617   vpassbig(a,d16384,2048);
1618 }
1619
1620 static void v32768(register WDL_FFT_REAL *a)
1621 {
1622   u8192((WDL_FFT_COMPLEX *)(a + 16384));
1623   v16384(a);
1624   vpassbig(a,d32768,4096);
1625 }
1626
1627
1628 void WDL_real_fft(WDL_FFT_REAL *buf, int len, int isInverse)
1629 {
1630   switch (len)
1631   {
1632     case 2: r2(buf); break;
1633 #define TMP(x) case x: if (!isInverse) r##x(buf); else v##x(buf); break;
1634     TMP(4)
1635     TMP(8)
1636     TMP(16)
1637     TMP(32)
1638     TMP(64)
1639     TMP(128)
1640     TMP(256)
1641     TMP(512)
1642     TMP(1024)
1643     TMP(2048)
1644     TMP(4096)
1645     TMP(8192)
1646     TMP(16384)
1647     TMP(32768)
1648 #undef TMP
1649   }
1650 }
1651
1652 #endif