Audio plugin host https://kx.studio/carla
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

juce_FFT.cpp 9.8KB

9 years ago
9 years ago
9 years ago
9 years ago
9 years ago
9 years ago
9 years ago
9 years ago
9 years ago
9 years ago
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305
  1. /*
  2. ==============================================================================
  3. This file is part of the JUCE library.
  4. Copyright (c) 2015 - ROLI Ltd.
  5. Permission is granted to use this software under the terms of either:
  6. a) the GPL v2 (or any later version)
  7. b) the Affero GPL v3
  8. Details of these licenses can be found at: www.gnu.org/licenses
  9. JUCE is distributed in the hope that it will be useful, but WITHOUT ANY
  10. WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  11. A PARTICULAR PURPOSE. See the GNU General Public License for more details.
  12. ------------------------------------------------------------------------------
  13. To release a closed-source product which uses JUCE, commercial licenses are
  14. available: visit www.juce.com for more information.
  15. ==============================================================================
  16. */
  17. // (For the moment, we'll implement a few local operators for this complex class - one
  18. // day we'll probably either have a juce complex class, or use the C++11 one)
  19. static FFT::Complex operator+ (FFT::Complex a, FFT::Complex b) noexcept { FFT::Complex c = { a.r + b.r, a.i + b.i }; return c; }
  20. static FFT::Complex operator- (FFT::Complex a, FFT::Complex b) noexcept { FFT::Complex c = { a.r - b.r, a.i - b.i }; return c; }
  21. static FFT::Complex operator* (FFT::Complex a, FFT::Complex b) noexcept { FFT::Complex c = { a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r }; return c; }
  22. static FFT::Complex& operator+= (FFT::Complex& a, FFT::Complex b) noexcept { a.r += b.r; a.i += b.i; return a; }
  23. //==============================================================================
  24. struct FFT::FFTConfig
  25. {
  26. FFTConfig (int sizeOfFFT, bool isInverse)
  27. : fftSize (sizeOfFFT), inverse (isInverse), twiddleTable ((size_t) sizeOfFFT)
  28. {
  29. for (int i = 0; i < fftSize; ++i)
  30. {
  31. const double phase = (isInverse ? 2.0 : -2.0) * double_Pi * i / fftSize;
  32. twiddleTable[i].r = (float) cos (phase);
  33. twiddleTable[i].i = (float) sin (phase);
  34. }
  35. const int root = (int) std::sqrt ((double) fftSize);
  36. int divisor = 4, n = fftSize;
  37. for (int i = 0; i < numElementsInArray (factors); ++i)
  38. {
  39. while ((n % divisor) != 0)
  40. {
  41. if (divisor == 2) divisor = 3;
  42. else if (divisor == 4) divisor = 2;
  43. else divisor += 2;
  44. if (divisor > root)
  45. divisor = n;
  46. }
  47. n /= divisor;
  48. jassert (divisor == 1 || divisor == 2 || divisor == 4);
  49. factors[i].radix = divisor;
  50. factors[i].length = n;
  51. }
  52. }
  53. void perform (const Complex* input, Complex* output) const noexcept
  54. {
  55. perform (input, output, 1, 1, factors);
  56. }
  57. const int fftSize;
  58. const bool inverse;
  59. struct Factor { int radix, length; };
  60. Factor factors[32];
  61. HeapBlock<Complex> twiddleTable;
  62. void perform (const Complex* input, Complex* output, const int stride, const int strideIn, const Factor* facs) const noexcept
  63. {
  64. const Factor factor (*facs++);
  65. Complex* const originalOutput = output;
  66. const Complex* const outputEnd = output + factor.radix * factor.length;
  67. if (stride == 1 && factor.radix <= 5)
  68. {
  69. for (int i = 0; i < factor.radix; ++i)
  70. perform (input + stride * strideIn * i, output + i * factor.length, stride * factor.radix, strideIn, facs);
  71. butterfly (factor, output, stride);
  72. return;
  73. }
  74. if (factor.length == 1)
  75. {
  76. do
  77. {
  78. *output++ = *input;
  79. input += stride * strideIn;
  80. }
  81. while (output < outputEnd);
  82. }
  83. else
  84. {
  85. do
  86. {
  87. perform (input, output, stride * factor.radix, strideIn, facs);
  88. input += stride * strideIn;
  89. output += factor.length;
  90. }
  91. while (output < outputEnd);
  92. }
  93. butterfly (factor, originalOutput, stride);
  94. }
  95. void butterfly (const Factor factor, Complex* data, const int stride) const noexcept
  96. {
  97. switch (factor.radix)
  98. {
  99. case 1: break;
  100. case 2: butterfly2 (data, stride, factor.length); return;
  101. case 4: butterfly4 (data, stride, factor.length); return;
  102. default: jassertfalse; break;
  103. }
  104. Complex* scratch = static_cast<Complex*> (alloca (sizeof (Complex) * (size_t) factor.radix));
  105. for (int i = 0; i < factor.length; ++i)
  106. {
  107. for (int k = i, q1 = 0; q1 < factor.radix; ++q1)
  108. {
  109. scratch[q1] = data[k];
  110. k += factor.length;
  111. }
  112. for (int k = i, q1 = 0; q1 < factor.radix; ++q1)
  113. {
  114. int twiddleIndex = 0;
  115. data[k] = scratch[0];
  116. for (int q = 1; q < factor.radix; ++q)
  117. {
  118. twiddleIndex += stride * k;
  119. if (twiddleIndex >= fftSize)
  120. twiddleIndex -= fftSize;
  121. data[k] += scratch[q] * twiddleTable[twiddleIndex];
  122. }
  123. k += factor.length;
  124. }
  125. }
  126. }
  127. void butterfly2 (Complex* data, const int stride, const int length) const noexcept
  128. {
  129. Complex* dataEnd = data + length;
  130. const Complex* tw = twiddleTable;
  131. for (int i = length; --i >= 0;)
  132. {
  133. const Complex s (*dataEnd * *tw);
  134. tw += stride;
  135. *dataEnd++ = *data - s;
  136. *data++ += s;
  137. }
  138. }
  139. void butterfly4 (Complex* data, const int stride, const int length) const noexcept
  140. {
  141. const int lengthX2 = length * 2;
  142. const int lengthX3 = length * 3;
  143. const Complex* twiddle1 = twiddleTable;
  144. const Complex* twiddle2 = twiddle1;
  145. const Complex* twiddle3 = twiddle1;
  146. for (int i = length; --i >= 0;)
  147. {
  148. const Complex s0 = data[length] * *twiddle1;
  149. const Complex s1 = data[lengthX2] * *twiddle2;
  150. const Complex s2 = data[lengthX3] * *twiddle3;
  151. const Complex s3 = s0 + s2;
  152. const Complex s4 = s0 - s2;
  153. const Complex s5 = *data - s1;
  154. *data += s1;
  155. data[lengthX2] = *data - s3;
  156. twiddle1 += stride;
  157. twiddle2 += stride * 2;
  158. twiddle3 += stride * 3;
  159. *data += s3;
  160. if (inverse)
  161. {
  162. data[length].r = s5.r - s4.i;
  163. data[length].i = s5.i + s4.r;
  164. data[lengthX3].r = s5.r + s4.i;
  165. data[lengthX3].i = s5.i - s4.r;
  166. }
  167. else
  168. {
  169. data[length].r = s5.r + s4.i;
  170. data[length].i = s5.i - s4.r;
  171. data[lengthX3].r = s5.r - s4.i;
  172. data[lengthX3].i = s5.i + s4.r;
  173. }
  174. ++data;
  175. }
  176. }
  177. JUCE_DECLARE_NON_COPYABLE_WITH_LEAK_DETECTOR (FFTConfig)
  178. };
  179. //==============================================================================
  180. FFT::FFT (int order, bool inverse) : config (new FFTConfig (1 << order, inverse)), size (1 << order) {}
  181. FFT::~FFT() {}
  182. void FFT::perform (const Complex* const input, Complex* const output) const noexcept
  183. {
  184. config->perform (input, output);
  185. }
  186. const size_t maxFFTScratchSpaceToAlloca = 256 * 1024;
  187. void FFT::performRealOnlyForwardTransform (float* d) const noexcept
  188. {
  189. const size_t scratchSize = 16 + sizeof (FFT::Complex) * (size_t) size;
  190. if (scratchSize < maxFFTScratchSpaceToAlloca)
  191. {
  192. performRealOnlyForwardTransform (static_cast<Complex*> (alloca (scratchSize)), d);
  193. }
  194. else
  195. {
  196. HeapBlock<char> heapSpace (scratchSize);
  197. performRealOnlyForwardTransform (reinterpret_cast<Complex*> (heapSpace.getData()), d);
  198. }
  199. }
  200. void FFT::performRealOnlyInverseTransform (float* d) const noexcept
  201. {
  202. const size_t scratchSize = 16 + sizeof (FFT::Complex) * (size_t) size;
  203. if (scratchSize < maxFFTScratchSpaceToAlloca)
  204. {
  205. performRealOnlyInverseTransform (static_cast<Complex*> (alloca (scratchSize)), d);
  206. }
  207. else
  208. {
  209. HeapBlock<char> heapSpace (scratchSize);
  210. performRealOnlyInverseTransform (reinterpret_cast<Complex*> (heapSpace.getData()), d);
  211. }
  212. }
  213. void FFT::performRealOnlyForwardTransform (Complex* scratch, float* d) const noexcept
  214. {
  215. // This can only be called on an FFT object that was created to do forward transforms.
  216. jassert (! config->inverse);
  217. for (int i = 0; i < size; ++i)
  218. {
  219. scratch[i].r = d[i];
  220. scratch[i].i = 0;
  221. }
  222. perform (scratch, reinterpret_cast<Complex*> (d));
  223. }
  224. void FFT::performRealOnlyInverseTransform (Complex* scratch, float* d) const noexcept
  225. {
  226. // This can only be called on an FFT object that was created to do inverse transforms.
  227. jassert (config->inverse);
  228. perform (reinterpret_cast<const Complex*> (d), scratch);
  229. const float scaleFactor = 1.0f / size;
  230. for (int i = 0; i < size; ++i)
  231. {
  232. d[i] = scratch[i].r * scaleFactor;
  233. d[i + size] = scratch[i].i * scaleFactor;
  234. }
  235. }
  236. void FFT::performFrequencyOnlyForwardTransform (float* d) const noexcept
  237. {
  238. performRealOnlyForwardTransform (d);
  239. const int twiceSize = size * 2;
  240. for (int i = 0; i < twiceSize; i += 2)
  241. {
  242. d[i / 2] = juce_hypot (d[i], d[i + 1]);
  243. if (i >= size)
  244. {
  245. d[i] = 0;
  246. d[i + 1] = 0;
  247. }
  248. }
  249. }