Add RealTimeConvolver implementation (moved from Befaco), added

attribution for pffft

LICENSE-dist.txt

@@ -416,3 +416,51 @@ THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR IMPLIED WA
  * [including the GNU Public Licence.]
  */
 
+
+# pffft
+
+Copyright (c) 2013  Julien Pommier ( pommier@modartt.com ) 
+
+Based on original fortran 77 code from FFTPACKv4 from NETLIB,
+authored by Dr Paul Swarztrauber of NCAR, in 1985.
+
+As confirmed by the NCAR fftpack software curators, the following
+FFTPACKv5 license applies to FFTPACKv4 sources. My changes are
+released under the same terms.
+
+FFTPACK license:
+
+http://www.cisl.ucar.edu/css/software/fftpack5/ftpk.html
+
+Copyright (c) 2004 the University Corporation for Atmospheric
+Research ("UCAR"). All rights reserved. Developed by NCAR's
+Computational and Information Systems Laboratory, UCAR,
+www.cisl.ucar.edu.
+
+Redistribution and use of the Software in source and binary forms,
+with or without modification, is permitted provided that the
+following conditions are met:
+
+- Neither the names of NCAR's Computational and Information Systems
+Laboratory, the University Corporation for Atmospheric Research,
+nor the names of its sponsors or contributors may be used to
+endorse or promote products derived from this Software without
+specific prior written permission.  
+
+- Redistributions of source code must retain the above copyright
+notices, this list of conditions, and the disclaimer below.
+
+- Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions, and the disclaimer below in the
+documentation and/or other materials provided with the
+distribution.
+
+THIS SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR COPYRIGHT
+HOLDERS BE LIABLE FOR ANY CLAIM, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES OR OTHER LIABILITY, WHETHER IN AN
+ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS WITH THE
+SOFTWARE.

include/dsp/fir.hpp

@@ -1,5 +1,6 @@
 #pragma once
 #include "dsp/functions.hpp"
+#include "pffft.h"
 
 
 namespace rack {
@@ -35,4 +36,112 @@ inline void boxcarFIR(float *x, int n, float cutoff) {
 	}
 }
 
+
+struct RealTimeConvolver {
+	// `kernelBlocks` number of contiguous FFT blocks of size `blockSize`
+	// indexed by [i * blockSize*2 + j]
+	float *kernelFfts = NULL;
+	float *inputFfts = NULL;
+	float *outputTail = NULL;
+	float *tmpBlock = NULL;
+	size_t blockSize;
+	size_t kernelBlocks = 0;
+	size_t inputPos = 0;
+	PFFFT_Setup *pffft;
+
+	/** `blockSize` is the size of each FFT block. It should be >=32 and a power of 2. */
+	RealTimeConvolver(size_t blockSize) {
+		this->blockSize = blockSize;
+		pffft = pffft_new_setup(blockSize*2, PFFFT_REAL);
+		outputTail = new float[blockSize]();
+		tmpBlock = new float[blockSize*2]();
+	}
+
+	~RealTimeConvolver() {
+		clear();
+		delete[] outputTail;
+		delete[] tmpBlock;
+		pffft_destroy_setup(pffft);
+	}
+
+	void clear() {
+		if (kernelFfts) {
+			pffft_aligned_free(kernelFfts);
+			kernelFfts = NULL;
+		}
+		if (inputFfts) {
+			pffft_aligned_free(inputFfts);
+			inputFfts = NULL;
+		}
+		kernelBlocks = 0;
+		inputPos = 0;
+	}
+
+	void setKernel(const float *kernel, size_t length) {
+		clear();
+
+		assert(kernel);
+		assert(length > 0);
+
+		// Round up to the nearest factor of `blockSize`
+		kernelBlocks = (length - 1) / blockSize + 1;
+
+		// Allocate blocks
+		kernelFfts = (float*) pffft_aligned_malloc(sizeof(float) * blockSize*2 * kernelBlocks);
+		inputFfts = (float*) pffft_aligned_malloc(sizeof(float) * blockSize*2 * kernelBlocks);
+		memset(inputFfts, 0, sizeof(float) * blockSize*2 * kernelBlocks);
+
+		for (size_t i = 0; i < kernelBlocks; i++) {
+			// Pad each block with zeros
+			memset(tmpBlock, 0, sizeof(float) * blockSize*2);
+			size_t len = min((int) blockSize, (int) (length - i*blockSize));
+			memcpy(tmpBlock, &kernel[i*blockSize], sizeof(float)*len);
+			// Compute fft
+			pffft_transform(pffft, tmpBlock, &kernelFfts[blockSize*2 * i], NULL, PFFFT_FORWARD);
+		}
+	}
+
+	/** Applies reverb to input
+	input and output must be size blockSize
+	*/
+	void processBlock(const float *input, float *output) {
+		if (kernelBlocks == 0) {
+			memset(output, 0, sizeof(float) * blockSize);
+			return;
+		}
+
+		// Step input position
+		inputPos = (inputPos + 1) % kernelBlocks;
+		// Pad block with zeros
+		memset(tmpBlock, 0, sizeof(float) * blockSize*2);
+		memcpy(tmpBlock, input, sizeof(float) * blockSize);
+		// Compute input fft
+		pffft_transform(pffft, tmpBlock, &inputFfts[blockSize*2 * inputPos], NULL, PFFFT_FORWARD);
+		// Create output fft
+		memset(tmpBlock, 0, sizeof(float) * blockSize*2);
+		// convolve input fft by kernel fft
+		// Note: This is the CPU bottleneck loop
+		for (size_t i = 0; i < kernelBlocks; i++) {
+			size_t pos = (inputPos - i + kernelBlocks) % kernelBlocks;
+			pffft_zconvolve_accumulate(pffft, &kernelFfts[blockSize*2 * i], &inputFfts[blockSize*2 * pos], tmpBlock, 1.0);
+		}
+		// Compute output
+		pffft_transform(pffft, tmpBlock, tmpBlock, NULL, PFFFT_BACKWARD);
+		// Add block tail from last output block
+		for (size_t i = 0; i < blockSize; i++) {
+			tmpBlock[i] += outputTail[i];
+		}
+		// Copy output block to output
+		for (size_t i = 0; i < blockSize; i++) {
+			// Scale based on FFT
+			output[i] = tmpBlock[i] / blockSize;
+		}
+		// Set tail
+		for (size_t i = 0; i < blockSize; i++) {
+			outputTail[i] = tmpBlock[i + blockSize];
+		}
+	}
+};
+
+
 } // namespace rack