mirror of
https://git.rwth-aachen.de/acs/public/villas/node/
synced 2025-03-09 00:00:00 +01:00
final fixes from Steffen
This commit is contained in:
Manuel Pitz
parent
700ec35392
commit
3dcdfe23d1
3 changed files with 155 additions and 167 deletions
|
|
@ -67,7 +67,7 @@ if(APPLE)
|
|||
endif()
|
||||
|
||||
add_definitions(-D_POSIX_C_SOURCE=200809L -D_GNU_SOURCE)
|
||||
add_compile_options(-Wall -Werror -fdiagnostics-color=auto)
|
||||
add_compile_options(-Wall -Wno-unknown-pragmas -Werror -fdiagnostics-color=auto)
|
||||
|
||||
# Check OS
|
||||
check_include_file("sys/eventfd.h" HAS_EVENTFD)
|
||||
|
|
|
|||
|
|
@ -34,9 +34,8 @@
|
|||
#include <villas/hook.hpp>
|
||||
#include <villas/path.h>
|
||||
#include <villas/sample.h>
|
||||
#include <villas/io.h>
|
||||
#include <villas/plugin.h>
|
||||
|
||||
/* Uncomment to enable dumper of memory windows */
|
||||
//#define DFT_MEM_DUMP
|
||||
|
||||
namespace villas {
|
||||
|
|
@ -71,54 +70,59 @@ protected:
|
|||
double frequency;
|
||||
double amplitude;
|
||||
double phase;
|
||||
double rocof;//rate of change of frequency
|
||||
double rocof; /**< Rate of change of frequency. */
|
||||
};
|
||||
|
||||
std::shared_ptr<Dumper> origSigSync;
|
||||
std::shared_ptr<Dumper> windowdSigSync;
|
||||
std::shared_ptr<Dumper> phasorPhase;
|
||||
std::shared_ptr<Dumper> phasorAmplitude;
|
||||
std::shared_ptr<Dumper> phasorFreq;
|
||||
std::shared_ptr<Dumper> ppsSigSync;
|
||||
|
||||
enum WindowType windowType;
|
||||
enum PaddingType paddingType;
|
||||
enum FreqEstimationType freqEstType;
|
||||
|
||||
struct format_type *format;
|
||||
|
||||
std::vector<std::vector<double>> smpMemory;
|
||||
std::vector<double> ppsMemory;//this is just temporary for debugging
|
||||
std::vector<std::vector<std::complex<double>>> dftMatrix;
|
||||
std::vector<std::vector<std::complex<double>>> dftResults;
|
||||
#ifdef DFT_MEM_DUMP
|
||||
std::vector<double> ppsMemory;
|
||||
#endif
|
||||
std::vector<std::vector<std::complex<double>>> matrix;
|
||||
std::vector<std::vector<std::complex<double>>> results;
|
||||
std::vector<double> filterWindowCoefficents;
|
||||
std::vector<std::vector<double>> absDftResults;
|
||||
std::vector<double> absDftFreqs;
|
||||
std::vector<std::vector<double>> absResults;
|
||||
std::vector<double> absFrequencies;
|
||||
|
||||
uint64_t dftCalcCount;
|
||||
uint64_t calcCount;
|
||||
unsigned sampleRate;
|
||||
double startFreqency;
|
||||
double startFrequency;
|
||||
double endFreqency;
|
||||
double frequencyResolution;
|
||||
unsigned dftRate;
|
||||
unsigned rate;
|
||||
unsigned ppsIndex;
|
||||
unsigned windowSize;
|
||||
unsigned windowMultiplier; /**< Multiplyer for the window to achieve frequency resolution */
|
||||
unsigned freqCount; /**< Number of requency bins that are calculated */
|
||||
bool syncDft;
|
||||
bool sync;
|
||||
|
||||
uint64_t smpMemPos;
|
||||
uint64_t lastSequence;
|
||||
|
||||
std::complex<double> omega;
|
||||
|
||||
double windowCorretionFactor;
|
||||
struct timespec lastDftCal;
|
||||
double nextDftCalc;
|
||||
double windowCorrectionFactor;
|
||||
struct timespec lastCalc;
|
||||
double nextCalc;
|
||||
|
||||
std::vector<int> signalIndex; /**< A list of signalIndex to do dft on */
|
||||
Phasor lastResult;
|
||||
|
||||
std::string dumperPrefix;
|
||||
bool dumperEnable;
|
||||
#ifdef DFT_MEM_DUMP
|
||||
Dumper origSigSync;
|
||||
Dumper windowdSigSync;
|
||||
Dumper ppsSigSync;
|
||||
#endif
|
||||
Dumper phasorRocof;
|
||||
Dumper phasorPhase;
|
||||
Dumper phasorAmplitude;
|
||||
Dumper phasorFreq;
|
||||
|
||||
public:
|
||||
DftHook(struct vpath *p, struct vnode *n, int fl, int prio, bool en = true) :
|
||||
Hook(p, n, fl, prio, en),
|
||||
|
|
@ -126,48 +130,44 @@ public:
|
|||
paddingType(PaddingType::ZERO),
|
||||
freqEstType(FreqEstimationType::NONE),
|
||||
smpMemory(),
|
||||
#ifdef DFT_MEM_DUMP
|
||||
ppsMemory(),
|
||||
dftMatrix(),
|
||||
dftResults(),
|
||||
#endif
|
||||
matrix(),
|
||||
results(),
|
||||
filterWindowCoefficents(),
|
||||
absDftResults(),
|
||||
absDftFreqs(),
|
||||
dftCalcCount(0),
|
||||
absResults(),
|
||||
absFrequencies(),
|
||||
calcCount(0),
|
||||
sampleRate(0),
|
||||
startFreqency(0),
|
||||
startFrequency(0),
|
||||
endFreqency(0),
|
||||
frequencyResolution(0),
|
||||
dftRate(0),
|
||||
rate(0),
|
||||
ppsIndex(0),
|
||||
windowSize(0),
|
||||
windowMultiplier(0),
|
||||
freqCount(0),
|
||||
syncDft(0),
|
||||
sync(0),
|
||||
smpMemPos(0),
|
||||
lastSequence(0),
|
||||
windowCorretionFactor(0),
|
||||
lastDftCal({0, 0}),
|
||||
nextDftCalc(0.0),
|
||||
windowCorrectionFactor(0),
|
||||
lastCalc({0, 0}),
|
||||
nextCalc(0.0),
|
||||
signalIndex(),
|
||||
lastResult({0,0,0,0})
|
||||
{
|
||||
logger = logging.get("hook:dft");
|
||||
|
||||
format = format_type_lookup("villas.human");
|
||||
|
||||
if (logger->level() <= SPDLOG_LEVEL_DEBUG) {
|
||||
lastResult({0,0,0,0}),
|
||||
dumperPrefix("/tmp/plot/"),
|
||||
dumperEnable(logger->level() <= SPDLOG_LEVEL_DEBUG),
|
||||
#ifdef DFT_MEM_DUMP
|
||||
origSigSync = std::make_shared<Dumper>("/tmp/plot/origSigSync");
|
||||
windowdSigSync = std::make_shared<Dumper>("/tmp/plot/windowdSigSync");
|
||||
ppsSigSync = std::make_shared<Dumper>("/tmp/plot/ppsSigSync");
|
||||
origSigSync(dumperPrefix + "origSigSync"),
|
||||
windowdSigSync(dumperPrefix + "windowdSigSync"),
|
||||
ppsSigSync(dumperPrefix + "ppsSigSync"),
|
||||
#endif
|
||||
origSigSync = std::make_shared<Dumper>("/tmp/plot/origSigSync");
|
||||
phasorPhase = std::make_shared<Dumper>("/tmp/plot/phasorPhase");
|
||||
phasorAmplitude = std::make_shared<Dumper>("/tmp/plot/phasorAmplitude");
|
||||
phasorFreq = std::make_shared<Dumper>("/tmp/plot/phasorFreq");
|
||||
|
||||
}
|
||||
}
|
||||
phasorRocof(dumperPrefix + "phasorRocof"),
|
||||
phasorPhase(dumperPrefix + "phasorPhase"),
|
||||
phasorAmplitude(dumperPrefix + "phasorAmplitude"),
|
||||
phasorFreq(dumperPrefix + "phasorFreq")
|
||||
{ }
|
||||
|
||||
virtual void prepare()
|
||||
{
|
||||
|
|
@ -191,8 +191,6 @@ public:
|
|||
vlist_push(&signals, amplSig);
|
||||
vlist_push(&signals, phaseSig);
|
||||
vlist_push(&signals, rocofSig);
|
||||
|
||||
|
||||
}
|
||||
|
||||
/* Initialize sample memory */
|
||||
|
|
@ -200,31 +198,33 @@ public:
|
|||
for (unsigned i = 0; i < signalIndex.size(); i++)
|
||||
smpMemory.emplace_back(windowSize, 0.0);
|
||||
|
||||
#ifdef DFT_MEM_DUMP
|
||||
/* Initialize temporary ppsMemory */
|
||||
ppsMemory.clear();
|
||||
ppsMemory.resize(windowSize, 0.0);
|
||||
#endif
|
||||
|
||||
/* Calculate how much zero padding ist needed for a needed resolution */
|
||||
windowMultiplier = ceil(((double) sampleRate / windowSize) / frequencyResolution);
|
||||
|
||||
freqCount = ceil((endFreqency - startFreqency) / frequencyResolution) + 1;
|
||||
freqCount = ceil((endFreqency - startFrequency) / frequencyResolution) + 1;
|
||||
|
||||
/* Initialize matrix of dft coeffients */
|
||||
dftMatrix.clear();
|
||||
matrix.clear();
|
||||
for (unsigned i = 0; i < freqCount; i++)
|
||||
dftMatrix.emplace_back(windowSize * windowMultiplier, 0.0);
|
||||
matrix.emplace_back(windowSize * windowMultiplier, 0.0);
|
||||
|
||||
/* Initalize dft results matrix */
|
||||
dftResults.clear();
|
||||
results.clear();
|
||||
for (unsigned i = 0; i < signalIndex.size(); i++){
|
||||
dftResults.emplace_back(freqCount, 0.0);
|
||||
absDftResults.emplace_back(freqCount, 0.0);
|
||||
results.emplace_back(freqCount, 0.0);
|
||||
absResults.emplace_back(freqCount, 0.0);
|
||||
}
|
||||
|
||||
filterWindowCoefficents.resize(windowSize);
|
||||
|
||||
for (unsigned i = 0; i < freqCount; i++) {
|
||||
absDftFreqs.emplace_back(startFreqency + i * frequencyResolution);
|
||||
absFrequencies.emplace_back(startFrequency + i * frequencyResolution);
|
||||
}
|
||||
|
||||
generateDftMatrix();
|
||||
|
|
@ -233,34 +233,38 @@ public:
|
|||
state = State::PREPARED;
|
||||
}
|
||||
|
||||
virtual void parse(json_t *cfg)
|
||||
virtual void parse(json_t *json)
|
||||
{
|
||||
const char *paddingTypeC = nullptr, *windowTypeC = nullptr, *freqEstimateTypeC = nullptr;
|
||||
int ret;
|
||||
json_error_t err;
|
||||
int windowSizeFactor = 1;
|
||||
|
||||
const char *paddingTypeC = nullptr;
|
||||
const char *windowTypeC = nullptr;
|
||||
const char *freqEstimateTypeC = nullptr;
|
||||
|
||||
json_error_t err;
|
||||
json_t *jsonChannelList = nullptr;
|
||||
|
||||
assert(state != State::STARTED);
|
||||
|
||||
Hook::parse(cfg);
|
||||
Hook::parse(json);
|
||||
|
||||
ret = json_unpack_ex(cfg, &err, 0, "{ s?: i, s?: F, s?: F, s?: F, s?: i , s?: i, s?: s, s?: s, s?: s, s?: b, s?: o}",
|
||||
ret = json_unpack_ex(json, &err, 0, "{ s?: i, s?: F, s?: F, s?: F, s?: i , s?: i, s?: s, s?: s, s?: s, s?: b, s?: o }",
|
||||
"sample_rate", &sampleRate,
|
||||
"start_freqency", &startFreqency,
|
||||
"start_freqency", &startFrequency,
|
||||
"end_freqency", &endFreqency,
|
||||
"frequency_resolution", &frequencyResolution,
|
||||
"dft_rate", &dftRate,
|
||||
"window_size", &windowSize,
|
||||
"dft_rate", &rate,
|
||||
"window_size_factor", &windowSizeFactor,
|
||||
"window_type", &windowTypeC,
|
||||
"padding_type", &paddingTypeC,
|
||||
"freq_estimate_type", &freqEstimateTypeC,
|
||||
"sync", &syncDft,
|
||||
"sync", &sync,
|
||||
"signal_index", &jsonChannelList,
|
||||
"pps_index", &ppsIndex
|
||||
);
|
||||
if (ret)
|
||||
throw ConfigError(cfg, err, "node-config-hook-dft");
|
||||
throw ConfigError(json, err, "node-config-hook-dft");
|
||||
|
||||
if (jsonChannelList != nullptr) {
|
||||
signalIndex.clear();
|
||||
|
|
@ -285,42 +289,31 @@ public:
|
|||
signalIndex.push_back(idx);
|
||||
}
|
||||
else
|
||||
logger->warn("Could not parse channel list. Please check documentation for syntax");
|
||||
throw ConfigError(jsonChannelList, "node-config-hook-dft-signal-index", "Could not parse channel list.");
|
||||
}
|
||||
else
|
||||
throw ConfigError(jsonChannelList, "node-config-node-signal", "No parameter signalIndex given.");
|
||||
|
||||
if (!windowSize) {
|
||||
windowSize = (int)(sampleRate * 1 / (double)dftRate);
|
||||
logger->info("Set windows size to {} samples which fits 1/dftRate {}s", windowSize, 1/(double)dftRate);
|
||||
|
||||
}
|
||||
windowSize = sampleRate * windowSizeFactor / (double) rate;
|
||||
logger->debug("Set windows size to {} samples which fits 1 / rate {}s", windowSize, 1.0 / rate);
|
||||
|
||||
if (!windowTypeC) {
|
||||
if (!windowTypeC)
|
||||
logger->info("No Window type given, assume no windowing");
|
||||
windowType = WindowType::NONE;
|
||||
}
|
||||
else if (strcmp(windowTypeC, "flattop") == 0)
|
||||
windowType = WindowType::FLATTOP;
|
||||
else if (strcmp(windowTypeC, "hamming") == 0)
|
||||
windowType = WindowType::HAMMING;
|
||||
else if (strcmp(windowTypeC, "hann") == 0)
|
||||
windowType = WindowType::HANN;
|
||||
else {
|
||||
logger->info("Window type {} not recognized, assume no windowing", windowTypeC);
|
||||
windowType = WindowType::NONE;
|
||||
}
|
||||
else
|
||||
throw ConfigError(json, "node-config-hook-dft-window-type", "Invalid window type: {}", windowTypeC);
|
||||
|
||||
if (!paddingTypeC) {
|
||||
if (!paddingTypeC)
|
||||
logger->info("No Padding type given, assume no zeropadding");
|
||||
paddingType = PaddingType::ZERO;
|
||||
}
|
||||
else if (strcmp(paddingTypeC, "signal_repeat") == 0)
|
||||
paddingType = PaddingType::SIG_REPEAT;
|
||||
else {
|
||||
logger->info("Padding type {} not recognized, assume zero padding", paddingTypeC);
|
||||
paddingType = PaddingType::ZERO;
|
||||
}
|
||||
else
|
||||
throw ConfigError(json, "node-config-hook-dft-padding-type", "Padding type {} not recognized", paddingTypeC);
|
||||
|
||||
if (!freqEstimateTypeC) {
|
||||
logger->info("No Frequency estimation type given, assume no none");
|
||||
|
|
@ -329,15 +322,22 @@ public:
|
|||
else if (strcmp(freqEstimateTypeC, "quadratic") == 0)
|
||||
freqEstType = FreqEstimationType::QUADRATIC;
|
||||
|
||||
if (endFreqency < 0 || endFreqency > sampleRate)
|
||||
throw ConfigError(cfg, err, "node-config-hook-dft", "End frequency must be smaller than sampleRate {}", sampleRate);
|
||||
|
||||
if (frequencyResolution > ((double)sampleRate/windowSize))
|
||||
throw ConfigError(cfg, err, "node-config-hook-dft", "The maximum frequency resolution with smaple_rate:{} and window_site:{} is {}", sampleRate, windowSize, ((double)sampleRate/windowSize));
|
||||
|
||||
state = State::PARSED;
|
||||
}
|
||||
|
||||
virtual void check()
|
||||
{
|
||||
assert(state == State::PARSED);
|
||||
|
||||
if (endFreqency < 0 || endFreqency > sampleRate)
|
||||
throw RuntimeError("End frequency must be smaller than sampleRate {}", sampleRate);
|
||||
|
||||
if (frequencyResolution > (double) sampleRate / windowSize)
|
||||
throw RuntimeError("The maximum frequency resolution with smaple_rate:{} and window_site:{} is {}", sampleRate, windowSize, ((double)sampleRate/windowSize));
|
||||
|
||||
state = State::CHECKED;
|
||||
}
|
||||
|
||||
virtual Hook::Reason process(struct sample *smp)
|
||||
{
|
||||
assert(state == State::STARTED);
|
||||
|
|
@ -345,40 +345,39 @@ public:
|
|||
for (unsigned i = 0; i < signalIndex.size(); i++)
|
||||
smpMemory[i][smpMemPos % windowSize] = smp->data[signalIndex[i]].f;
|
||||
|
||||
// Debugging for pps signal this should only be temporary
|
||||
if (ppsSigSync)
|
||||
ppsMemory[smpMemPos % windowSize] = smp->data[ppsIndex].f;
|
||||
|
||||
#ifdef DFT_MEM_DUMP
|
||||
ppsMemory[smpMemPos % windowSize] = smp->data[ppsIndex].f;
|
||||
#endif
|
||||
smpMemPos++;
|
||||
|
||||
bool runDft = false;
|
||||
if (syncDft) {
|
||||
bool run = false;
|
||||
if (sync) {
|
||||
double smpNsec = smp->ts.origin.tv_sec * 1e9 + smp->ts.origin.tv_nsec;
|
||||
|
||||
if (smpNsec > nextDftCalc) {
|
||||
runDft = true;
|
||||
nextDftCalc = (( smp->ts.origin.tv_sec ) + ( ((dftCalcCount % dftRate) + 1) / (double)dftRate )) * 1e9;
|
||||
if (smpNsec > nextCalc) {
|
||||
run = true;
|
||||
nextCalc = (smp->ts.origin.tv_sec + (((calcCount % rate) + 1) / (double) rate)) * 1e9;
|
||||
}
|
||||
}
|
||||
|
||||
if (runDft) {
|
||||
lastDftCal = smp->ts.origin;
|
||||
if (run) {
|
||||
lastCalc = smp->ts.origin;
|
||||
|
||||
// Debugging for pps signal this should only be temporary
|
||||
if (ppsSigSync) {
|
||||
double tmpPPSWindow[windowSize];
|
||||
#ifdef DFT_MEM_DUMP
|
||||
double tmpPPSWindow[windowSize];
|
||||
|
||||
for (unsigned i = 0; i< windowSize; i++)
|
||||
tmpPPSWindow[i] = ppsMemory[(i + smpMemPos) % windowSize];
|
||||
for (unsigned i = 0; i< windowSize; i++)
|
||||
tmpPPSWindow[i] = ppsMemory[(i + smpMemPos) % windowSize];
|
||||
|
||||
if (dumperEnable)
|
||||
ppsSigSync.writeDataBinary(windowSize, tmpPPSWindow);
|
||||
#endif
|
||||
|
||||
ppsSigSync->writeDataBinary(windowSize, tmpPPSWindow);
|
||||
}
|
||||
|
||||
#pragma omp parallel for
|
||||
for (unsigned i = 0; i < signalIndex.size(); i++) {
|
||||
Phasor currentResult = {0,0,0,0};
|
||||
|
||||
calculateDft(PaddingType::ZERO, smpMemory[i], dftResults[i], smpMemPos);
|
||||
|
||||
calculateDft(PaddingType::ZERO, smpMemory[i], results[i], smpMemPos);
|
||||
|
||||
unsigned maxPos = 0;
|
||||
|
||||
|
|
@ -386,10 +385,10 @@ public:
|
|||
int multiplier = paddingType == PaddingType::ZERO
|
||||
? 1
|
||||
: windowMultiplier;
|
||||
absDftResults[i][j] = abs(dftResults[i][j]) * 2 / (windowSize * windowCorretionFactor * multiplier);
|
||||
if (currentResult.amplitude < absDftResults[i][j]) {
|
||||
currentResult.frequency = absDftFreqs[j];
|
||||
currentResult.amplitude = absDftResults[i][j];
|
||||
absResults[i][j] = abs(results[i][j]) * 2 / (windowSize * windowCorrectionFactor * multiplier);
|
||||
if (currentResult.amplitude < absResults[i][j]) {
|
||||
currentResult.frequency = absFrequencies[j];
|
||||
currentResult.amplitude = absResults[i][j];
|
||||
maxPos = j;
|
||||
}
|
||||
}
|
||||
|
|
@ -398,10 +397,10 @@ public:
|
|||
if (maxPos < 1 || maxPos >= freqCount - 1)
|
||||
logger->warn("Maximum frequency bin lies on window boundary. Using non-estimated results!");
|
||||
else {
|
||||
Point a = { absDftFreqs[maxPos - 1], absDftResults[i][maxPos - 1] };
|
||||
Point b = { absDftFreqs[maxPos + 0], absDftResults[i][maxPos + 0] };
|
||||
Point c = { absDftFreqs[maxPos + 1], absDftResults[i][maxPos + 1] };
|
||||
|
||||
Point a = { absFrequencies[maxPos - 1], absResults[i][maxPos - 1] };
|
||||
Point b = { absFrequencies[maxPos + 0], absResults[i][maxPos + 0] };
|
||||
Point c = { absFrequencies[maxPos + 1], absResults[i][maxPos + 1] };
|
||||
|
||||
Point estimate = quadraticEstimation(a, b, c, maxPos);
|
||||
currentResult.frequency = estimate.x;
|
||||
currentResult.amplitude = estimate.y;
|
||||
|
|
@ -412,32 +411,25 @@ public:
|
|||
|
||||
smp->data[i * 4 + 0].f = currentResult.frequency; /* Frequency */
|
||||
smp->data[i * 4 + 1].f = (currentResult.amplitude / pow(2, 0.5)); /* Amplitude */
|
||||
smp->data[i * 4 + 2].f = atan2(dftResults[i][maxPos].imag(), dftResults[i][maxPos].real()); /* Phase */
|
||||
smp->data[i * 4 + 3].f = (currentResult.frequency - lastResult.frequency) / (double)dftRate; /* RoCof */
|
||||
smp->data[i * 4 + 2].f = atan2(results[i][maxPos].imag(), results[i][maxPos].real()); /* Phase */
|
||||
smp->data[i * 4 + 3].f = (currentResult.frequency - lastResult.frequency) / (double)rate; /* RoCof */
|
||||
|
||||
}
|
||||
|
||||
lastResult = currentResult;
|
||||
}
|
||||
|
||||
//the following is a debug output and currently only for channel 0
|
||||
if (windowSize * 5 < smpMemPos){
|
||||
if (phasorFreq)
|
||||
phasorFreq->writeDataBinary(1, &(smp->data[0 * 4 + 0].f));
|
||||
|
||||
if (phasorPhase)
|
||||
phasorPhase->writeDataBinary(1, &(smp->data[0 * 4 + 2].f));
|
||||
|
||||
if (phasorAmplitude)
|
||||
phasorAmplitude->writeDataBinary(1, &(smp->data[0 * 4 + 1].f));
|
||||
|
||||
if (origSigSync)
|
||||
origSigSync->writeDataBinary(1, &(smp->data[0 * 4 + 3].f));
|
||||
// The following is a debug output and currently only for channel 0
|
||||
if (dumperEnable && windowSize * 5 < smpMemPos){
|
||||
phasorFreq.writeDataBinary(1, &(smp->data[0 * 4 + 0].f));
|
||||
phasorPhase.writeDataBinary(1, &(smp->data[0 * 4 + 2].f));
|
||||
phasorAmplitude.writeDataBinary(1, &(smp->data[0 * 4 + 1].f));
|
||||
phasorRocof.writeDataBinary(1, &(smp->data[0 * 4 + 3].f));
|
||||
}
|
||||
|
||||
smp->length = windowSize < smpMemPos ? signalIndex.size() * 4 : 0;
|
||||
|
||||
dftCalcCount++;
|
||||
calcCount++;
|
||||
}
|
||||
|
||||
if (smp->sequence - lastSequence > 1)
|
||||
|
|
@ -445,7 +437,7 @@ public:
|
|||
|
||||
lastSequence = smp->sequence;
|
||||
|
||||
if(runDft && windowSize < smpMemPos)
|
||||
if (run && windowSize < smpMemPos)
|
||||
return Reason::OK;
|
||||
|
||||
return Reason::SKIP_SAMPLE;
|
||||
|
|
@ -459,17 +451,16 @@ public:
|
|||
using namespace std::complex_literals;
|
||||
|
||||
omega = exp((-2i * M_PI) / (double)(windowSize * windowMultiplier));
|
||||
unsigned startBin = floor(startFreqency / frequencyResolution);
|
||||
unsigned startBin = floor(startFrequency / frequencyResolution);
|
||||
|
||||
for (unsigned i = 0; i < freqCount ; i++) {
|
||||
for (unsigned j = 0 ; j < windowSize * windowMultiplier ; j++)
|
||||
dftMatrix[i][j] = pow(omega, (i + startBin) * j);
|
||||
matrix[i][j] = pow(omega, (i + startBin) * j);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* This function calculates the discrete furie transform of the input signal
|
||||
* This function calculates the discrete furie transform of the input signal
|
||||
*/
|
||||
void calculateDft(enum PaddingType padding, std::vector<double> &ringBuffer, std::vector<std::complex<double>> &results, unsigned ringBufferPos)
|
||||
{
|
||||
|
|
@ -481,20 +472,16 @@ public:
|
|||
tmpSmpWindow[i] = ringBuffer[(i + ringBufferPos) % windowSize];
|
||||
|
||||
#ifdef DFT_MEM_DUMP
|
||||
|
||||
if (origSigSync)
|
||||
origSigSync->writeDataBinary(windowSize, tmpSmpWindow);
|
||||
|
||||
if (dumperEnable)
|
||||
origSigSync.writeDataBinary(windowSize, tmpSmpWindow);
|
||||
#endif
|
||||
|
||||
for (unsigned i = 0; i < windowSize; i++)
|
||||
tmpSmpWindow[i] *= filterWindowCoefficents[i];
|
||||
|
||||
#ifdef DFT_MEM_DUMP
|
||||
|
||||
if (windowdSigSync)
|
||||
windowdSigSync->writeDataBinary(windowSize, tmpSmpWindow);
|
||||
|
||||
if (dumperEnable)
|
||||
windowdSigSync.writeDataBinary(windowSize, tmpSmpWindow);
|
||||
#endif
|
||||
|
||||
for (unsigned i = 0; i < freqCount; i++) {
|
||||
|
|
@ -503,12 +490,12 @@ public:
|
|||
for (unsigned j = 0; j < windowSize * windowMultiplier; j++) {
|
||||
if (padding == PaddingType::ZERO) {
|
||||
if (j < (windowSize))
|
||||
results[i] += tmpSmpWindow[j] * dftMatrix[i][j];
|
||||
results[i] += tmpSmpWindow[j] * matrix[i][j];
|
||||
else
|
||||
results[i] += 0;
|
||||
}
|
||||
else if (padding == PaddingType::SIG_REPEAT) /* Repeat samples */
|
||||
results[i] += tmpSmpWindow[j % windowSize] * dftMatrix[i][j];
|
||||
results[i] += tmpSmpWindow[j % windowSize] * matrix[i][j];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -526,7 +513,7 @@ public:
|
|||
+ 0.277263158 * cos(4 * M_PI * i / (windowSize))
|
||||
- 0.083578947 * cos(6 * M_PI * i / (windowSize))
|
||||
+ 0.006947368 * cos(8 * M_PI * i / (windowSize));
|
||||
windowCorretionFactor += filterWindowCoefficents[i];
|
||||
windowCorrectionFactor += filterWindowCoefficents[i];
|
||||
}
|
||||
break;
|
||||
|
||||
|
|
@ -538,7 +525,7 @@ public:
|
|||
|
||||
for (unsigned i = 0; i < windowSize; i++) {
|
||||
filterWindowCoefficents[i] = a0 - (1 - a0) * cos(2 * M_PI * i / (windowSize));
|
||||
windowCorretionFactor += filterWindowCoefficents[i];
|
||||
windowCorrectionFactor += filterWindowCoefficents[i];
|
||||
}
|
||||
|
||||
break;
|
||||
|
|
@ -547,33 +534,33 @@ public:
|
|||
default:
|
||||
for (unsigned i = 0; i < windowSize; i++) {
|
||||
filterWindowCoefficents[i] = 1;
|
||||
windowCorretionFactor += filterWindowCoefficents[i];
|
||||
windowCorrectionFactor += filterWindowCoefficents[i];
|
||||
}
|
||||
break;
|
||||
}
|
||||
|
||||
windowCorretionFactor /= windowSize;
|
||||
windowCorrectionFactor /= windowSize;
|
||||
}
|
||||
|
||||
/**
|
||||
* This function is calculating the mximum based on a quadratic interpolation
|
||||
*
|
||||
*
|
||||
* This function is based on the following paper:
|
||||
* https://mgasior.web.cern.ch/pap/biw2004.pdf
|
||||
* https://dspguru.com/dsp/howtos/how-to-interpolate-fft-peak/
|
||||
* *
|
||||
* *
|
||||
* In particular equation 10
|
||||
*/
|
||||
Point quadraticEstimation(const Point &a, const Point &b, const Point &c, unsigned maxFBin)
|
||||
{
|
||||
// Frequency estimation
|
||||
double maxBinEst = (double) maxFBin + (c.y - a.y) / (2 * (2 * b.y - a.y - c.y));
|
||||
double y_Fmax = startFreqency + maxBinEst * frequencyResolution; // convert bin to frequency
|
||||
double y_Fmax = startFrequency + maxBinEst * frequencyResolution; // convert bin to frequency
|
||||
|
||||
// Amplitude estimation
|
||||
double f = (a.x * (b.y - c.y) + b.x * (c.y - a.y) + c.x * (a.y - b.y)) / ((a.x - b.x) * (a.x - c.x) * (c.x - b.x));
|
||||
double g = (pow(a.x, 2) * (b.y - c.y) + pow(b.x, 2) * (c.y - a.y) + pow(c.x, 2) * (a.y - b.y)) / ((a.x - b.x) * (a.x - c.x) * (b.x - c.x));
|
||||
double h = (pow(a.x, 2) * (b.x * c.y - c.x * b.y) + a.x * (pow(c.x, 2) * b.y - pow(b.x,2) * c.y)+ b.x * c.x * a.y * (b.x - c.x)) / ((a.x - b.x) * (a.x - c.x) * (b.x - c.x));
|
||||
double h = (pow(a.x, 2) * (b.x * c.y - c.x * b.y) + a.x * (pow(c.x, 2) * b.y - pow(b.x,2) * c.y)+ b.x * c.x * a.y * (b.x - c.x)) / ((a.x - b.x) * (a.x - c.x) * (b.x - c.x));
|
||||
double i = f * pow(y_Fmax,2) + g * y_Fmax + h;
|
||||
|
||||
return { y_Fmax, i };
|
||||
|
|
|
|||
|
|
@ -100,7 +100,8 @@ COPY . /villas/
|
|||
RUN mkdir -p /villas/build
|
||||
WORKDIR /villas/build
|
||||
RUN --security=insecure \
|
||||
cmake -DCMAKE_INSTALL_PREFIX=${PREFIX} \
|
||||
cmake -DWITH_OPENMP=OFF \
|
||||
-DCMAKE_INSTALL_PREFIX=${PREFIX} \
|
||||
-DCMAKE_PREFIX_PATH=${PREFIX} .. && \
|
||||
make -j8 install
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue