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-16];
1122 }
1123
1124 #endif
1125
1126 void WDL_fft_init()
1127 {
1128   static int ffttabinit;
1129   if (!ffttabinit)
1130   {
1131     int i, offs;
1132         ffttabinit=1;
1133
1134 #define fft_gen(x,y) __fft_gen(x,sizeof(x)/sizeof(x[0]),y)
1135     fft_gen(d16,1);
1136     fft_gen(d32,1);
1137     fft_gen(d64,1);
1138     fft_gen(d128,1);
1139     fft_gen(d256,1);
1140     fft_gen(d512,1);
1141     fft_gen(d1024,0);
1142     fft_gen(d2048,0);
1143     fft_gen(d4096,0);
1144     fft_gen(d8192,0);
1145     fft_gen(d16384,0);
1146     fft_gen(d32768,0);
1147 #undef fft_gen
1148
1149 #ifndef WDL_FFT_NO_PERMUTE
1150           offs = 0;
1151           for (i = 16; i <= 32768; i *= 2)
1152     {
1153                   idx_perm_calc(offs, i);
1154                   offs += i;
1155           }
1156 #endif
1157
1158   }
1159 }
1160
1161 void WDL_fft(WDL_FFT_COMPLEX *buf, int len, int isInverse)
1162 {
1163   switch (len)
1164   {
1165     case 2: c2(buf); break;
1166 #define TMP(x) case x: if (!isInverse) c##x(buf); else u##x(buf); break;
1167     TMP(4)
1168     TMP(8)
1169     TMP(16)
1170     TMP(32)
1171     TMP(64)
1172     TMP(128)
1173     TMP(256)
1174     TMP(512)
1175     TMP(1024)
1176     TMP(2048)
1177     TMP(4096)
1178     TMP(8192)
1179     TMP(16384)
1180     TMP(32768)
1181 #undef TMP
1182   }
1183 }
1184
1185
1186 #if 0
1187 /* n multiple of 4, n >= 8 */
1188 void WDL_fft_realmul(WDL_FFT_REAL *a,WDL_FFT_REAL *b,int n)
1189 {
1190   if (n<8 || (n&3)) return;
1191   mulr4(a,b);
1192   WDL_fft_complexmul((WDL_FFT_COMPLEX *)(a + 4),(WDL_FFT_COMPLEX *)(b + 4),(n - 4) / 2);
1193 }
1194
1195
1196
1197
1198 //////////// begin WDL_FFT_REAL modes
1199 static void r2(register WDL_FFT_REAL *a)
1200 {
1201   register WDL_FFT_REAL t1, t2;
1202
1203   t1 = a[0] + a[1];
1204   t2 = a[0] - a[1];
1205   a[0] = t1;
1206   a[1] = t2;
1207 }
1208
1209 static void r4(register WDL_FFT_REAL *a)
1210 {
1211   register WDL_FFT_REAL t1, t2, t3, t4, t6;
1212
1213   t3 = a[0] + a[1];
1214   t4 = a[2] + a[3];
1215   t1 = a[0] - a[1];
1216   t2 = a[2] - a[3];
1217   t6 = t3 - t4;
1218   t3 += t4;
1219   a[2] = t1;
1220   a[3] = t2;
1221   a[0] = t3;
1222   a[1] = t6;
1223 }
1224
1225 static void r8(register WDL_FFT_REAL *a)
1226 {
1227   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8;
1228
1229   t2 = a[0] + a[1];
1230   t8 = a[4] + a[5];
1231   t3 = a[2] - a[3];
1232   t6 = t2 - t8;
1233   t2 += t8;
1234   t1 = a[2] + a[3];
1235   t7 = a[6] + a[7];
1236   a[2] = t6;
1237   t5 = t1 - t7;
1238   t1 += t7;
1239   t4 = a[0] - a[1];
1240   a[3] = t5;
1241   t8 = t2 - t1;
1242   t2 += t1;
1243   t7 = a[6] - a[7];
1244   a[1] = t8;
1245   t6 = t3 - t7;
1246   t3 += t7;
1247   a[0] = t2;
1248   t6 *= sqrthalf;
1249   t8 = a[4] - a[5];
1250   t3 *= sqrthalf;
1251   t1 = t4 - t6;
1252   t4 += t6;
1253   t2 = t8 - t3;
1254   t8 += t3;
1255   a[6] = t1;
1256   a[4] = t4;
1257   a[7] = t2;
1258   a[5] = t8;
1259 }
1260
1261 /* a[0...8n-1], w[0...2n-1]; n even, n >= 4 */
1262 static void rpass(register WDL_FFT_REAL *a,register const WDL_FFT_COMPLEX *w,register unsigned int n)
1263 {
1264   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6;
1265   register WDL_FFT_REAL *b;
1266   register unsigned int k;
1267
1268   b = a + 4 * n;
1269   k = n - 2;
1270
1271   RZERO(a[0],a[1],b[0],b[1]);
1272   R(a[2],a[3],b[2],b[3],w[0].re,w[0].im);
1273   R(a[4],a[5],b[4],b[5],w[1].re,w[1].im);
1274   R(a[6],a[7],b[6],b[7],w[2].re,w[2].im);
1275
1276   for (;;) {
1277     R(a[8],a[9],b[8],b[9],w[3].re,w[3].im);
1278     R(a[10],a[11],b[10],b[11],w[4].re,w[4].im);
1279     R(a[12],a[13],b[12],b[13],w[5].re,w[5].im);
1280     R(a[14],a[15],b[14],b[15],w[6].re,w[6].im);
1281     if (!(k -= 2)) break;
1282     a += 8;
1283     b += 8;
1284     w += 4;
1285   }
1286 }
1287
1288 static void r16(register WDL_FFT_REAL *a)
1289 {
1290   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6;
1291
1292   RZERO(a[0],a[1],a[8],a[9]);
1293   R(a[2],a[3],a[10],a[11],d16[0].re,d16[0].im);
1294   R(a[4],a[5],a[12],a[13],d16[1].re,d16[1].im);
1295   R(a[6],a[7],a[14],a[15],d16[2].re,d16[2].im);
1296   r8(a);
1297   c4((WDL_FFT_COMPLEX *)(a + 8));
1298 }
1299
1300 static void r32(register WDL_FFT_REAL *a)
1301 {
1302   rpass(a,d32,4);
1303   r16(a);
1304   c8((WDL_FFT_COMPLEX *)(a + 16));
1305 }
1306
1307 static void r64(register WDL_FFT_REAL *a)
1308 {
1309   rpass(a,d64,8);
1310   r32(a);
1311   c16((WDL_FFT_COMPLEX *)(a + 32));
1312 }
1313
1314 static void r128(register WDL_FFT_REAL *a)
1315 {
1316   rpass(a,d128,16);
1317   r64(a);
1318   c32((WDL_FFT_COMPLEX *)(a + 64));
1319 }
1320
1321 static void r256(register WDL_FFT_REAL *a)
1322 {
1323   rpass(a,d256,32);
1324   r128(a);
1325   c64((WDL_FFT_COMPLEX *)(a + 128));
1326 }
1327
1328 static void r512(register WDL_FFT_REAL *a)
1329 {
1330   rpass(a,d512,64);
1331   r256(a);
1332   c128((WDL_FFT_COMPLEX *)(a + 256));
1333 }
1334
1335
1336 /* a[0...8n-1], w[0...n-1]; n even, n >= 8 */
1337 static void rpassbig(register WDL_FFT_REAL *a,register const WDL_FFT_COMPLEX *w,register unsigned int n)
1338 {
1339   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6;
1340   register WDL_FFT_REAL *b;
1341   register unsigned int k;
1342
1343   b = a + 4 * n;
1344
1345   RZERO(a[0],a[1],b[0],b[1]);
1346   R(a[2],a[3],b[2],b[3],w[0].re,w[0].im);
1347
1348   k = n - 2;
1349   do {
1350     R(a[4],a[5],b[4],b[5],w[1].re,w[1].im);
1351     R(a[6],a[7],b[6],b[7],w[2].re,w[2].im);
1352     a += 4;
1353     b += 4;
1354     w += 2;
1355   } while (k -= 2);
1356
1357   RHALF(a[4],a[5],b[4],b[5]);
1358   R(a[6],a[7],b[6],b[7],w[0].im,w[0].re);
1359
1360   k = n - 2;
1361   do {
1362     R(a[8],a[9],b[8],b[9],w[-1].im,w[-1].re);
1363     R(a[10],a[11],b[10],b[11],w[-2].im,w[-2].re);
1364     a += 4;
1365     b += 4;
1366     w -= 2;
1367   } while (k -= 2);
1368 }
1369
1370
1371 static void r1024(register WDL_FFT_REAL *a)
1372 {
1373   rpassbig(a,d1024,128);
1374   r512(a);
1375   c256((WDL_FFT_COMPLEX *)(a + 512));
1376 }
1377
1378 static void r2048(register WDL_FFT_REAL *a)
1379 {
1380   rpassbig(a,d2048,256);
1381   r1024(a);
1382   c512((WDL_FFT_COMPLEX *)(a + 1024));
1383 }
1384
1385
1386 static void r4096(register WDL_FFT_REAL *a)
1387 {
1388   rpassbig(a,d4096,512);
1389   r2048(a);
1390   c1024((WDL_FFT_COMPLEX *)(a + 2048));
1391 }
1392
1393 static void r8192(register WDL_FFT_REAL *a)
1394 {
1395   rpassbig(a,d8192,1024);
1396   r4096(a);
1397   c2048((WDL_FFT_COMPLEX *)(a + 4096));
1398 }
1399
1400 static void r16384(register WDL_FFT_REAL *a)
1401 {
1402   rpassbig(a,d16384,2048);
1403   r8192(a);
1404   c4096((WDL_FFT_COMPLEX *)(a + 8192));
1405 }
1406
1407 static void r32768(register WDL_FFT_REAL *a)
1408 {
1409   rpassbig(a,d32768,4096);
1410   r16384(a);
1411   c8192((WDL_FFT_COMPLEX *)(a + 16384));
1412 }
1413
1414
1415
1416 static void v4(register WDL_FFT_REAL *a)
1417 {
1418   register WDL_FFT_REAL t1, t3, t5, t6;
1419
1420   t5 = a[0] + a[1];
1421   t6 = a[0] - a[1];
1422   t1 = t5 + a[2];
1423   t5 -= a[2];
1424   t3 = t6 + a[3];
1425   t6 -= a[3];
1426   a[0] = t1;
1427   a[1] = t5;
1428   a[2] = t3;
1429   a[3] = t6;
1430 }
1431
1432 static void v8(register WDL_FFT_REAL *a)
1433 {
1434   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8;
1435
1436   t5 = a[0] + a[1];
1437   t2 = a[4] + a[6];
1438   t8 = t5 + a[2];
1439   t5 -= a[2];
1440   t1 = a[0] - a[1];
1441   t7 = t8 + t2;
1442   t8 -= t2;
1443   t3 = a[4] - a[6];
1444   a[0] = t7;
1445   t6 = a[5] + a[7];
1446   a[1] = t8;
1447   t7 = t5 + t6;
1448   t5 -= t6;
1449   t4 = a[5] - a[7];
1450   a[4] = t7;
1451   t6 = t4 - t3;
1452   t3 += t4;
1453   a[5] = t5;
1454   t3 *= sqrthalf;
1455   t2 = t1 + a[3];
1456   t1 -= a[3];
1457   t6 *= sqrthalf;
1458   t7 = t2 - t3;
1459   t3 += t2;
1460   t8 = t1 - t6;
1461   t6 += t1;
1462   a[3] = t7;
1463   a[7] = t8;
1464   a[2] = t3;
1465   a[6] = t6;
1466 }
1467
1468 /* a[0...8n-1], w[0...2n-1]; n even, n >= 4 */
1469 static void vpass(register WDL_FFT_REAL *a,register const WDL_FFT_COMPLEX *w,register unsigned int n)
1470 {
1471   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6;
1472   register WDL_FFT_REAL *b;
1473   register unsigned int k;
1474
1475   b = a + 4 * n;
1476   k = n - 2;
1477
1478   VZERO(a[0],a[1],b[0],b[1]);
1479   V(a[2],a[3],b[2],b[3],w[0].re,w[0].im);
1480   V(a[4],a[5],b[4],b[5],w[1].re,w[1].im);
1481   V(a[6],a[7],b[6],b[7],w[2].re,w[2].im);
1482
1483   for (;;) {
1484     V(a[8],a[9],b[8],b[9],w[3].re,w[3].im);
1485     V(a[10],a[11],b[10],b[11],w[4].re,w[4].im);
1486     V(a[12],a[13],b[12],b[13],w[5].re,w[5].im);
1487     V(a[14],a[15],b[14],b[15],w[6].re,w[6].im);
1488     if (!(k -= 2)) break;
1489     a += 8;
1490     b += 8;
1491     w += 4;
1492   }
1493 }
1494
1495 static void v16(register WDL_FFT_REAL *a)
1496 {
1497   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6;
1498
1499   u4((WDL_FFT_COMPLEX *)(a + 8));
1500   v8(a);
1501   VZERO(a[0],a[1],a[8],a[9]);
1502   V(a[2],a[3],a[10],a[11],d16[0].re,d16[0].im);
1503   V(a[4],a[5],a[12],a[13],d16[1].re,d16[1].im);
1504   V(a[6],a[7],a[14],a[15],d16[2].re,d16[2].im);
1505 }
1506
1507 static void v32(register WDL_FFT_REAL *a)
1508 {
1509   u8((WDL_FFT_COMPLEX *)(a + 16));
1510   v16(a);
1511   vpass(a,d32,4);
1512 }
1513
1514 static void v64(register WDL_FFT_REAL *a)
1515 {
1516   u16((WDL_FFT_COMPLEX *)(a + 32));
1517   v32(a);
1518   vpass(a,d64,8);
1519 }
1520
1521 static void v128(register WDL_FFT_REAL *a)
1522 {
1523   u32((WDL_FFT_COMPLEX *)(a + 64));
1524   v64(a);
1525   vpass(a,d128,16);
1526 }
1527
1528 static void v256(register WDL_FFT_REAL *a)
1529 {
1530   u64((WDL_FFT_COMPLEX *)(a + 128));
1531   v128(a);
1532   vpass(a,d256,32);
1533 }
1534
1535 static void v512(register WDL_FFT_REAL *a)
1536 {
1537   u128((WDL_FFT_COMPLEX *)(a + 256));
1538   v256(a);
1539   vpass(a,d512,64);
1540 }
1541
1542
1543
1544 /* a[0...8n-1], w[0...n-1]; n even, n >= 8 */
1545 static void vpassbig(register WDL_FFT_REAL *a,register const WDL_FFT_COMPLEX *w,register unsigned int n)
1546 {
1547   register WDL_FFT_REAL t1, t2, t3, t4, t5, t6;
1548   register WDL_FFT_REAL *b;
1549   register unsigned int k;
1550
1551   b = a + 4 * n;
1552
1553   VZERO(a[0],a[1],b[0],b[1]);
1554   V(a[2],a[3],b[2],b[3],w[0].re,w[0].im);
1555
1556   k = n - 2;
1557   do {
1558     V(a[4],a[5],b[4],b[5],w[1].re,w[1].im);
1559     V(a[6],a[7],b[6],b[7],w[2].re,w[2].im);
1560     a += 4;
1561     b += 4;
1562     w += 2;
1563   } while (k -= 2);
1564
1565   VHALF(a[4],a[5],b[4],b[5]);
1566   V(a[6],a[7],b[6],b[7],w[0].im,w[0].re);
1567
1568   k = n - 2;
1569   do {
1570     V(a[8],a[9],b[8],b[9],w[-1].im,w[-1].re);
1571     V(a[10],a[11],b[10],b[11],w[-2].im,w[-2].re);
1572     a += 4;
1573     b += 4;
1574     w -= 2;
1575   } while (k -= 2);
1576 }
1577
1578
1579 static void v1024(register WDL_FFT_REAL *a)
1580 {
1581   u256((WDL_FFT_COMPLEX *)(a + 512));
1582   v512(a);
1583   vpassbig(a,d1024,128);
1584 }
1585
1586 static void v2048(register WDL_FFT_REAL *a)
1587 {
1588   u512((WDL_FFT_COMPLEX *)(a + 1024));
1589   v1024(a);
1590   vpassbig(a,d2048,256);
1591 }
1592
1593
1594 static void v4096(register WDL_FFT_REAL *a)
1595 {
1596   u1024((WDL_FFT_COMPLEX *)(a + 2048));
1597   v2048(a);
1598   vpassbig(a,d4096,512);
1599 }
1600
1601 static void v8192(register WDL_FFT_REAL *a)
1602 {
1603   u2048((WDL_FFT_COMPLEX *)(a + 4096));
1604   v4096(a);
1605   vpassbig(a,d8192,1024);
1606 }
1607
1608 static void v16384(register WDL_FFT_REAL *a)
1609 {
1610   u4096((WDL_FFT_COMPLEX *)(a + 8192));
1611   v8192(a);
1612   vpassbig(a,d16384,2048);
1613 }
1614
1615 static void v32768(register WDL_FFT_REAL *a)
1616 {
1617   u8192((WDL_FFT_COMPLEX *)(a + 16384));
1618   v16384(a);
1619   vpassbig(a,d32768,4096);
1620 }
1621
1622
1623 void WDL_real_fft(WDL_FFT_REAL *buf, int len, int isInverse)
1624 {
1625   switch (len)
1626   {
1627     case 2: r2(buf); break;
1628 #define TMP(x) case x: if (!isInverse) r##x(buf); else v##x(buf); break;
1629     TMP(4)
1630     TMP(8)
1631     TMP(16)
1632     TMP(32)
1633     TMP(64)
1634     TMP(128)
1635     TMP(256)
1636     TMP(512)
1637     TMP(1024)
1638     TMP(2048)
1639     TMP(4096)
1640     TMP(8192)
1641     TMP(16384)
1642     TMP(32768)
1643 #undef TMP
1644   }
1645 }
1646
1647 #endif