Initial commit

.gitignore

@@ -0,0 +1 @@
+build/

CMakeLists.txt

@@ -0,0 +1,16 @@
+cmake_minimum_required(VERSION 2.8)
+project(SlidingDFT_Test CXX)
+
+include_directories(.)
+add_executable(SlidingDFT_Test sliding_dft.hpp testbench/main.cpp testbench/test1_data.h)
+
+# ensure C++11
+include(CheckCXXCompilerFlag)
+CHECK_CXX_COMPILER_FLAG("-std=c++11" COMPILER_SUPPORTS_CXX11)
+if(COMPILER_SUPPORTS_CXX11)
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+endif()
+
+# on modern CMake, do this to ensure C++11
+set_property(TARGET SlidingDFT_Test PROPERTY CXX_STANDARD 11)
+set_property(TARGET SlidingDFT_Test PROPERTY CXX_STANDARD_REQUIRED ON)

README.md

@@ -0,0 +1,74 @@
+Sliding discrete Fourier transform (C++)
+====
+
+This code efficiently computes discrete Fourier transforms (DFTs) from a
+continuous sequence of input values. It is a recursive algorithm that updates
+the DFT when each new time-domain measurement arrives, effectively applying a
+sliding window over the last *N* samples. This implementation applies the
+Hanning window in order to minimise spectral leakage.
+
+The update step has computational complexity *O(N)*. If a new DFT is required
+every *M* samples, and *M* < log2(*N*), then this approach is more efficient
+that recalculating the DFT from scratch each time.
+
+This is a header-only C++ library. Simply copy sliding_dft.hpp into your
+project, and use it as follows:
+
+    // Use double precision arithmetic and a 512-length DFT
+    static SlidingDFT<double, 512> dft;
+    // avoid allocating on the stack because the object is large
+
+    // When a new time sample arrives, update the DFT with:
+    dft.update(x);
+
+    // After at least 512 samples have been processed:
+    std::complex<double> DC_bin = dft.dft[0];
+
+Your application should call update() as each time domain sample arrives. Output
+data is an array of `std::complex` values in the `dft` field. The length of this
+array is the length of the DFT.
+
+The output data is not valid until at least *N* samples have been processed. You
+can detect this using the `is_data_valid()` method, or by storing the return
+value of the `update()` method.
+
+This is a header-only C++ library. Simply copy sliding_dft.hpp into your
+project. The included CMakeLists.txt is for building the testbench.
+
+Implementation details
+----
+
+See references [1, 2] for an overview of sliding DFT algorithms. A damping
+factor is used to improve stability in the face of numerical rounding errors. If
+you experience stability issues, reduce `dft.damping_factor`. It should be
+slightly less than one.
+
+Windowing is done using a Hanning window, computed in the frequency domain [1].
+
+[1] E. Jacobsen and R. Lyons, “The Sliding DFT,” IEEE Signal Process. Mag., vol. 20, no. 2, pp. 74–80, Mar. 2003.
+
+[2] E. Jacobsen and R. Lyons, “An Update to the Sliding DFT,” IEEE Signal Process. Mag., vol. 21, no. 1, pp. 110-111, 2004.
+
+
+MIT License
+----
+
+Copyright (c) 2016 Bronson Philippa
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE 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
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, 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 IN THE
+SOFTWARE.

sliding_dft.hpp

