feat : adding asynchronous sliding DFT (need to be tested)

Signed-off-by: Vincenzo Petrolo <vincenzo@kernel-space.org>

.gitmodules

@@ -0,0 +1,6 @@
+[submodule "SlidingDFT"]
+	path = SlidingDFT
+	url = https://github.com/bronsonp/SlidingDFT.git
+[submodule "Desktop_Interface/SlidingDFT"]
+	path = Desktop_Interface/SlidingDFT
+	url = https://github.com/bronsonp/SlidingDFT.git

Desktop_Interface/SlidingDFT

@@ -0,0 +1 @@
+Subproject commit 13e0121253c73f2e85e418a3ae960463c0fbf31f

Desktop_Interface/asyncdft.cpp

@@ -0,0 +1,89 @@
+#include "asyncdft.h"
+#include <iostream>
+
+#define DBG 1
+
+AsyncDFT::AsyncDFT()
+{
+    /*Creating the main thread, which will manage everything*/
+    stopping = false;
+    manager = std::thread(&AsyncDFT::threadManager, this);
+}
+
+AsyncDFT::~AsyncDFT()
+{
+    stopping = true;
+    dft_mutex.unlock();
+    while (!manager.joinable());
+    manager.join();
+#if DBG
+    std::cout << "Joined manager thread [DESTRUCTOR]" << std::endl;
+#endif
+}
+
+void AsyncDFT::threadManager()
+{
+    while(stopping == false) {
+            dft_mutex.lock();
+            std::thread t = std::thread(&AsyncDFT::updateDFT, this);
+#if DBG
+            std::cout << "Starting new thread [MANAGER]" << std::endl;
+#endif
+            while(!t.joinable());
+            t.join();
+#if DBG
+            std::cout << "Joined a thread [MANAGER]" << std::endl;
+#endif
+    }
+}
+
+void AsyncDFT::updateDFT()
+{
+#if DBG
+    std::cout << "Im a thread and i gained lock, now calculating [THREAD]" << std::endl;
+#endif
+    if (!waiting_jobs.empty()) {
+        short sample = waiting_jobs.front();
+        waiting_jobs.pop();
+        dft.update(sample);
+    }
+#if DBG
+    std::cout << "Im a thread and i released lock, exiting [THREAD]" << std::endl;
+#endif
+}
+
+void AsyncDFT::addSample(short sample)
+{
+    /*Adding to the waiting jobs the sample*/
+    waiting_jobs.push(sample);
+    dft_mutex.unlock();
+}
+
+QVector<double> AsyncDFT::getPowerSpectrum()
+{
+    /*Before doing anything, check if sliding DFT is computable*/
+    if (!dft.is_data_valid()) {
+        throw std::exception();
+    }
+
+    QVector<double> pow_spectr(n_samples, 0);
+    double maximum = -1;
+
+    /*Computing Power Spectrum*/
+    for (int i = 0; i < n_samples; i++) {
+        pow_spectr[i] = sqrt(dft.dft[i].real()*dft.dft[i].real() +
+                             dft.dft[i].imag()*dft.dft[i].imag());
+
+        if (pow_spectr[i] > maximum) {
+            maximum = pow_spectr[i];
+        }
+    }
+
+    for (int i = 0; i < n_samples; i++) {
+        pow_spectr[i] /= maximum;
+        pow_spectr[i] *= 100;
+    }
+
+    /*Returning normalized spectrum*/
+    return pow_spectr;
+}

Desktop_Interface/asyncdft.h

@@ -0,0 +1,36 @@
+#ifndef ASYNCDFT_H
+#define ASYNCDFT_H
+#include "sliding_dft.hpp"
+#include <thread>
+#include <QVector>
+#include <mutex>
+#include <queue>
+
+class AsyncDFT
+{
+public:
+    AsyncDFT();
+    ~AsyncDFT();
+    static const int n_samples = 350000;
+
+    /* Raise exception if not ready yet*/
+    QVector<double> getPowerSpectrum();
+
+    void addSample(short sample);
+
+private:
+    SlidingDFT<double, n_samples> dft;
+    std::thread dft_computing;
+    std::queue<short> waiting_jobs;
+
+    int pending_jobs;
+    std::mutex dft_mutex;
+
+    std::thread manager;
+    bool stopping;
+    void threadManager(); //threaded
+    void updateDFT(); //threaded
+
+};
+
+#endif // ASYNCDFT_H

Desktop_Interface/isobuffer.cpp

@@ -57,6 +57,8 @@ void isoBuffer::insertIntoBuffer(short item)
         m_back = 0;
     }
 
+    async_dft.addSample(item);
+
     checkTriggered();
 }
 
@@ -442,3 +444,28 @@ double isoBuffer::getDelayedTriggerPoint(double delay)
 
     return 0;
 }
