Move DFT/TDR calculation into dedicated thread, limit update rate

This commit is contained in:
Jan Käberich 2021-12-30 15:23:07 +01:00
parent fd786c4176
commit 8e47d14192
5 changed files with 299 additions and 148 deletions

View File

@ -6,6 +6,8 @@
#include "ui_dftdialog.h" #include "ui_dftdialog.h"
#include "ui_dftexplanationwidget.h" #include "ui_dftexplanationwidget.h"
#include <QDebug>
using namespace std; using namespace std;
Math::DFT::DFT() Math::DFT::DFT()
@ -13,9 +15,22 @@ Math::DFT::DFT()
automaticDC = true; automaticDC = true;
DCfreq = 1000000000.0; DCfreq = 1000000000.0;
destructing = false;
thread = new DFTThread(*this);
thread->start(TDRThread::Priority::LowestPriority);
connect(&window, &WindowFunction::changed, this, &DFT::updateDFT); connect(&window, &WindowFunction::changed, this, &DFT::updateDFT);
} }
Math::DFT::~DFT()
{
// tell thread to exit
destructing = true;
semphr.release();
thread->wait();
delete thread;
}
TraceMath::DataType Math::DFT::outputType(TraceMath::DataType inputType) TraceMath::DataType Math::DFT::outputType(TraceMath::DataType inputType)
{ {
if(inputType == DataType::Time) { if(inputType == DataType::Time) {
@ -120,16 +135,47 @@ void Math::DFT::inputSamplesChanged(unsigned int begin, unsigned int end)
// not the end, do nothing // not the end, do nothing
return; return;
} }
double DC = DCfreq; // trigger calculation in thread
semphr.release();
success();
}
void Math::DFT::updateDFT()
{
if(dataType != DataType::Invalid) {
inputSamplesChanged(0, input->rData().size());
}
}
Math::DFTThread::DFTThread(Math::DFT &dft)
: dft(dft)
{
}
void Math::DFTThread::run()
{
qDebug() << "DFT thread starting";
while(1) {
dft.semphr.acquire();
// clear possible additional semaphores
dft.semphr.tryAcquire(dft.semphr.available());
if(dft.destructing) {
// TDR object about to be deleted, exit thread
qDebug() << "DFT thread exiting";
return;
}
qDebug() << "DFT thread calculating";
double DC = dft.DCfreq;
TDR *tdr = nullptr; TDR *tdr = nullptr;
if(automaticDC) { if(dft.automaticDC) {
// find the last operation that transformed from the frequency domain to the time domain // find the last operation that transformed from the frequency domain to the time domain
auto in = input; auto in = dft.input;
while(in->getInput()->getDataType() != DataType::Frequency) { while(in->getInput()->getDataType() != DFT::DataType::Frequency) {
in = input->getInput(); in = dft.input->getInput();
} }
switch(in->getType()) { switch(in->getType()) {
case Type::TDR: { case DFT::Type::TDR: {
tdr = static_cast<TDR*>(in); tdr = static_cast<TDR*>(in);
if(tdr->getMode() == TDR::Mode::Lowpass) { if(tdr->getMode() == TDR::Mode::Lowpass) {
DC = 0; DC = 0;
@ -145,28 +191,28 @@ void Math::DFT::inputSamplesChanged(unsigned int begin, unsigned int end)
break; break;
} }
} }
auto samples = input->rData().size(); auto samples = dft.input->rData().size();
auto timeSpacing = input->rData()[1].x - input->rData()[0].x; auto timeSpacing = dft.input->rData()[1].x - dft.input->rData()[0].x;
vector<complex<double>> timeDomain(samples); vector<complex<double>> timeDomain(samples);
for(unsigned int i=0;i<samples;i++) { for(unsigned int i=0;i<samples;i++) {
timeDomain.at(i) = input->rData()[i].y; timeDomain.at(i) = dft.input->rData()[i].y;
} }
Fft::shift(timeDomain, false); Fft::shift(timeDomain, false);
window.apply(timeDomain); dft.window.apply(timeDomain);
Fft::shift(timeDomain, true); Fft::shift(timeDomain, true);
Fft::transform(timeDomain, false); Fft::transform(timeDomain, false);
// shift DC bin into the middle // shift DC bin into the middle
Fft::shift(timeDomain, false); Fft::shift(timeDomain, false);
double binSpacing = 1.0 / (timeSpacing * timeDomain.size()); double binSpacing = 1.0 / (timeSpacing * timeDomain.size());
data.clear(); dft.data.clear();
int DCbin = timeDomain.size() / 2, startBin = 0; int DCbin = timeDomain.size() / 2, startBin = 0;
if(DC > 0) { if(DC > 0) {
data.resize(timeDomain.size()); dft.data.resize(timeDomain.size());
} else { } else {
startBin = (timeDomain.size()+1) / 2; startBin = (timeDomain.size()+1) / 2;
data.resize(timeDomain.size()/2); dft.data.resize(timeDomain.size()/2);
} }
// reverse effect of frequency domain window function from TDR (if available) // reverse effect of frequency domain window function from TDR (if available)
@ -176,16 +222,9 @@ void Math::DFT::inputSamplesChanged(unsigned int begin, unsigned int end)
for(int i = startBin;(unsigned int) i<timeDomain.size();i++) { for(int i = startBin;(unsigned int) i<timeDomain.size();i++) {
auto freq = (i - DCbin) * binSpacing + DC; auto freq = (i - DCbin) * binSpacing + DC;
data[i - startBin].x = round(freq); dft.data[i - startBin].x = round(freq);
data[i - startBin].y = timeDomain.at(i); dft.data[i - startBin].y = timeDomain.at(i);
} }
emit outputSamplesChanged(0, data.size()); emit dft.outputSamplesChanged(0, dft.data.size());
success();
}
void Math::DFT::updateDFT()
{
if(dataType != DataType::Invalid) {
inputSamplesChanged(0, input->rData().size());
} }
} }

View File

@ -4,12 +4,31 @@
#include "tracemath.h" #include "tracemath.h"
#include "windowfunction.h" #include "windowfunction.h"
#include <QThread>
#include <QSemaphore>
namespace Math { namespace Math {
class DFT;
class DFTThread : public QThread
{
Q_OBJECT
public:
DFTThread(DFT &dft);
~DFTThread(){};
private:
void run() override;
DFT &dft;
};
class DFT : public TraceMath class DFT : public TraceMath
{ {
friend class DFTThread;
Q_OBJECT
public: public:
DFT(); DFT();
~DFT();
DataType outputType(DataType inputType) override; DataType outputType(DataType inputType) override;
QString description() override; QString description() override;
@ -29,6 +48,9 @@ private:
bool automaticDC; bool automaticDC;
double DCfreq; double DCfreq;
WindowFunction window; WindowFunction window;
DFTThread *thread;
bool destructing;
QSemaphore semphr;
}; };
} }

View File

@ -19,9 +19,22 @@ TDR::TDR()
stepResponse = true; stepResponse = true;
mode = Mode::Lowpass; mode = Mode::Lowpass;
destructing = false;
thread = new TDRThread(*this);
thread->start(TDRThread::Priority::LowestPriority);
connect(&window, &WindowFunction::changed, this, &TDR::updateTDR); connect(&window, &WindowFunction::changed, this, &TDR::updateTDR);
} }
TDR::~TDR()
{
// tell thread to exit
destructing = true;
semphr.release();
thread->wait();
delete thread;
}
TraceMath::DataType TDR::outputType(TraceMath::DataType inputType) TraceMath::DataType TDR::outputType(TraceMath::DataType inputType)
{ {
if(inputType == DataType::Frequency) { if(inputType == DataType::Frequency) {
@ -190,85 +203,8 @@ void TDR::inputSamplesChanged(unsigned int begin, unsigned int end)
// not the end, do nothing // not the end, do nothing
return; return;
} }
vector<complex<double>> frequencyDomain; // trigger calculation in thread
auto stepSize = (input->rData().back().x - input->rData().front().x) / (input->rData().size() - 1); semphr.release();
if(mode == Mode::Lowpass) {
if(stepResponse) {
auto steps = input->rData().size();
auto firstStep = input->rData().front().x;
// frequency points need to be evenly spaced all the way to DC
if(firstStep == 0) {
// zero as first step would result in infinite number of points, skip and start with second
firstStep = input->rData()[1].x;
steps--;
}
if(firstStep * steps != input->rData().back().x) {
// data is not available with correct frequency spacing, calculate required steps
steps = input->rData().back().x / firstStep;
stepSize = firstStep;
}
frequencyDomain.resize(2 * steps + 1);
// copy frequencies, use the flipped conjugate for negative part
for(unsigned int i = 1;i<=steps;i++) {
auto S = input->getInterpolatedSample(stepSize * i).y;
frequencyDomain[steps - i] = conj(S);
frequencyDomain[steps + i] = S;
}
if(automaticDC) {
// use simple extrapolation from lowest two points to extract DC value
auto abs_DC = 2.0 * abs(frequencyDomain[steps + 1]) - abs(frequencyDomain[steps + 2]);
auto phase_DC = 2.0 * arg(frequencyDomain[steps + 1]) - arg(frequencyDomain[steps + 2]);
frequencyDomain[steps] = polar(abs_DC, phase_DC);
} else {
frequencyDomain[steps] = manualDC;
}
} else {
auto steps = input->rData().size();
unsigned int offset = 0;
if(input->rData().front().x == 0) {
// DC measurement is inaccurate, skip
steps--;
offset++;
}
// no step response required, can use frequency values as they are. No extra extrapolated DC value here -> 2 values less than with step response
frequencyDomain.resize(2 * steps - 1);
frequencyDomain[steps - 1] = input->rData()[offset].y;
for(unsigned int i = 1;i<steps;i++) {
auto S = input->rData()[i + offset].y;
frequencyDomain[steps - i - 1] = conj(S);
frequencyDomain[steps + i - 1] = S;
}
}
} else {
// bandpass mode
// Can use input data directly, no need to extend with complex conjugate
frequencyDomain.resize(input->rData().size());
for(unsigned int i=0;i<input->rData().size();i++) {
frequencyDomain[i] = input->rData()[i].y;
}
}
window.apply(frequencyDomain);
Fft::shift(frequencyDomain, true);
auto fft_bins = frequencyDomain.size();
const double fs = 1.0 / (stepSize * fft_bins);
Fft::transform(frequencyDomain, true);
data.clear();
data.resize(fft_bins);
for(unsigned int i = 0;i<fft_bins;i++) {
data[i].x = fs * i;
data[i].y = frequencyDomain[i] / (double) fft_bins;
}
if(stepResponse && mode == Mode::Lowpass) {
updateStepResponse(true);
} else {
updateStepResponse(false);
}
emit outputSamplesChanged(0, data.size());
success(); success();
} else { } else {
// not enough input data // not enough input data
@ -295,3 +231,105 @@ TDR::Mode TDR::getMode() const
{ {
return mode; return mode;
} }
TDRThread::TDRThread(TDR &tdr)
: tdr(tdr)
{
}
void TDRThread::run()
{
qDebug() << "TDR thread starting";
while(1) {
tdr.semphr.acquire();
// clear possible additional semaphores
tdr.semphr.tryAcquire(tdr.semphr.available());
if(tdr.destructing) {
// TDR object about to be deleted, exit thread
qDebug() << "TDR thread exiting";
return;
}
qDebug() << "TDR thread calculating";
// perform calculation
vector<complex<double>> frequencyDomain;
auto stepSize = (tdr.input->rData().back().x - tdr.input->rData().front().x) / (tdr.input->rData().size() - 1);
if(tdr.mode == TDR::Mode::Lowpass) {
if(tdr.stepResponse) {
auto steps = tdr.input->rData().size();
auto firstStep = tdr.input->rData().front().x;
// frequency points need to be evenly spaced all the way to DC
if(firstStep == 0) {
// zero as first step would result in infinite number of points, skip and start with second
firstStep = tdr.input->rData()[1].x;
steps--;
}
if(firstStep * steps != tdr.input->rData().back().x) {
// data is not available with correct frequency spacing, calculate required steps
steps = tdr.input->rData().back().x / firstStep;
stepSize = firstStep;
}
frequencyDomain.resize(2 * steps + 1);
// copy frequencies, use the flipped conjugate for negative part
for(unsigned int i = 1;i<=steps;i++) {
auto S = tdr.input->getInterpolatedSample(stepSize * i).y;
frequencyDomain[steps - i] = conj(S);
frequencyDomain[steps + i] = S;
}
if(tdr.automaticDC) {
// use simple extrapolation from lowest two points to extract DC value
auto abs_DC = 2.0 * abs(frequencyDomain[steps + 1]) - abs(frequencyDomain[steps + 2]);
auto phase_DC = 2.0 * arg(frequencyDomain[steps + 1]) - arg(frequencyDomain[steps + 2]);
frequencyDomain[steps] = polar(abs_DC, phase_DC);
} else {
frequencyDomain[steps] = tdr.manualDC;
}
} else {
auto steps = tdr.input->rData().size();
unsigned int offset = 0;
if(tdr.input->rData().front().x == 0) {
// DC measurement is inaccurate, skip
steps--;
offset++;
}
// no step response required, can use frequency values as they are. No extra extrapolated DC value here -> 2 values less than with step response
frequencyDomain.resize(2 * steps - 1);
frequencyDomain[steps - 1] = tdr.input->rData()[offset].y;
for(unsigned int i = 1;i<steps;i++) {
auto S = tdr.input->rData()[i + offset].y;
frequencyDomain[steps - i - 1] = conj(S);
frequencyDomain[steps + i - 1] = S;
}
}
} else {
// bandpass mode
// Can use input data directly, no need to extend with complex conjugate
frequencyDomain.resize(tdr.input->rData().size());
for(unsigned int i=0;i<tdr.input->rData().size();i++) {
frequencyDomain[i] = tdr.input->rData()[i].y;
}
}
tdr.window.apply(frequencyDomain);
Fft::shift(frequencyDomain, true);
auto fft_bins = frequencyDomain.size();
const double fs = 1.0 / (stepSize * fft_bins);
Fft::transform(frequencyDomain, true);
tdr.data.clear();
tdr.data.resize(fft_bins);
for(unsigned int i = 0;i<fft_bins;i++) {
tdr.data[i].x = fs * i;
tdr.data[i].y = frequencyDomain[i] / (double) fft_bins;
}
if(tdr.stepResponse && tdr.mode == TDR::Mode::Lowpass) {
tdr.updateStepResponse(true);
} else {
tdr.updateStepResponse(false);
}
emit tdr.outputSamplesChanged(0, tdr.data.size());
}
}

View File

@ -4,12 +4,31 @@
#include "tracemath.h" #include "tracemath.h"
#include "windowfunction.h" #include "windowfunction.h"
#include <QThread>
#include <QSemaphore>
namespace Math { namespace Math {
class TDR;
class TDRThread : public QThread
{
Q_OBJECT
public:
TDRThread(TDR &tdr);
~TDRThread(){};
private:
void run() override;
TDR &tdr;
};
class TDR : public TraceMath class TDR : public TraceMath
{ {
friend class TDRThread;
Q_OBJECT
public: public:
TDR(); TDR();
~TDR();
DataType outputType(DataType inputType) override; DataType outputType(DataType inputType) override;
QString description() override; QString description() override;
@ -39,6 +58,9 @@ private:
bool stepResponse; bool stepResponse;
bool automaticDC; bool automaticDC;
std::complex<double> manualDC; std::complex<double> manualDC;
TDRThread *thread;
bool destructing;
QSemaphore semphr;
}; };
} }

View File

@ -346,13 +346,35 @@ void Math::TimeGateGraph::paintEvent(QPaintEvent *event)
p.setBackground(QBrush(pref.Graphs.Color.background)); p.setBackground(QBrush(pref.Graphs.Color.background));
p.fillRect(0, 0, width(), height(), QBrush(pref.Graphs.Color.background)); p.fillRect(0, 0, width(), height(), QBrush(pref.Graphs.Color.background));
if(!gate->getInput() || !gate->getInput()->rData().size()) {
// no data yet, nothing to plot
return;
}
// plot trace // plot trace
auto pen = QPen(Qt::green, 1); auto pen = QPen(Qt::green, 1);
pen.setCosmetic(true); pen.setCosmetic(true);
pen.setStyle(Qt::SolidLine); pen.setStyle(Qt::SolidLine);
p.setPen(pen); p.setPen(pen);
for(unsigned int i=1;i<input.size();i++) {
auto last = input[i-1]; auto minX = input.front().x;
auto maxX = input.back().x;
int plotLeft = 0;
int plotRight = size().width();
int plotTop = 0;
int plotBottom = size().height();
QPoint p1, p2;
// limit amount of displayed points to keep GUI snappy
auto increment = input.size() / 500;
if(!increment) {
increment = 1;
}
for(unsigned int i=increment;i<input.size();i+=increment) {
auto last = input[i-increment];
auto now = input[i]; auto now = input[i];
auto y_last = Util::SparamTodB(last.y); auto y_last = Util::SparamTodB(last.y);
@ -363,8 +385,12 @@ void Math::TimeGateGraph::paintEvent(QPaintEvent *event)
} }
// scale to plot coordinates // scale to plot coordinates
auto p1 = plotValueToPixel(last.x, y_last); p1.setX(Util::Scale<double>(last.x, minX, maxX, plotLeft, plotRight));
auto p2 = plotValueToPixel(now.x, y_now); p1.setY(Util::Scale<double>(y_last, -120, 20, plotBottom, plotTop));
p2.setX(Util::Scale<double>(now.x, minX, maxX, plotLeft, plotRight));
p2.setY(Util::Scale<double>(y_now, -120, 20, plotBottom, plotTop));
// auto p1 = plotValueToPixel(last.x, y_last);
// auto p2 = plotValueToPixel(now.x, y_now);
// draw line // draw line
p.drawLine(p1, p2); p.drawLine(p1, p2);
} }
@ -372,11 +398,11 @@ void Math::TimeGateGraph::paintEvent(QPaintEvent *event)
auto filter = gate->rFilter(); auto filter = gate->rFilter();
pen = QPen(Qt::red, 1); pen = QPen(Qt::red, 1);
p.setPen(pen); p.setPen(pen);
for(unsigned int i=1;i<filter.size();i++) { for(unsigned int i=increment;i<filter.size();i+=increment) {
auto x_last = input[i-1].x; auto x_last = input[i-increment].x;
auto x_now = input[i].x; auto x_now = input[i].x;
auto f_last = Util::SparamTodB(filter[i-1]); auto f_last = Util::SparamTodB(filter[i-increment]);
auto f_now = Util::SparamTodB(filter[i]); auto f_now = Util::SparamTodB(filter[i]);
if(std::isnan(f_last) || std::isnan(f_now) || std::isinf(f_last) || std::isinf(f_now)) { if(std::isnan(f_last) || std::isnan(f_now) || std::isinf(f_last) || std::isinf(f_now)) {
@ -384,8 +410,12 @@ void Math::TimeGateGraph::paintEvent(QPaintEvent *event)
} }
// scale to plot coordinates // scale to plot coordinates
auto p1 = plotValueToPixel(x_last, f_last); p1.setX(Util::Scale<double>(x_last, minX, maxX, plotLeft, plotRight));
auto p2 = plotValueToPixel(x_now, f_now); p1.setY(Util::Scale<double>(f_last, -120, 20, plotBottom, plotTop));
p2.setX(Util::Scale<double>(x_now, minX, maxX, plotLeft, plotRight));
p2.setY(Util::Scale<double>(f_now, -120, 20, plotBottom, plotTop));
// auto p1 = plotValueToPixel(x_last, f_last);
// auto p2 = plotValueToPixel(x_now, f_now);
// draw line // draw line
p.drawLine(p1, p2); p.drawLine(p1, p2);