@@ -0,0 +1,167 @@
+/**
+Sliding discrete Fourier transform (C++)
+====
+
+This code efficiently computes discrete Fourier transforms (DFTs) from a
+continuous sequence of input values. It is a recursive algorithm that updates
+the DFT when each new time-domain measurement arrives, effectively applying a
+sliding window over the last *N* samples. This implementation applies the
+Hanning window in order to minimise spectral leakage.
+
+The update step has computational complexity *O(N)*. If a new DFT is required
+every *M* samples, and *M* < log2(*N*), then this approach is more efficient
+that recalculating the DFT from scratch each time.
+
+This is a header-only C++ library. Simply copy sliding_dft.hpp into your
+project, and use it as follows:
+
+	// Use double precision arithmetic and a 512-length DFT
+	static SlidingDFT<double, 512> dft;
+	// avoid allocating on the stack because the object is large
+
+	// When a new time sample arrives, update the DFT with:
+	dft.update(x);
+
+	// After at least 512 samples have been processed:
+	std::complex<double> DC_bin = dft.dft[0];
+
+Your application should call update() as each time domain sample arrives. Output
+data is an array of `std::complex` values in the `dft` field. The length of this
+array is the length of the DFT.
+
+The output data is not valid until at least *N* samples have been processed. You
+can detect this using the `is_data_valid()` method, or by storing the return
+value of the `update()` method.
+
+This is a header-only C++ library. Simply copy sliding_dft.hpp into your
+project. The included CMakeLists.txt is for building the testbench.
+
+Implementation details
+----
+
+See references [1, 2] for an overview of sliding DFT algorithms. A damping
+factor is used to improve stability in the face of numerical rounding errors. If
+you experience stability issues, reduce `dft.damping_factor`. It should be
+slightly less than one.
+
+Windowing is done using a Hanning window, computed in the frequency domain [1].
+
+[1] E. Jacobsen and R. Lyons, “The Sliding DFT,” IEEE Signal Process. Mag., vol. 20, no. 2, pp. 74–80, Mar. 2003.
+
+[2] E. Jacobsen and R. Lyons, “An Update to the Sliding DFT,” IEEE Signal Process. Mag., vol. 21, no. 1, pp. 110-111, 2004.
+
+
+MIT License
+----
+
+Copyright (c) 2016 Bronson Philippa
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE 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
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, 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 IN THE
+SOFTWARE.
+*/
+
+#pragma once
+
+#define _USE_MATH_DEFINES
+#include <complex>
+#include <math.h>
+
+template <class NumberFormat, size_t DFT_Length>
+class SlidingDFT
+{
+private:
+	/// Are the frequency domain values valid? (i.e. have at elast DFT_Length data
+	/// points been seen?)
+	bool data_valid = false;
+
+	/// Time domain samples are stored in this circular buffer.
+	NumberFormat x[DFT_Length] = { 0 };
+
+	/// Index of the next item in the buffer to be used. Equivalently, the number
+	/// of samples that have been seen so far modulo DFT_Length.
+	size_t x_index = 0;
+
+	/// Twiddle factors for the update algorithm
+	std::complex<NumberFormat> twiddle[DFT_Length];
+
+	/// Frequency domain values (unwindowed!)
+	std::complex<NumberFormat> S[DFT_Length];
+
+public:
+	/// Frequency domain values (windowed)
+	std::complex<NumberFormat> dft[DFT_Length];
+
+	/// A damping factor introduced into the recursive DFT algorithm to guarantee
+	/// stability.
+	NumberFormat damping_factor = std::nexttoward((NumberFormat)1, (NumberFormat)0);
+
+	/// Constructor
+	SlidingDFT()
+	{
+		const std::complex<NumberFormat> j(0.0, 1.0);
+		const NumberFormat N = DFT_Length;
+
+		// Compute the twiddle factors, and zero the x and S arrays
+		for (size_t k = 0; k < DFT_Length; k++) {
+			NumberFormat factor = (NumberFormat)(2.0 * M_PI) * k / N;
+			this->twiddle[k] = std::exp(j * factor);
+			this->S[k] = 0;
+			this->x[k] = 0;
+		}
+	}
+
+	/// Determine whether the output data is valid
+	bool is_data_valid()
+	{
+		return this->data_valid;
+	}
+
+	/// Update the calculation with a new sample
+	/// Returns true if the data are valid (because enough samples have been
+	/// presented), or false if the data are invalid.
+	bool update(NumberFormat new_x)
+	{
+		// Update the storage of the time domain values
+		const NumberFormat old_x = this->x[this->x_index];
+		this->x[this->x_index] = new_x;
+
+		// Update the DFT
+		const NumberFormat r = this->damping_factor;
+		const NumberFormat r_to_N = pow(r, (NumberFormat)DFT_Length);
+		for (size_t k = 0; k < DFT_Length; k++) {
+			this->S[k] = this->twiddle[k] * (r * this->S[k] - r_to_N * old_x + new_x);
+		}
+
+		// Apply the Hanning window
+		this->dft[0] = (NumberFormat)0.5*this->S[0] - (NumberFormat)0.25*(this->S[DFT_Length - 1] + this->S[1]);
+		for (size_t k = 1; k < (DFT_Length - 1); k++) {
+			this->dft[k] = (NumberFormat)0.5*this->S[k] - (NumberFormat)0.25*(this->S[k - 1] + this->S[k + 1]);
+		}
+		this->dft[DFT_Length - 1] = (NumberFormat)0.5*this->S[DFT_Length - 1] - (NumberFormat)0.25*(this->S[DFT_Length - 2] + this->S[0]);
+
+		// Increment the counter
+		this->x_index++;
+		if (this->x_index >= DFT_Length) {
+			this->data_valid = true;
+			this->x_index = 0;
+		}
+
+		// Done.
+		return this->data_valid;
+	}
+};

testbench/generate_test1_data.m