+
+QVector<double> isoBuffer::getDFTPowerSpectrum()
+{
+    try {
+        return async_dft.getPowerSpectrum();
+    }  catch (std::exception) {
+        /*If DFT is not ready, returning a 0s array*/
+        std::cout << "DFT not yet ready " <<  std::endl;
+        return QVector<double>(async_dft.n_samples, 0);
+    }
+}
+
+QVector<double> isoBuffer::getFrequenciyWindow()
+{
+    int max_freq = 62500;
+    double delta_freq = ((double) 350000)/ ((double) async_dft.n_samples);
+    int tot = max_freq/delta_freq + 1;
+    QVector<double> f(tot);
+
+    for (int i = 0; i < tot; i++) {
+        f[i] = i*delta_freq;
+    }
+
+    return f;
+}

Desktop_Interface/isobuffer.h

@@ -16,6 +16,7 @@
 #include "xmega.h"
 #include "desktop_settings.h"
 #include "genericusbdriver.h"
+#include "asyncdft.h"
 
 class isoDriver;
 class uartStyleDecoder;
@@ -65,7 +66,6 @@ public:
 	void writeBuffer_short(short* data, int len);
 
 	std::unique_ptr<short[]> readBuffer(double sampleWindow, int numSamples, bool singleBit, double delayOffset);
-
 //	file I/O
 private:
 	void outputSampleToFile(double averageSample);
@@ -86,6 +86,10 @@ public:
     void setTriggerType(TriggerType newType);
     void setTriggerLevel(double voltageLevel, uint16_t top, bool acCoupled);
     double getDelayedTriggerPoint(double delay);
+    // DFT
+    AsyncDFT async_dft;
+    QVector<double> getDFTPowerSpectrum();
+    QVector<double> getFrequenciyWindow();
 
 // ---- MEMBER VARIABLES ----
 
@@ -129,6 +133,8 @@ private:
 	unsigned int m_currentColumn = 0;
 
 	isoDriver* m_virtualParent;
+    //DFT
+
 signals:
 	void fileIOinternalDisable();
 public slots:

Desktop_Interface/isodriver.cpp

@@ -4,6 +4,7 @@
 #include "platformspecific.h"
 #include <math.h>
 #include "daqloadprompt.h"
+#include "sliding_dft.hpp"
 #include <iostream>
 #include <omp.h>
 
@@ -34,19 +35,6 @@ isoDriver::isoDriver(QWidget *parent) : QLabel(parent)
     slowTimer->setTimerType(Qt::PreciseTimer);
     slowTimer->start(MULTIMETER_PERIOD);
     connect(slowTimer, SIGNAL(timeout()), this, SLOT(slowTimerTick()));
-
-    /*Creating DFT plan*/
-    fftw_init_threads();
-    fftw_plan_with_nthreads(omp_get_max_threads());
-    std::cout << "Starting with " << omp_get_max_threads() << "threads" << std::endl;
-    this->N = 1<<17;
-    this->N *= omp_get_max_threads();
-    this->in_buffer = fftw_alloc_real(N);
-    this->out_buffer = fftw_alloc_complex(N);
-    std::cout << in_buffer << " " << out_buffer << " " << N<<  std::endl;
-    this->plan = fftw_plan_dft_r2c_1d(N,in_buffer, out_buffer,0);
-    std::cout << plan <<  std::endl;
-
 }
 
 void isoDriver::setDriver(genericUsbDriver *newDriver){
@@ -625,64 +613,6 @@ void isoDriver::setTriggerMode(int newMode)
     triggerStateChanged();
 }
 
