From db6d6d68b46b6707f48fd8810bbbd6122c11ae2e Mon Sep 17 00:00:00 2001
From: Chris Robinson <chris.kcat@gmail.com>
Date: Sat, 7 Nov 2020 17:59:03 -0800
Subject: Avoid bitshifts for loop counting

---
 common/alcomplex.cpp | 32 +++++++++++++++++++++-----------
 1 file changed, 21 insertions(+), 11 deletions(-)

(limited to 'common/alcomplex.cpp')

diff --git a/common/alcomplex.cpp b/common/alcomplex.cpp
index 458ac99e..8a823b01 100644
--- a/common/alcomplex.cpp
+++ b/common/alcomplex.cpp
@@ -8,33 +8,41 @@
 #include <cstddef>
 #include <utility>
 
+#include "alnumeric.h"
 #include "math_defs.h"
 
 
 void complex_fft(const al::span<std::complex<double>> buffer, const double sign)
 {
     const size_t fftsize{buffer.size()};
-    /* Bit-reversal permutation applied to a sequence of FFTSize items */
-    for(size_t i{1u};i < fftsize-1;i++)
+    /* Get the number of bits used for indexing. Simplifies bit-reversal and
+     * the main loop count.
+     */
+    const size_t log2_size{static_cast<size_t>(CountTrailingZeros(fftsize))};
+
+    /* Bit-reversal permutation applied to a sequence of fftsize items. */
+    for(size_t idx{1u};idx < fftsize-1;++idx)
     {
-        size_t j{0u};
-        for(size_t imask{i + fftsize};imask;imask >>= 1)
-            j = (j<<1) + (imask&1);
-        j >>= 1;
+        size_t revidx{0u}, imask{idx};
+        for(size_t i{0};i < log2_size;++i)
+        {
+            revidx = (revidx<<1) | (imask&1);
+            imask >>= 1;
+        }
 
-        if(i < j)
-            std::swap(buffer[i], buffer[j]);
+        if(idx < revidx)
+            std::swap(buffer[idx], buffer[revidx]);
     }
 
     /* Iterative form of Danielson-Lanczos lemma */
-    size_t step{2u};
-    for(size_t i{1u};i < fftsize;i<<=1, step<<=1)
+    size_t step2{1u};
+    for(size_t i{0};i < log2_size;++i)
     {
-        const size_t step2{step >> 1};
         const double arg{al::MathDefs<double>::Pi() / static_cast<double>(step2)};
 
         const std::complex<double> w{std::cos(arg), std::sin(arg)*sign};
         std::complex<double> u{1.0, 0.0};
+        const size_t step{step2 << 1};
         for(size_t j{0};j < step2;j++)
         {
             for(size_t k{j};k < fftsize;k+=step)
@@ -46,6 +54,8 @@ void complex_fft(const al::span<std::complex<double>> buffer, const double sign)
 
             u *= w;
         }
+
+        step2 <<= 1;
     }
 }
 
-- 
cgit v1.2.3