From 11ac7454dbc1b905d0d89880b93554a18c126ae4 Mon Sep 17 00:00:00 2001 From: Vas Crabb Date: Wed, 4 Mar 2020 22:41:15 +1100 Subject: Maintenance: * Removed empty nl_examples from dist.mak * Added copyright acknowledgements and full text of licenses to binary distribution * Fixed up the list of third-party libraries * Moved WDL fft.c to 3rdparty --- 3rdparty/README.md | 12 +- 3rdparty/wdlfft/fft.c | 1199 +++++++++++++++++++++++++++++++++++++++++++++++++ 3rdparty/wdlfft/fft.h | 77 ++++ 3 files changed, 1283 insertions(+), 5 deletions(-) create mode 100644 3rdparty/wdlfft/fft.c create mode 100644 3rdparty/wdlfft/fft.h (limited to '3rdparty') diff --git a/3rdparty/README.md b/3rdparty/README.md index f60fb7c168b..00616b7773d 100644 --- a/3rdparty/README.md +++ b/3rdparty/README.md @@ -24,11 +24,11 @@ genie - [The BSD 3-Clause License](http://opensource.org/licenses/BSD-3-Clause) glm - [The MIT License (MIT)](http://opensource.org/licenses/MIT) -libflac - [The BSD 3-Clause License](http://opensource.org/licenses/BSD-3-Clause) +libflac - [The BSD 3-Clause License](http://opensource.org/licenses/BSD-3-Clause) (codecs), [The GNU General Public License, version 2](https://opensource.org/licenses/GPL-2.0) libjpeg - [Custom BSD-like](https://github.com/numenta/nupic/blob/master/external/licenses/LICENSE.libjpeg-6b.txt) -linenoise-ng - [The BSD 3-Clause License](http://opensource.org/licenses/BSD-3-Clause) +linenoise - [The BSD 2-Clause License](http://opensource.org/licenses/BSD-2-Clause) lsqlite3 - [The MIT License (MIT)](http://opensource.org/licenses/MIT) @@ -40,7 +40,7 @@ lua-zlib - [The MIT License (MIT)](http://opensource.org/licenses/MIT) luafilesystem - [The MIT License (MIT)](http://opensource.org/licenses/MIT) -lzma - [The GNU Lesser General Public License](http://opensource.org/licenses/LGPL-2.1) +lzma - Custom public domain license minimp3 - [Creative Commons Zero v1.0 Universal](https://creativecommons.org/publicdomain/zero/1.0/) @@ -52,7 +52,7 @@ portmidi - [The MIT License (MIT)](http://opensource.org/licenses/MIT) pugixml - [The MIT License (MIT)](http://opensource.org/licenses/MIT) -rapidjson - [The BSD 3-Clause License](http://opensource.org/licenses/BSD-3-Clause) +rapidjson - [The MIT License (MIT)](http://opensource.org/licenses/MIT), [The BSD 3-Clause License](http://opensource.org/licenses/BSD-3-Clause) (msinttypes), JSON license (JSON_checker) sdl2 - [zlib license](http://opensource.org/licenses/Zlib) @@ -62,12 +62,14 @@ softfloat3 - [The BSD 3-Clause License](http://opensource.org/licenses/BSD-3-Cla sol2 - [The MIT License (MIT)](http://opensource.org/licenses/MIT) -sqlite3 - Public Domain +sqlite3 - Custom public domain license tap-windows6.h - Dual-licensed [The GNU General Public License, version 2](https://opensource.org/licenses/GPL-2.0) and [The MIT License (MIT)](http://opensource.org/licenses/MIT) utf8proc - [The MIT License (MIT)](http://opensource.org/licenses/MIT) +wdlfft - [zlib license](http://opensource.org/licenses/Zlib) + winpcap - [The BSD 3-Clause License](http://opensource.org/licenses/BSD-3-Clause) zlib - [zlib license](http://opensource.org/licenses/Zlib) diff --git a/3rdparty/wdlfft/fft.c b/3rdparty/wdlfft/fft.c new file mode 100644 index 00000000000..8b15d6ca35b --- /dev/null +++ b/3rdparty/wdlfft/fft.c @@ -0,0 +1,1199 @@ +/* + WDL - fft.cpp + Copyright (C) 2006 and later Cockos Incorporated + Copyright 1999 D. J. Bernstein + + This software is provided 'as-is', without any express or implied + warranty. In no event will the authors be held liable for any damages + arising from the use of this software. + + Permission is granted to anyone to use this software for any purpose, + including commercial applications, and to alter it and redistribute it + freely, subject to the following restrictions: + + 1. The origin of this software must not be misrepresented; you must not + claim that you wrote the original software. If you use this software + in a product, an acknowledgment in the product documentation would be + appreciated but is not required. + 2. Altered source versions must be plainly marked as such, and must not be + misrepresented as being the original software. + 3. This notice may not be removed or altered from any source distribution. + + + + This file implements the WDL FFT library. These routines are based on the + DJBFFT library, which are Copyright 1999 D. J. Bernstein, djb@pobox.com + + The DJB FFT web page is: http://cr.yp.to/djbfft.html + +*/ + + +// this is based on djbfft + +#include +#include "fft.h" + + +#define FFT_MAXBITLEN 15 + +#ifdef _MSC_VER +#define inline __inline +#endif + +#define PI 3.1415926535897932384626433832795 + +static WDL_FFT_COMPLEX d16[3]; +static WDL_FFT_COMPLEX d32[7]; +static WDL_FFT_COMPLEX d64[15]; +static WDL_FFT_COMPLEX d128[31]; +static WDL_FFT_COMPLEX d256[63]; +static WDL_FFT_COMPLEX d512[127]; +static WDL_FFT_COMPLEX d1024[127]; +static WDL_FFT_COMPLEX d2048[255]; +static WDL_FFT_COMPLEX d4096[511]; +static WDL_FFT_COMPLEX d8192[1023]; +static WDL_FFT_COMPLEX d16384[2047]; +static WDL_FFT_COMPLEX d32768[4095]; + + +#define sqrthalf (d16[1].re) + +#define VOL *(volatile WDL_FFT_REAL *)& + +#define TRANSFORM(a0,a1,a2,a3,wre,wim) { \ + t6 = a2.re; \ + t1 = a0.re - t6; \ + t6 += a0.re; \ + a0.re = t6; \ + t3 = a3.im; \ + t4 = a1.im - t3; \ + t8 = t1 - t4; \ + t1 += t4; \ + t3 += a1.im; \ + a1.im = t3; \ + t5 = wre; \ + t7 = t8 * t5; \ + t4 = t1 * t5; \ + t8 *= wim; \ + t2 = a3.re; \ + t3 = a1.re - t2; \ + t2 += a1.re; \ + a1.re = t2; \ + t1 *= wim; \ + t6 = a2.im; \ + t2 = a0.im - t6; \ + t6 += a0.im; \ + a0.im = t6; \ + t6 = t2 + t3; \ + t2 -= t3; \ + t3 = t6 * wim; \ + t7 -= t3; \ + a2.re = t7; \ + t6 *= t5; \ + t6 += t8; \ + a2.im = t6; \ + t5 *= t2; \ + t5 -= t1; \ + a3.im = t5; \ + t2 *= wim; \ + t4 += t2; \ + a3.re = t4; \ + } + +#define TRANSFORMHALF(a0,a1,a2,a3) { \ + t1 = a2.re; \ + t5 = a0.re - t1; \ + t1 += a0.re; \ + a0.re = t1; \ + t4 = a3.im; \ + t8 = a1.im - t4; \ + t1 = t5 - t8; \ + t5 += t8; \ + t4 += a1.im; \ + a1.im = t4; \ + t3 = a3.re; \ + t7 = a1.re - t3; \ + t3 += a1.re; \ + a1.re = t3; \ + t8 = a2.im; \ + t6 = a0.im - t8; \ + t2 = t6 + t7; \ + t6 -= t7; \ + t8 += a0.im; \ + a0.im = t8; \ + t4 = t6 + t5; \ + t3 = sqrthalf; \ + t4 *= t3; \ + a3.re = t4; \ + t6 -= t5; \ + t6 *= t3; \ + a3.im = t6; \ + t7 = t1 - t2; \ + t7 *= t3; \ + a2.re = t7; \ + t2 += t1; \ + t2 *= t3; \ + a2.im = t2; \ + } + +#define TRANSFORMZERO(a0,a1,a2,a3) { \ + t5 = a2.re; \ + t1 = a0.re - t5; \ + t5 += a0.re; \ + a0.re = t5; \ + t8 = a3.im; \ + t4 = a1.im - t8; \ + t7 = a3.re; \ + t6 = t1 - t4; \ + a2.re = t6; \ + t1 += t4; \ + a3.re = t1; \ + t8 += a1.im; \ + a1.im = t8; \ + t3 = a1.re - t7; \ + t7 += a1.re; \ + a1.re = t7; \ + t6 = a2.im; \ + t2 = a0.im - t6; \ + t7 = t2 + t3; \ + a2.im = t7; \ + t2 -= t3; \ + a3.im = t2; \ + t6 += a0.im; \ + a0.im = t6; \ + } + +#define UNTRANSFORM(a0,a1,a2,a3,wre,wim) { \ + t6 = VOL wre; \ + t1 = VOL a2.re; \ + t1 *= t6; \ + t8 = VOL wim; \ + t3 = VOL a2.im; \ + t3 *= t8; \ + t2 = VOL a2.im; \ + t4 = VOL a2.re; \ + t5 = VOL a3.re; \ + t5 *= t6; \ + t7 = VOL a3.im; \ + t1 += t3; \ + t7 *= t8; \ + t5 -= t7; \ + t3 = t5 + t1; \ + t5 -= t1; \ + t2 *= t6; \ + t6 *= a3.im; \ + t4 *= t8; \ + t2 -= t4; \ + t8 *= a3.re; \ + t6 += t8; \ + t1 = a0.re - t3; \ + t3 += a0.re; \ + a0.re = t3; \ + t7 = a1.im - t5; \ + t5 += a1.im; \ + a1.im = t5; \ + t4 = t2 - t6; \ + t6 += t2; \ + t8 = a1.re - t4; \ + t4 += a1.re; \ + a1.re = t4; \ + t2 = a0.im - t6; \ + t6 += a0.im; \ + a0.im = t6; \ + a2.re = t1; \ + a3.im = t7; \ + a3.re = t8; \ + a2.im = t2; \ + } + + +#define UNTRANSFORMHALF(a0,a1,a2,a3) { \ + t6 = sqrthalf; \ + t1 = a2.re; \ + t2 = a2.im - t1; \ + t2 *= t6; \ + t1 += a2.im; \ + t1 *= t6; \ + t4 = a3.im; \ + t3 = a3.re - t4; \ + t3 *= t6; \ + t4 += a3.re; \ + t4 *= t6; \ + t8 = t3 - t1; \ + t7 = t2 - t4; \ + t1 += t3; \ + t2 += t4; \ + t4 = a1.im - t8; \ + a3.im = t4; \ + t8 += a1.im; \ + a1.im = t8; \ + t3 = a1.re - t7; \ + a3.re = t3; \ + t7 += a1.re; \ + a1.re = t7; \ + t5 = a0.re - t1; \ + a2.re = t5; \ + t1 += a0.re; \ + a0.re = t1; \ + t6 = a0.im - t2; \ + a2.im = t6; \ + t2 += a0.im; \ + a0.im = t2; \ + } + +#define UNTRANSFORMZERO(a0,a1,a2,a3) { \ + t2 = a3.im; \ + t3 = a2.im - t2; \ + t2 += a2.im; \ + t1 = a2.re; \ + t4 = a3.re - t1; \ + t1 += a3.re; \ + t5 = a0.re - t1; \ + a2.re = t5; \ + t6 = a0.im - t2; \ + a2.im = t6; \ + t7 = a1.re - t3; \ + a3.re = t7; \ + t8 = a1.im - t4; \ + a3.im = t8; \ + t1 += a0.re; \ + a0.re = t1; \ + t2 += a0.im; \ + a0.im = t2; \ + t3 += a1.re; \ + a1.re = t3; \ + t4 += a1.im; \ + a1.im = t4; \ + } + +static void c2(register WDL_FFT_COMPLEX *a) +{ + register WDL_FFT_REAL t1; + + t1 = a[1].re; + a[1].re = a[0].re - t1; + a[0].re += t1; + + t1 = a[1].im; + a[1].im = a[0].im - t1; + a[0].im += t1; +} + +static inline void c4(register WDL_FFT_COMPLEX *a) +{ + register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8; + + t5 = a[2].re; + t1 = a[0].re - t5; + t7 = a[3].re; + t5 += a[0].re; + t3 = a[1].re - t7; + t7 += a[1].re; + t8 = t5 + t7; + a[0].re = t8; + t5 -= t7; + a[1].re = t5; + t6 = a[2].im; + t2 = a[0].im - t6; + t6 += a[0].im; + t5 = a[3].im; + a[2].im = t2 + t3; + t2 -= t3; + a[3].im = t2; + t4 = a[1].im - t5; + a[3].re = t1 + t4; + t1 -= t4; + a[2].re = t1; + t5 += a[1].im; + a[0].im = t6 + t5; + t6 -= t5; + a[1].im = t6; +} + +static void c8(register WDL_FFT_COMPLEX *a) +{ + register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8; + + t7 = a[4].im; + t4 = a[0].im - t7; + t7 += a[0].im; + a[0].im = t7; + + t8 = a[6].re; + t5 = a[2].re - t8; + t8 += a[2].re; + a[2].re = t8; + + t7 = a[6].im; + a[6].im = t4 - t5; + t4 += t5; + a[4].im = t4; + + t6 = a[2].im - t7; + t7 += a[2].im; + a[2].im = t7; + + t8 = a[4].re; + t3 = a[0].re - t8; + t8 += a[0].re; + a[0].re = t8; + + a[4].re = t3 - t6; + t3 += t6; + a[6].re = t3; + + t7 = a[5].re; + t3 = a[1].re - t7; + t7 += a[1].re; + a[1].re = t7; + + t8 = a[7].im; + t6 = a[3].im - t8; + t8 += a[3].im; + a[3].im = t8; + t1 = t3 - t6; + t3 += t6; + + t7 = a[5].im; + t4 = a[1].im - t7; + t7 += a[1].im; + a[1].im = t7; + + t8 = a[7].re; + t5 = a[3].re - t8; + t8 += a[3].re; + a[3].re = t8; + + t2 = t4 - t5; + t4 += t5; + + t6 = t1 - t4; + t8 = sqrthalf; + t6 *= t8; + a[5].re = a[4].re - t6; + t1 += t4; + t1 *= t8; + a[5].im = a[4].im - t1; + t6 += a[4].re; + a[4].re = t6; + t1 += a[4].im; + a[4].im = t1; + + t5 = t2 - t3; + t5 *= t8; + a[7].im = a[6].im - t5; + t2 += t3; + t2 *= t8; + a[7].re = a[6].re - t2; + t2 += a[6].re; + a[6].re = t2; + t5 += a[6].im; + a[6].im = t5; + + c4(a); +} + +static void c16(register WDL_FFT_COMPLEX *a) +{ + register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8; + + TRANSFORMZERO(a[0],a[4],a[8],a[12]); + TRANSFORM(a[1],a[5],a[9],a[13],d16[0].re,d16[0].im); + TRANSFORMHALF(a[2],a[6],a[10],a[14]); + TRANSFORM(a[3],a[7],a[11],a[15],d16[0].im,d16[0].re); + c4(a + 8); + c4(a + 12); + + c8(a); +} + +/* a[0...8n-1], w[0...2n-2]; n >= 2 */ +static void cpass(register WDL_FFT_COMPLEX *a,register const WDL_FFT_COMPLEX *w,register unsigned int n) +{ + register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8; + register WDL_FFT_COMPLEX *a1; + register WDL_FFT_COMPLEX *a2; + register WDL_FFT_COMPLEX *a3; + + a2 = a + 4 * n; + a1 = a + 2 * n; + a3 = a2 + 2 * n; + --n; + + TRANSFORMZERO(a[0],a1[0],a2[0],a3[0]); + TRANSFORM(a[1],a1[1],a2[1],a3[1],w[0].re,w[0].im); + + for (;;) { + TRANSFORM(a[2],a1[2],a2[2],a3[2],w[1].re,w[1].im); + TRANSFORM(a[3],a1[3],a2[3],a3[3],w[2].re,w[2].im); + if (!--n) break; + a += 2; + a1 += 2; + a2 += 2; + a3 += 2; + w += 2; + } +} + +static void c32(register WDL_FFT_COMPLEX *a) +{ + cpass(a,d32,4); + c8(a + 16); + c8(a + 24); + c16(a); +} + +static void c64(register WDL_FFT_COMPLEX *a) +{ + cpass(a,d64,8); + c16(a + 32); + c16(a + 48); + c32(a); +} + +static void c128(register WDL_FFT_COMPLEX *a) +{ + cpass(a,d128,16); + c32(a + 64); + c32(a + 96); + c64(a); +} + +static void c256(register WDL_FFT_COMPLEX *a) +{ + cpass(a,d256,32); + c64(a + 128); + c64(a + 192); + c128(a); +} + +static void c512(register WDL_FFT_COMPLEX *a) +{ + cpass(a,d512,64); + c128(a + 384); + c128(a + 256); + c256(a); +} + +/* a[0...8n-1], w[0...n-2]; n even, n >= 4 */ +static void cpassbig(register WDL_FFT_COMPLEX *a,register const WDL_FFT_COMPLEX *w,register unsigned int n) +{ + register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8; + register WDL_FFT_COMPLEX *a1; + register WDL_FFT_COMPLEX *a2; + register WDL_FFT_COMPLEX *a3; + register unsigned int k; + + a2 = a + 4 * n; + a1 = a + 2 * n; + a3 = a2 + 2 * n; + k = n - 2; + + TRANSFORMZERO(a[0],a1[0],a2[0],a3[0]); + TRANSFORM(a[1],a1[1],a2[1],a3[1],w[0].re,w[0].im); + a += 2; + a1 += 2; + a2 += 2; + a3 += 2; + + do { + TRANSFORM(a[0],a1[0],a2[0],a3[0],w[1].re,w[1].im); + TRANSFORM(a[1],a1[1],a2[1],a3[1],w[2].re,w[2].im); + a += 2; + a1 += 2; + a2 += 2; + a3 += 2; + w += 2; + } while (k -= 2); + + TRANSFORMHALF(a[0],a1[0],a2[0],a3[0]); + TRANSFORM(a[1],a1[1],a2[1],a3[1],w[0].im,w[0].re); + a += 2; + a1 += 2; + a2 += 2; + a3 += 2; + + k = n - 2; + do { + TRANSFORM(a[0],a1[0],a2[0],a3[0],w[-1].im,w[-1].re); + TRANSFORM(a[1],a1[1],a2[1],a3[1],w[-2].im,w[-2].re); + a += 2; + a1 += 2; + a2 += 2; + a3 += 2; + w -= 2; + } while (k -= 2); +} + + +static void c1024(register WDL_FFT_COMPLEX *a) +{ + cpassbig(a,d1024,128); + c256(a + 768); + c256(a + 512); + c512(a); +} + +static void c2048(register WDL_FFT_COMPLEX *a) +{ + cpassbig(a,d2048,256); + c512(a + 1536); + c512(a + 1024); + c1024(a); +} + +static void c4096(register WDL_FFT_COMPLEX *a) +{ + cpassbig(a,d4096,512); + c1024(a + 3072); + c1024(a + 2048); + c2048(a); +} + +static void c8192(register WDL_FFT_COMPLEX *a) +{ + cpassbig(a,d8192,1024); + c2048(a + 6144); + c2048(a + 4096); + c4096(a); +} + +static void c16384(register WDL_FFT_COMPLEX *a) +{ + cpassbig(a,d16384,2048); + c4096(a + 8192 + 4096); + c4096(a + 8192); + c8192(a); +} + +static void c32768(register WDL_FFT_COMPLEX *a) +{ + cpassbig(a,d32768,4096); + c8192(a + 16384 + 8192); + c8192(a + 16384); + c16384(a); +} + + +/* n even, n > 0 */ +void WDL_fft_complexmul(WDL_FFT_COMPLEX *a,WDL_FFT_COMPLEX *b,int n) +{ + register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8; + if (n<2 || (n&1)) return; + + do { + t1 = a[0].re * b[0].re; + t2 = a[0].im * b[0].im; + t3 = a[0].im * b[0].re; + t4 = a[0].re * b[0].im; + t5 = a[1].re * b[1].re; + t6 = a[1].im * b[1].im; + t7 = a[1].im * b[1].re; + t8 = a[1].re * b[1].im; + t1 -= t2; + t3 += t4; + t5 -= t6; + t7 += t8; + a[0].re = t1; + a[1].re = t5; + a[0].im = t3; + a[1].im = t7; + a += 2; + b += 2; + } while (n -= 2); +} + +void WDL_fft_complexmul2(WDL_FFT_COMPLEX *c, WDL_FFT_COMPLEX *a, WDL_FFT_COMPLEX *b, int n) +{ + register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8; + if (n<2 || (n&1)) return; + + do { + t1 = a[0].re * b[0].re; + t2 = a[0].im * b[0].im; + t3 = a[0].im * b[0].re; + t4 = a[0].re * b[0].im; + t5 = a[1].re * b[1].re; + t6 = a[1].im * b[1].im; + t7 = a[1].im * b[1].re; + t8 = a[1].re * b[1].im; + t1 -= t2; + t3 += t4; + t5 -= t6; + t7 += t8; + c[0].re = t1; + c[1].re = t5; + c[0].im = t3; + c[1].im = t7; + a += 2; + b += 2; + c += 2; + } while (n -= 2); +} +void WDL_fft_complexmul3(WDL_FFT_COMPLEX *c, WDL_FFT_COMPLEX *a, WDL_FFT_COMPLEX *b, int n) +{ + register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8; + if (n<2 || (n&1)) return; + + do { + t1 = a[0].re * b[0].re; + t2 = a[0].im * b[0].im; + t3 = a[0].im * b[0].re; + t4 = a[0].re * b[0].im; + t5 = a[1].re * b[1].re; + t6 = a[1].im * b[1].im; + t7 = a[1].im * b[1].re; + t8 = a[1].re * b[1].im; + t1 -= t2; + t3 += t4; + t5 -= t6; + t7 += t8; + c[0].re += t1; + c[1].re += t5; + c[0].im += t3; + c[1].im += t7; + a += 2; + b += 2; + c += 2; + } while (n -= 2); +} + + +static inline void u4(register WDL_FFT_COMPLEX *a) +{ + register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8; + + t1 = VOL a[1].re; + t3 = a[0].re - t1; + t6 = VOL a[2].re; + t1 += a[0].re; + t8 = a[3].re - t6; + t6 += a[3].re; + a[0].re = t1 + t6; + t1 -= t6; + a[2].re = t1; + + t2 = VOL a[1].im; + t4 = a[0].im - t2; + t2 += a[0].im; + t5 = VOL a[3].im; + a[1].im = t4 + t8; + t4 -= t8; + a[3].im = t4; + + t7 = a[2].im - t5; + t5 += a[2].im; + a[1].re = t3 + t7; + t3 -= t7; + a[3].re = t3; + a[0].im = t2 + t5; + t2 -= t5; + a[2].im = t2; +} + +static void u8(register WDL_FFT_COMPLEX *a) +{ + register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8; + + u4(a); + + t1 = a[5].re; + a[5].re = a[4].re - t1; + t1 += a[4].re; + + t3 = a[7].re; + a[7].re = a[6].re - t3; + t3 += a[6].re; + + t8 = t3 - t1; + t1 += t3; + + t6 = a[2].im - t8; + t8 += a[2].im; + a[2].im = t8; + + t5 = a[0].re - t1; + a[4].re = t5; + t1 += a[0].re; + a[0].re = t1; + + t2 = a[5].im; + a[5].im = a[4].im - t2; + t2 += a[4].im; + + t4 = a[7].im; + a[7].im = a[6].im - t4; + t4 += a[6].im; + + a[6].im = t6; + + t7 = t2 - t4; + t2 += t4; + + t3 = a[2].re - t7; + a[6].re = t3; + t7 += a[2].re; + a[2].re = t7; + + t6 = a[0].im - t2; + a[4].im = t6; + t2 += a[0].im; + a[0].im = t2; + + t6 = sqrthalf; + + t1 = a[5].re; + t2 = a[5].im - t1; + t2 *= t6; + t1 += a[5].im; + t1 *= t6; + t4 = a[7].im; + t3 = a[7].re - t4; + t3 *= t6; + t4 += a[7].re; + t4 *= t6; + + t8 = t3 - t1; + t1 += t3; + t7 = t2 - t4; + t2 += t4; + + t4 = a[3].im - t8; + a[7].im = t4; + t5 = a[1].re - t1; + a[5].re = t5; + t3 = a[3].re - t7; + a[7].re = t3; + t6 = a[1].im - t2; + a[5].im = t6; + + t8 += a[3].im; + a[3].im = t8; + t1 += a[1].re; + a[1].re = t1; + t7 += a[3].re; + a[3].re = t7; + t2 += a[1].im; + a[1].im = t2; +} + +static void u16(register WDL_FFT_COMPLEX *a) +{ + register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8; + + u8(a); + u4(a + 8); + u4(a + 12); + + UNTRANSFORMZERO(a[0],a[4],a[8],a[12]); + UNTRANSFORMHALF(a[2],a[6],a[10],a[14]); + UNTRANSFORM(a[1],a[5],a[9],a[13],d16[0].re,d16[0].im); + UNTRANSFORM(a[3],a[7],a[11],a[15],d16[0].im,d16[0].re); +} + +/* a[0...8n-1], w[0...2n-2] */ +static void upass(register WDL_FFT_COMPLEX *a,register const WDL_FFT_COMPLEX *w,register unsigned int n) +{ + register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8; + register WDL_FFT_COMPLEX *a1; + register WDL_FFT_COMPLEX *a2; + register WDL_FFT_COMPLEX *a3; + + a2 = a + 4 * n; + a1 = a + 2 * n; + a3 = a2 + 2 * n; + n -= 1; + + UNTRANSFORMZERO(a[0],a1[0],a2[0],a3[0]); + UNTRANSFORM(a[1],a1[1],a2[1],a3[1],w[0].re,w[0].im); + + for (;;) { + UNTRANSFORM(a[2],a1[2],a2[2],a3[2],w[1].re,w[1].im); + UNTRANSFORM(a[3],a1[3],a2[3],a3[3],w[2].re,w[2].im); + if (!--n) break; + a += 2; + a1 += 2; + a2 += 2; + a3 += 2; + w += 2; + } +} + +static void u32(register WDL_FFT_COMPLEX *a) +{ + u16(a); + u8(a + 16); + u8(a + 24); + upass(a,d32,4); +} + +static void u64(register WDL_FFT_COMPLEX *a) +{ + u32(a); + u16(a + 32); + u16(a + 48); + upass(a,d64,8); +} + +static void u128(register WDL_FFT_COMPLEX *a) +{ + u64(a); + u32(a + 64); + u32(a + 96); + upass(a,d128,16); +} + +static void u256(register WDL_FFT_COMPLEX *a) +{ + u128(a); + u64(a + 128); + u64(a + 192); + upass(a,d256,32); +} + +static void u512(register WDL_FFT_COMPLEX *a) +{ + u256(a); + u128(a + 256); + u128(a + 384); + upass(a,d512,64); +} + + +/* a[0...8n-1], w[0...n-2]; n even, n >= 4 */ +static void upassbig(register WDL_FFT_COMPLEX *a,register const WDL_FFT_COMPLEX *w,register unsigned int n) +{ + register WDL_FFT_REAL t1, t2, t3, t4, t5, t6, t7, t8; + register WDL_FFT_COMPLEX *a1; + register WDL_FFT_COMPLEX *a2; + register WDL_FFT_COMPLEX *a3; + register unsigned int k; + + a2 = a + 4 * n; + a1 = a + 2 * n; + a3 = a2 + 2 * n; + k = n - 2; + + UNTRANSFORMZERO(a[0],a1[0],a2[0],a3[0]); + UNTRANSFORM(a[1],a1[1],a2[1],a3[1],w[0].re,w[0].im); + a += 2; + a1 += 2; + a2 += 2; + a3 += 2; + + do { + UNTRANSFORM(a[0],a1[0],a2[0],a3[0],w[1].re,w[1].im); + UNTRANSFORM(a[1],a1[1],a2[1],a3[1],w[2].re,w[2].im); + a += 2; + a1 += 2; + a2 += 2; + a3 += 2; + w += 2; + } while (k -= 2); + + UNTRANSFORMHALF(a[0],a1[0],a2[0],a3[0]); + UNTRANSFORM(a[1],a1[1],a2[1],a3[1],w[0].im,w[0].re); + a += 2; + a1 += 2; + a2 += 2; + a3 += 2; + + k = n - 2; + do { + UNTRANSFORM(a[0],a1[0],a2[0],a3[0],w[-1].im,w[-1].re); + UNTRANSFORM(a[1],a1[1],a2[1],a3[1],w[-2].im,w[-2].re); + a += 2; + a1 += 2; + a2 += 2; + a3 += 2; + w -= 2; + } while (k -= 2); +} + + + +static void u1024(register WDL_FFT_COMPLEX *a) +{ + u512(a); + u256(a + 512); + u256(a + 768); + upassbig(a,d1024,128); +} + +static void u2048(register WDL_FFT_COMPLEX *a) +{ + u1024(a); + u512(a + 1024); + u512(a + 1536); + upassbig(a,d2048,256); +} + + +static void u4096(register WDL_FFT_COMPLEX *a) +{ + u2048(a); + u1024(a + 2048); + u1024(a + 3072); + upassbig(a,d4096,512); +} + +static void u8192(register WDL_FFT_COMPLEX *a) +{ + u4096(a); + u2048(a + 4096); + u2048(a + 6144); + upassbig(a,d8192,1024); +} + +static void u16384(register WDL_FFT_COMPLEX *a) +{ + u8192(a); + u4096(a + 8192); + u4096(a + 8192 + 4096); + upassbig(a,d16384,2048); +} + +static void u32768(register WDL_FFT_COMPLEX *a) +{ + u16384(a); + u8192(a + 16384); + u8192(a + 16384 + 8192 ); + upassbig(a,d32768,4096); +} + + +static void __fft_gen(WDL_FFT_COMPLEX *buf, const WDL_FFT_COMPLEX *buf2, int sz, int isfull) +{ + int x; + double div=PI*0.25/(sz+1); + + if (isfull) div*=2.0; + + for (x = 0; x < sz; x ++) + { + if (!(x & 1) || !buf2) + { + buf[x].re = (WDL_FFT_REAL) cos((x+1)*div); + buf[x].im = (WDL_FFT_REAL) sin((x+1)*div); + } + else + { + buf[x].re = buf2[x >> 1].re; + buf[x].im = buf2[x >> 1].im; + } + } +} + +#ifndef WDL_FFT_NO_PERMUTE + +static unsigned int fftfreq_c(unsigned int i,unsigned int n) +{ + unsigned int m; + + if (n <= 2) return i; + + m = n >> 1; + if (i < m) return fftfreq_c(i,m) << 1; + + i -= m; + m >>= 1; + if (i < m) return (fftfreq_c(i,m) << 2) + 1; + i -= m; + return ((fftfreq_c(i,m) << 2) - 1) & (n - 1); +} + +static int _idxperm[2<> 1, quart = half >> 1, eighth = quart >> 1; + const int *permute = WDL_fft_permute_tab(half); + unsigned int i, j; + + WDL_FFT_COMPLEX *p, *q, tw, sum, diff; + WDL_FFT_REAL tw1, tw2; + + if (!isInverse) + { + WDL_fft((WDL_FFT_COMPLEX*)buf, half, isInverse); + r2(buf); + } + else + { + v2(buf); + } + + /* Source: http://www.katjaas.nl/realFFT/realFFT2.html */ + + for (i = 1; i < quart; ++i) + { + p = (WDL_FFT_COMPLEX*)buf + permute[i]; + q = (WDL_FFT_COMPLEX*)buf + permute[half - i]; + +/* tw.re = cos(2*PI * i / len); + tw.im = sin(2*PI * i / len); */ + + if (i < eighth) + { + j = i - 1; + tw.re = d[j].re; + tw.im = d[j].im; + } + else if (i > eighth) + { + j = quart - i - 1; + tw.re = d[j].im; + tw.im = d[j].re; + } + else + { + tw.re = tw.im = sqrthalf; + } + + if (!isInverse) tw.re = -tw.re; + + sum.re = p->re + q->re; + sum.im = p->im + q->im; + diff.re = p->re - q->re; + diff.im = p->im - q->im; + + tw1 = tw.re * sum.im + tw.im * diff.re; + tw2 = tw.im * sum.im - tw.re * diff.re; + + p->re = sum.re - tw1; + p->im = diff.im - tw2; + q->re = sum.re + tw1; + q->im = -(diff.im + tw2); + } + + p = (WDL_FFT_COMPLEX*)buf + permute[i]; + p->re *= 2; + p->im *= -2; + + if (isInverse) WDL_fft((WDL_FFT_COMPLEX*)buf, half, isInverse); +} + +void WDL_real_fft(WDL_FFT_REAL* buf, int len, int isInverse) +{ + switch (len) + { + case 2: if (!isInverse) r2(buf); else v2(buf); break; + case 4: case 8: two_for_one(buf, 0, len, isInverse); break; +#define TMP(x) case x: two_for_one(buf, d##x, len, isInverse); break; + TMP(16) + TMP(32) + TMP(64) + TMP(128) + TMP(256) + TMP(512) + TMP(1024) + TMP(2048) + TMP(4096) + TMP(8192) + TMP(16384) + TMP(32768) +#undef TMP + } +} diff --git a/3rdparty/wdlfft/fft.h b/3rdparty/wdlfft/fft.h new file mode 100644 index 00000000000..d2be245fb04 --- /dev/null +++ b/3rdparty/wdlfft/fft.h @@ -0,0 +1,77 @@ +/* + WDL - fft.h + Copyright (C) 2006 and later Cockos Incorporated + + This software is provided 'as-is', without any express or implied + warranty. In no event will the authors be held liable for any damages + arising from the use of this software. + + Permission is granted to anyone to use this software for any purpose, + including commercial applications, and to alter it and redistribute it + freely, subject to the following restrictions: + + 1. The origin of this software must not be misrepresented; you must not + claim that you wrote the original software. If you use this software + in a product, an acknowledgment in the product documentation would be + appreciated but is not required. + 2. Altered source versions must be plainly marked as such, and must not be + misrepresented as being the original software. + 3. This notice may not be removed or altered from any source distribution. + + + + This file defines the interface to the WDL FFT library. These routines are based on the + DJBFFT library, which are Copyright 1999 D. J. Bernstein, djb@pobox.com + + The DJB FFT web page is: http://cr.yp.to/djbfft.html + +*/ + +#ifndef _WDL_FFT_H_ +#define _WDL_FFT_H_ + +#ifdef __cplusplus +extern "C" { +#endif + +#ifndef WDL_FFT_REALSIZE +#define WDL_FFT_REALSIZE 4 +#endif + +#if WDL_FFT_REALSIZE == 4 +typedef float WDL_FFT_REAL; +#elif WDL_FFT_REALSIZE == 8 +typedef double WDL_FFT_REAL; +#else +#error invalid FFT item size +#endif + +typedef struct { + WDL_FFT_REAL re; + WDL_FFT_REAL im; +} WDL_FFT_COMPLEX; + +extern void WDL_fft_init(); + +extern void WDL_fft_complexmul(WDL_FFT_COMPLEX *dest, WDL_FFT_COMPLEX *src, int len); +extern void WDL_fft_complexmul2(WDL_FFT_COMPLEX *dest, WDL_FFT_COMPLEX *src, WDL_FFT_COMPLEX *src2, int len); +extern void WDL_fft_complexmul3(WDL_FFT_COMPLEX *destAdd, WDL_FFT_COMPLEX *src, WDL_FFT_COMPLEX *src2, int len); + +/* Expects WDL_FFT_COMPLEX input[0..len-1] scaled by 1.0/len, returns +WDL_FFT_COMPLEX output[0..len-1] order by WDL_fft_permute(len). */ +extern void WDL_fft(WDL_FFT_COMPLEX *, int len, int isInverse); + +/* Expects WDL_FFT_REAL input[0..len-1] scaled by 0.5/len, returns +WDL_FFT_COMPLEX output[0..len/2-1], for len >= 4 order by +WDL_fft_permute(len/2). Note that output[len/2].re is stored in +output[0].im. */ +extern void WDL_real_fft(WDL_FFT_REAL *, int len, int isInverse); + +extern int WDL_fft_permute(int fftsize, int idx); +extern int *WDL_fft_permute_tab(int fftsize); + +#ifdef __cplusplus +}; +#endif + +#endif \ No newline at end of file -- cgit v1.2.3-70-g09d2