-QVector<double> isoDriver::getDFTAmplitude(QVector<double> input)
-{
-    /*Zero-padding for better resolution of DFT*/
-    QVector<double> amplitude(N/2+1,0);
-    static int count = 100;
-    double cur_maximum = -1;
-    for (int i = 0; i < input.length(); i++) {
-        in_buffer[i] = input[i];
-    }
-    fftw_execute(plan);
-    amplitude[0] = sqrt(out_buffer[0][0]*out_buffer[0][0] + out_buffer[0][1]*out_buffer[0][1]);  /* DC component */
-
-    cur_maximum = (amplitude[0] > cur_maximum ) ? amplitude[0] : cur_maximum;
-
-    for (int k = 1; k < (N+1)/2; ++k) {  /* (k < N/2 rounded up) */
-         amplitude[k] = sqrt(out_buffer[k][0]*out_buffer[k][0] + out_buffer[k][1]*out_buffer[k][1]);
-
-         cur_maximum = (amplitude[k] > cur_maximum ) ? amplitude[k] : cur_maximum;
-    }
-    if (N % 2 == 0) { /* N is even */
-         amplitude[N/2] = sqrt(out_buffer[N/2][0]*out_buffer[N/2][0]);  /* Nyquist freq. */
-
-         cur_maximum = (amplitude[N/2] > cur_maximum ) ? amplitude[N/2] : cur_maximum;
-
-    }
-
-    if (cur_maximum < maximum) {
-        count--;
-    } else {
-        /*current maximum is higher resetting count to 10*/
-        maximum = cur_maximum;
-        count = 100;
-    }
-
-    if (count == 0) {
-        /*current maximum is lower for 10 consecutive samples*/
-        maximum = cur_maximum;
-        count = 100;
-    }
-
-
-    return amplitude;
-}
-
-QVector<double> isoDriver::getFrequencies()
-{
-    int max_freq = 62500;
-    double delta_freq = ((double) 375000)/ ((double) N);
-    int tot = max_freq/delta_freq + 1;
-    QVector<double> f(tot);
-
-    for (int i = 0; i < tot; i++) {
-        f[i] = i*delta_freq;
-    }
-    return f;
-}
-
-
 //0 for off, 1 for ana, 2 for dig, -1 for ana750, -2 for file
 void isoDriver::frameActionGeneric(char CH1_mode, char CH2_mode)  
 {
@@ -739,19 +669,6 @@ void isoDriver::frameActionGeneric(char CH1_mode, char CH2_mode)
         if(CH2_mode) readData375_CH2 = internalBuffer375_CH2->readBuffer(display.window,GRAPH_SAMPLES,CH2_mode==2, display.delay + triggerDelay);
         if(CH1_mode == -1) readData750 = internalBuffer750->readBuffer(display.window,GRAPH_SAMPLES,false, display.delay + triggerDelay);
         if(CH1_mode == -2) readDataFile = internalBufferFile->readBuffer(display.window,GRAPH_SAMPLES,false, display.delay);
-    } else {
-        /*Don't allow moving frequency spectrum right or left
-         * by overwriting display window and delay before reading
-         * the buffer each time.
-         * @TODO improve this limitation.
-        */
-        double const_displ_window = ((double)N)/(internalBuffer375_CH1->m_samplesPerSecond);
-        double const_displ_delay = 0;
-        display.delay = const_displ_delay;
-        display.window = const_displ_window;
-        readData375_CH1 = internalBuffer375_CH1->readBuffer(const_displ_window,N,CH1_mode==2, const_displ_delay);
-        if(CH2_mode) readData375_CH2 = internalBuffer375_CH2->readBuffer(const_displ_window,N,CH2_mode==2,const_displ_delay);
-        if(CH1_mode == -1) readData750 = internalBuffer750->readBuffer(const_displ_window,N,false, const_displ_delay);
     }
 
     QVector<double> x(GRAPH_SAMPLES), CH1(GRAPH_SAMPLES), CH2(GRAPH_SAMPLES);
@@ -811,15 +728,22 @@ void isoDriver::frameActionGeneric(char CH1_mode, char CH2_mode)
         axes->yAxis->setRange(ymin, ymax);
     } else{
         if (spectrum) { /*If frequency spectrum mode*/
-            QVector<double> amplitude = getDFTAmplitude(CH1);
+            try {
-            QVector<double> f = getFrequencies();
+                QVector<double> amplitude = this->internalBuffer375_CH1->getDFTPowerSpectrum();
-            axes->graph(0)->setData(f,amplitude);
+                QVector<double> f = this->internalBuffer375_CH1->getFrequenciyWindow();
-            if(CH2_mode) {
+                axes->graph(0)->setData(f,amplitude);
-                amplitude = getDFTAmplitude(CH2);
+#if 0
-                axes->graph(1)->setData(f,amplitude);
+                if(CH2_mode) {
+                    amplitude = getDFTAmplitude(CH2);
+                    axes->graph(1)->setData(f,amplitude);
+                }
+#endif
+                axes->xAxis->setRange(f.length(), 0);
+                axes->yAxis->setRange(100,0);
+            }  catch (std::exception) {
+                std::cout << "Cannot yet get correct value for DFT" << std::endl;
             }
-            axes->xAxis->setRange(f.length(), 0);
+
-            axes->yAxis->setRange(maximum,0);
         } else {
             axes->graph(0)->setData(x,CH1);
             if(CH2_mode) axes->graph(1)->setData(x,CH2);

Desktop_Interface/isodriver.h

@@ -5,8 +5,6 @@
 #include <QLabel>
 #include <QDebug>
 #include <QVector>
-#include <fftw3.h>
-
 #include "qcustomplot.h"
 #include "genericusbdriver.h"
 #include "desktop_settings.h"
@@ -126,12 +124,6 @@ private:
     bool firstFrame = true;
     bool hexDisplay_CH1 = false;
     bool hexDisplay_CH2 = false;
-    //DFT
-    fftw_plan plan;
-    double *in_buffer;
-    fftw_complex *out_buffer;
-    int N;
-    double maximum = -1;
 
 
     //Generic Functions
@@ -195,9 +187,6 @@ private:
     uint8_t deviceMode_prev;
     //DAQ
     double daqLoad_startTime, daqLoad_endTime;
-    //DFT
-    QVector<double> getDFTAmplitude(QVector<double> input);
-    QVector<double> getFrequencies();
 
 signals:
     void setGain(double newGain);

Desktop_Interface/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;
+	}
+};

SlidingDFT

@@ -0,0 +1 @@
+Subproject commit 13e0121253c73f2e85e418a3ae960463c0fbf31f