FFT processor

LibSource/FFTProcessor.h

@@ -0,0 +1,124 @@
+#ifndef __FFT_PROCESSOR_H__
+#define __FFT_PROCESSOR_H__
+
+#include "SignalProcessor.h"
+#include "Window.h"
+
+template <class FrequencyDomainProcessor, size_t fft_size, size_t window_size>
+class FFTProcessor : public SignalProcessor {
+public:
+    FFTProcessor(FastFourierTransform* fft, FrequencyDomainProcessor* processor,
+        FloatArray timedomain, ComplexFloatArray freqdomain, FloatArray in_overlap,
+        FloatArray out_overlap, Window in_win, Window out_win)
+        : fft(fft)
+        , index(0)
+        , timedomain(timedomain)
+        , freqdomain(freqdomain)
+        , in_overlap(in_overlap)
+        , out_overlap(out_overlap)
+        , in_win(in_win)
+        , out_win(out_win)
+        , processor(processor) {
+    }
+    FrequencyDomainProcessor* getFrequencyDomainProcessor() {
+        return processor;
+    }
+    /**
+     * We don't support sample based processing in this class
+     **/
+    float process(float input) override {
+        return 0.0;
+    }
+    void process(FloatArray input, FloatArray output) override{
+        const size_t block_size = input.getSize();
+
+        // Insert twice to prepare contiguous copy of data when we reach buffer end
+        in_overlap.insert(input, index, block_size);
+        in_overlap.insert(input, (index + window_size) % (window_size * 2), block_size);
+
+        // When end is reach, we rewind to middle of the buffer
+        if (index >= window_size * 2 - block_size) {
+            index = window_size;
+        }
+        else {
+            index += block_size;
+        }
+
+        // Process data only when there's enough input to fill first window
+        if (index >= window_size) {
+            // FFT
+            in_win.apply(in_overlap.getData() + index - window_size, timedomain);
+            fft->fft(timedomain, freqdomain);
+            ComplexFloatArray bins = freqdomain.subArray(0, fft_size / 2);
+            // Time domain processing
+            processor->process(bins, bins);
+            // IFFT
+            fft->ifft(freqdomain, timedomain);
+            out_win.apply(timedomain);
+            // Overlap accumulation
+            const size_t out_index = index - window_size;
+            output.copyFrom(timedomain.subArray(0, block_size));
+            FloatArray t = out_overlap.subArray(out_index, block_size);
+            output.add(t);
+            t.setAll(0);
+            out_overlap
+                .subArray(out_index + block_size, window_size - out_index - block_size)
+                .add(timedomain.subArray(
+                    block_size, window_size - out_index - block_size));
+            out_overlap.subArray(0, out_index)
+                .add(timedomain.subArray(window_size - out_index, out_index));
+        }
+    }
+    template <typename... Args>
+    static FFTProcessor* create(size_t block_size,
+        Window::WindowType in_win_type = Window::HannWindow,
+        Window::WindowType out_win_type = Window::HannWindow, Args&&... args) {
+        FrequencyDomainProcessor* processor =
+            FrequencyDomainProcessor::create(std::forward<Args>(args)...);
+        Window in_win = Window::create(in_win_type, window_size);
+        Window out_win = Window::create(out_win_type, window_size);
+
+        float norm[block_size];
+        memset(norm, 0, block_size * sizeof(float));
+
+        // Accumulate shifted multiplications of both windows as norm
+        for (size_t i = 0; i < window_size; i++) {
+            norm[i % block_size] += in_win[i] * out_win[i];
+        }
+
+        // Divide output window by the norm for every block size, compensating amplitude changes
+        float* t = out_win.getData();
+        for (size_t i = 0; i < window_size / block_size; i++) {
+            for (size_t j = 0; j < block_size; j++) {
+                *t++ /= norm[j];
+            }
+        }
+
+        return new FFTProcessor(FastFourierTransform::create(fft_size), processor,
+            FloatArray::create(fft_size), ComplexFloatArray::create(fft_size),
+            FloatArray::create(window_size * 2),
+            FloatArray::create(window_size), in_win, out_win);
+    }
+    static void destroy(FFTProcessor* fft) {
+        FloatArray::destroy(fft->timedomain);
+        ComplexFloatArray::destroy(fft->freqdomain);
+        FloatArray::destroy(fft->in_overlap);
+        FloatArray::destroy(fft->out_overlap);
+        Window::destroy(fft->in_win);
+        Window::destroy(fft->out_win);
+        FrequencyDomainProcessor::destroy(fft->processor);
+        FastFourierTransform::destroy(fft->fft);
+        delete fft;
+    }
+
+protected:
+    size_t index;
+    FastFourierTransform* fft;
+    FloatArray timedomain;
+    ComplexFloatArray freqdomain;
+    Window in_win, out_win;
+    FloatArray in_overlap, out_overlap;
+    FrequencyDomainProcessor* processor;
+};
+
+#endif