@@ -0,0 +1,85 @@
+%% Generate test data
+FFT_Length = 64;
+t = (0:(1/100):1.2)';
+% FFT_Length = 8;
+% t = linspace(0, 10, 16)';
+v = 0.5*randn(size(t)) + 2*sin(2*pi*(0.1+t/0.07).*t);
+plot(t, v)
+
+% Calculate the FFT
+s = zeros(numel(v) - FFT_Length, FFT_Length);
+window = hanning(FFT_Length, 'periodic');
+num_FFTs = numel(v) - FFT_Length + 1;
+for i = 1:num_FFTs
+   s(i, :) = fft(window .* v(i+(1:FFT_Length)-1));
+end
+
+% Plot the amplitude
+figure()
+surf(1:num_FFTs, 1:FFT_Length, (abs(s)'), 'LineStyle', 'none')
+xlabel('time');
+ylabel('freq');
+
+%% Write the data file
+fid = fopen('test1_data.h', 'wt');
+% fid = 1;
+
+fprintf(fid, 'static const size_t test1_FFT_LENGTH = %i;\n', FFT_Length);
+fprintf(fid, 'static const size_t test1_TIME_LENGTH = %i;\n\n', numel(t));
+
+fprintf(fid, 'static const unsigned long long test1_time_domain_data [test1_TIME_LENGTH] = {\n    ');
+for i = 1:numel(v)
+    fprintf(fid, '0x%s', num2hex(v(i)));
+    if i < numel(v)
+        fprintf(fid, ',');
+    end
+    if mod(i-1, 8) == 7
+        fprintf(fid, '\n    ');
+    end
+end
+fprintf(fid, '};\n');
+fprintf(fid, 'static const double *test1_time_domain = (const double *)test1_time_domain_data;\n\n');
+
+fprintf(fid, 'static const unsigned long long test1_dft_real_data [test1_TIME_LENGTH-test1_FFT_LENGTH+1][test1_FFT_LENGTH] = {\n');
+for i = 1:num_FFTs
+    fprintf(fid, '    { // time index = %i\n        ', i-1);
+    for j = 1:FFT_Length
+        fprintf(fid, '0x%s', num2hex(real(s(i, j))));
+        if j < FFT_Length
+            fprintf(fid, ',');
+            if mod(j-1, 8) == 7
+                fprintf(fid, '\n        ');
+            end
+        end
+    end
+    fprintf(fid, '\n    }');
+    if i < num_FFTs
+        fprintf(fid, ',');
+    end
+    fprintf(fid, '\n');
+end
+fprintf(fid, '};\n');
+fprintf(fid, 'static const double (*test1_dft_real)[test1_FFT_LENGTH] = (const double (*)[test1_FFT_LENGTH])test1_dft_real_data;\n\n');
+
+fprintf(fid, 'static const unsigned long long test1_dft_imag_data [test1_TIME_LENGTH-test1_FFT_LENGTH+1][test1_FFT_LENGTH] = {\n');
+for i = 1:num_FFTs
+    fprintf(fid, '    { // time index = %i\n        ', i-1);
+    for j = 1:FFT_Length
+        fprintf(fid, '0x%s', num2hex(imag(s(i, j))));
+        if j < FFT_Length
+            fprintf(fid, ',');
+            if mod(j-1, 8) == 7
+                fprintf(fid, '\n        ');
+            end
+        end
+    end
+    fprintf(fid, '\n    }');
+    if i < num_FFTs
+        fprintf(fid, ',');
+    end
+    fprintf(fid, '\n');
+end
+fprintf(fid, '};\n');
+fprintf(fid, 'static const double (*test1_dft_imag)[test1_FFT_LENGTH] = (const double (*)[test1_FFT_LENGTH])test1_dft_imag_data;\n\n');
+
+fclose(fid);

testbench/main.cpp

@@ -0,0 +1,43 @@
+#include <iostream>
+#include <complex>
+
+#include "sliding_dft.hpp"
+#include "testbench/test1_data.h"
+
+using std::complex;
+
+template<class T>
+void test(const char *name)
+{
+	static SlidingDFT<T, test1_FFT_LENGTH> dft;
+	double error_sum = 0;
+	size_t num_comparisons = 0;
+
+	for (size_t i = 0; i < test1_TIME_LENGTH; i++) {
+		dft.update((T)test1_time_domain[i]);
+		if (dft.is_data_valid()) {
+			for (size_t j = 0; j < test1_FFT_LENGTH; j++) {
+				complex<double> calculated = dft.dft[j];
+				complex<double> exact(test1_dft_real[i - test1_FFT_LENGTH + 1][j], test1_dft_imag[i - test1_FFT_LENGTH + 1][j]);
+				double error = abs(calculated - exact);
+				/*printf("[%2i,%2i] Calculated (%10.3e,%10.3e); expected (%10.3e,%10.3e), error=%11.6e\n", i, j,
+					   calculated.real(), calculated.imag(),
+					   exact.real(), exact.imag(),
+					   error);*/
+				error_sum += error;
+				num_comparisons++;
+			}
+		}
+	}
+	double mean_error = error_sum / num_comparisons;
+	printf("%s: compared %zu complex numbers with an average error of %e\n", name, num_comparisons, mean_error);
+}
+
+
+int main(int argc, const char* argv[])
+{
+	test<float>("Single precision");
+	test<double>("Double precision");
+
+	return 0;
+}