diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 00000000..45fb05db --- /dev/null +++ b/.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 diff --git a/Desktop_Interface/SlidingDFT b/Desktop_Interface/SlidingDFT new file mode 160000 index 00000000..13e01212 --- /dev/null +++ b/Desktop_Interface/SlidingDFT @@ -0,0 +1 @@ +Subproject commit 13e0121253c73f2e85e418a3ae960463c0fbf31f diff --git a/Desktop_Interface/asyncdft.cpp b/Desktop_Interface/asyncdft.cpp new file mode 100644 index 00000000..f3e81dc6 --- /dev/null +++ b/Desktop_Interface/asyncdft.cpp @@ -0,0 +1,89 @@ +#include "asyncdft.h" +#include + +#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 AsyncDFT::getPowerSpectrum() +{ + /*Before doing anything, check if sliding DFT is computable*/ + if (!dft.is_data_valid()) { + throw std::exception(); + } + + QVector 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; +} diff --git a/Desktop_Interface/asyncdft.h b/Desktop_Interface/asyncdft.h new file mode 100644 index 00000000..5469777a --- /dev/null +++ b/Desktop_Interface/asyncdft.h @@ -0,0 +1,36 @@ +#ifndef ASYNCDFT_H +#define ASYNCDFT_H +#include "sliding_dft.hpp" +#include +#include +#include +#include + +class AsyncDFT +{ +public: + AsyncDFT(); + ~AsyncDFT(); + static const int n_samples = 350000; + + /* Raise exception if not ready yet*/ + QVector getPowerSpectrum(); + + void addSample(short sample); + +private: + SlidingDFT dft; + std::thread dft_computing; + std::queue waiting_jobs; + + int pending_jobs; + std::mutex dft_mutex; + + std::thread manager; + bool stopping; + void threadManager(); //threaded + void updateDFT(); //threaded + +}; + +#endif // ASYNCDFT_H diff --git a/Desktop_Interface/isobuffer.cpp b/Desktop_Interface/isobuffer.cpp index 5c986feb..cc532a91 100644 --- a/Desktop_Interface/isobuffer.cpp +++ b/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 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(async_dft.n_samples, 0); + } +} + +QVector isoBuffer::getFrequenciyWindow() +{ + int max_freq = 62500; + double delta_freq = ((double) 350000)/ ((double) async_dft.n_samples); + int tot = max_freq/delta_freq + 1; + QVector f(tot); + + for (int i = 0; i < tot; i++) { + f[i] = i*delta_freq; + } + + return f; +} diff --git a/Desktop_Interface/isobuffer.h b/Desktop_Interface/isobuffer.h index 07c506d2..52bf0b5a 100644 --- a/Desktop_Interface/isobuffer.h +++ b/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 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 getDFTPowerSpectrum(); + QVector getFrequenciyWindow(); // ---- MEMBER VARIABLES ---- @@ -129,6 +133,8 @@ private: unsigned int m_currentColumn = 0; isoDriver* m_virtualParent; + //DFT + signals: void fileIOinternalDisable(); public slots: diff --git a/Desktop_Interface/isodriver.cpp b/Desktop_Interface/isodriver.cpp index 8750b7f0..e167cacf 100644 --- a/Desktop_Interface/isodriver.cpp +++ b/Desktop_Interface/isodriver.cpp @@ -4,6 +4,7 @@ #include "platformspecific.h" #include #include "daqloadprompt.h" +#include "sliding_dft.hpp" #include #include @@ -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 isoDriver::getDFTAmplitude(QVector input) -{ - /*Zero-padding for better resolution of DFT*/ - QVector 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 isoDriver::getFrequencies() -{ - int max_freq = 62500; - double delta_freq = ((double) 375000)/ ((double) N); - int tot = max_freq/delta_freq + 1; - QVector 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 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 amplitude = getDFTAmplitude(CH1); - QVector f = getFrequencies(); - axes->graph(0)->setData(f,amplitude); - if(CH2_mode) { - amplitude = getDFTAmplitude(CH2); - axes->graph(1)->setData(f,amplitude); + try { + QVector amplitude = this->internalBuffer375_CH1->getDFTPowerSpectrum(); + QVector f = this->internalBuffer375_CH1->getFrequenciyWindow(); + axes->graph(0)->setData(f,amplitude); +#if 0 + 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); diff --git a/Desktop_Interface/isodriver.h b/Desktop_Interface/isodriver.h index 1f4df124..c6757b12 100644 --- a/Desktop_Interface/isodriver.h +++ b/Desktop_Interface/isodriver.h @@ -5,8 +5,6 @@ #include #include #include -#include - #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 getDFTAmplitude(QVector input); - QVector getFrequencies(); signals: void setGain(double newGain); diff --git a/Desktop_Interface/sliding_dft.hpp b/Desktop_Interface/sliding_dft.hpp new file mode 100644 index 00000000..63f30102 --- /dev/null +++ b/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 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 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 +#include + +template +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 twiddle[DFT_Length]; + + /// Frequency domain values (unwindowed!) + std::complex S[DFT_Length]; + +public: + /// Frequency domain values (windowed) + std::complex 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 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; + } +}; diff --git a/SlidingDFT b/SlidingDFT new file mode 160000 index 00000000..13e01212 --- /dev/null +++ b/SlidingDFT @@ -0,0 +1 @@ +Subproject commit 13e0121253c73f2e85e418a3ae960463c0fbf31f