LibSource/OpenWareLibrary.h

@@ -19,6 +19,7 @@
 #include "ExponentialDecayEnvelope.h"
 #include "FastFourierTransform.h"
 #include "FeedbackProcessor.h"
+#include "FFTProcessor.h"
 #include "FractionalCircularBuffer.h"
 #include "FirFilter.h"
 #include "FloatArray.h"

LibSource/Window.h

@@ -7,7 +7,8 @@
 
 /*
  * Window provides static methods to generate and apply window functions:
- * rectangular, hann, hanning, hamming, triangular
+ * rectangular, hann, hamming, triangular, sine, gauss, blackman, blackman_harris,
+ * nuttall, blackman_nuttall, bohman, flattop, lanczos, parzen, welch
  *   or
  * window(WindowType type, float* window, int size);
  * a static method to apply a window to a signal (nothing but pairwise multiplication)
@@ -16,6 +17,8 @@
  * (as computing a triangular window is very cheap, 
  * this solution can be handy as it requires less memory (no window array required))
  * applyTriangularWindow(float *signal, int size)
+ * 
+ * Gauss and Tukey windows accept an extra parameter that modifies their shape
  */
 class Window : public FloatArray, SignalProcessor {
 private:
@@ -24,7 +27,20 @@ class Window : public FloatArray, SignalProcessor {
   typedef enum WindowType {
     HammingWindow,
     HannWindow,
-    HanningWindow,
+    SineWindow,
+    BartlettHannWindow,
+    GaussWindow,
+    TukeyWindow,
+    BlackmanWindow,
+    BlackmanHarrisWindow,
+    NuttallWindow,
+    BlackmanNuttallWindow,
+    BohmanWindow,
+    FlatTopWindow,
+    LanczosWindow,
+    ParzenWindow,
+    WelchWindow,
+    BartlettWindow,
     TriangularWindow,
     RectangularWindow // no window
   } WindowType;
@@ -43,27 +59,68 @@ class Window : public FloatArray, SignalProcessor {
     return value;
   }
   void process(FloatArray input, FloatArray output){
-    Window::applyWindow(input, getData(), output, getSize());    
+    Window::applyWindow(input, getData(), output, getSize());
   }
   static Window create(WindowType type, int size){
     Window win(new float[size], size);
     win.window(type, win, size);
     return win;
   }
   static Window create(int size){
-    Window win(new float[size], size);
-    return win;
+    return Window(new float[size], size);
+  }
+  static void destroy(Window win) {
+    delete[] win.getData();
   }
   static void window(WindowType type, float *window, int size){
     switch(type){
     case HannWindow:
-    case HanningWindow:
       hann(window, size);
       break;
     case HammingWindow:
       hamming(window, size);
       break;
+    case SineWindow:
+      sine(window, size);
+      break;
+    case BartlettHannWindow:
+      bartlett_hann(window, size);
+      break;
+    case GaussWindow:
+      gauss(window, size);
+      break;
+    case TukeyWindow:
+      tukey(window, size);
+      break;
+    case BlackmanWindow:
+      blackman(window, size);
+      break;
+    case BlackmanHarrisWindow:
+      blackman_harris(window, size);
+      break;
+    case NuttallWindow:
+      nuttall(window, size);
+      break;
+    case BlackmanNuttallWindow:
+      blackman_nuttall(window, size);
+      break;
+    case BohmanWindow:
+      bohman(window, size);
+      break;
+    case FlatTopWindow:
+      flat_top(window, size);
+      break;
+    case LanczosWindow:
+      lanczos(window, size);
+      break;
+    case ParzenWindow:
+      parzen(window, size);
+      break;
+    case WelchWindow:
+      welch(window, size);
+      break;
     case TriangularWindow:
+    case BartlettWindow:
       triangular(window, size);
       break;
     case RectangularWindow:
@@ -84,6 +141,97 @@ class Window : public FloatArray, SignalProcessor {
     for(int n=0; n<size; n++)
       window[n] = 0.54-0.46*cosf((float)n/(size-1)*2*M_PI);
   }
+  static void sine(float* window, int size) {
+    for(int n=0; n<size; n++)
+      window[n] = sinf(M_PI * n / (size - 1));
+  }
+  static void bartlett_hann(float* window, int size) {
+    for(int n=0; n<size; n++)
+      window[n] = 0.62 - 0.48 * std::abs((float)n / (size - 1) - 0.5)
+        + 0.38 * cosf(M_PI * 2 * (float(n) / (size - 1) - 0.5));
+  }
+  static void gauss(float *window, int size, float q = 0.5){
+    float a = (size - 1) / 2;
+    for(int n=0; n<size; n++) {
+      float t = (n - a) / (q * a);
+      window[n] = expf(-t * t / 2);
+    }
+  }
+  static void tukey(float *window, int size, float q = 0.5){
+    float a = q * 0.5 * (size - 1);
+    for(int n=0; n<size; n++) {
+      if (n <= a) {
+        window[n] =  0.5 * (1.0 + cosf(M_PI * (float(n) /a - 1.0)));
+      }
+      else if (n >= 1.0 - a) {
+        window[n] = 0.5 * (1.0 + cosf(M_PI * (float(n) / a - 1.0)));
+      }
+      else {
+        window[n] = 1;
+      }
+    }
+  }
+  static void blackman(float *window, int size){
+    for(int n=0; n<size; n++) {
+      float a = M_PI * 2 * n / (size - 1);
+      window[n] = 0.42 - 0.5 * cosf(a) + 0.08 * cosf(a * 2);
+    }
+  }
+  static void blackman_harris(float *window, int size){
+    for(int n=0; n<size; n++)
+      window[n] = 0.35875 - (0.48829 * cosf((2 * M_PI * n)/(size - 1))) +
+        (0.14128 * cosf((4 * M_PI * n) / (size - 1))) - (0.01168 * cosf((4 * M_PI * n) / (size - 1)));
+  }
+  static void nuttall(float *window, int size){
+    for(int n=0; n<size; n++) {
+      float x = M_PI * 2 * n / (size - 1);
+      window[n] = 0.355768 - 0.487396 * cosf(x * 2) + 0.144232 * cosf(x * 4) - 0.012604 * cosf(x * 6);
+    }
+  }
+  static void blackman_nuttall(float *window, int size){
+    for(int n=0; n<size; n++){
+      float x = float(n) / (size - 1) * M_PI;
+      window[n] = 0.3635819 - 0.4891775 * cosf(x * 2) + 0.1365995 * cosf(x * 4) - 0.0106411 * cosf(x * 6);
+    }
+  }
+  static void bohman(float* window, int size) {
+    for(int n=0; n<size; n++) {
+      float x = float(n) / (size - 1) * 2 - 1;
+      window[n] = (1.0 - std::abs(x)) * cosf(M_PI * std::abs(x)) + M_1_PI * sinf(M_PI * std::abs(x));
+    }
+  }
+  static void flat_top(float* window, int size) {
+    for(int n=0; n<size; n++) {
+      float x = float(n) / (size - 1) * M_PI;
+      window[n] = 0.21557895 - 0.41663158 * cosf(x * 2)
+        + 0.277263158 * cosf(x * 4) - 0.083578947 * cosf(x * 6)
+        + 0.006947368 * cosf(x * 8);
+    }
+  }
+  static void parzen(float* window, int size) {
+    for(int n=0; n<size; n++) {
+      int n1 = n - size / 2;
+      if (std::abs(n1) < size / 4) {
+        float a = n1 * 2.0 / size;
+        window[n] = 1.0 - 7.0 * a * a * (1.0 - std::abs(n1) * 2.0 / size);
+      }
+      else {
+        float a = (1.0 - std::abs(n1) * 2.0 / size);
+        window[n] = 2.0 * a * a * a;
+      }
+    }
+  }
+  static void welch(float* window, int size) {
+    for(int n=0; n<size; n++) {
+      float a = float(size - 1) / 2;
+      float b = (float(n) - a) / a;
+      window[n] = 1.0 - a * a;
+    }
+  }
+  static void lanczos(float* window, int size) {
+    for(int n=0; n<size; n++)
+      window[n] = sinf(M_PI * 2 * n / (size - 1) - M_PI) / (M_PI * 2 * n / (size - 1) - M_PI);
+  }
   static void triangular(float *window, int size){
     float rampStep = 1/(size/2.0f);
     float ramp = 0;