diff --git a/LibSource/FFTProcessor.h b/LibSource/FFTProcessor.h new file mode 100644 index 00000000..72308fe2 --- /dev/null +++ b/LibSource/FFTProcessor.h @@ -0,0 +1,124 @@ +#ifndef __FFT_PROCESSOR_H__ +#define __FFT_PROCESSOR_H__ + +#include "SignalProcessor.h" +#include "Window.h" + +template +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 + 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)...); + 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 diff --git a/LibSource/OpenWareLibrary.h b/LibSource/OpenWareLibrary.h index a94871c1..cd3339b5 100644 --- a/LibSource/OpenWareLibrary.h +++ b/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" diff --git a/LibSource/Window.h b/LibSource/Window.h index f0f6d7de..62a04d85 100644 --- a/LibSource/Window.h +++ b/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,7 +59,7 @@ 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); @@ -51,19 +67,60 @@ class Window : public FloatArray, SignalProcessor { 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= 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