diff --git a/CMakeLists.txt b/CMakeLists.txt index adad61ddc..d641d24ac 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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) diff --git a/lib/hooks/dft.cpp b/lib/hooks/dft.cpp index 4c88c534c..009a2bec8 100644 --- a/lib/hooks/dft.cpp +++ b/lib/hooks/dft.cpp @@ -34,9 +34,8 @@ #include #include #include -#include -#include +/* 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 origSigSync; - std::shared_ptr windowdSigSync; - std::shared_ptr phasorPhase; - std::shared_ptr phasorAmplitude; - std::shared_ptr phasorFreq; - std::shared_ptr ppsSigSync; - enum WindowType windowType; enum PaddingType paddingType; enum FreqEstimationType freqEstType; - struct format_type *format; - std::vector> smpMemory; - std::vector ppsMemory;//this is just temporary for debugging - std::vector>> dftMatrix; - std::vector>> dftResults; +#ifdef DFT_MEM_DUMP + std::vector ppsMemory; +#endif + std::vector>> matrix; + std::vector>> results; std::vector filterWindowCoefficents; - std::vector> absDftResults; - std::vector absDftFreqs; + std::vector> absResults; + std::vector 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 omega; - double windowCorretionFactor; - struct timespec lastDftCal; - double nextDftCalc; + double windowCorrectionFactor; + struct timespec lastCalc; + double nextCalc; std::vector 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("/tmp/plot/origSigSync"); - windowdSigSync = std::make_shared("/tmp/plot/windowdSigSync"); - ppsSigSync = std::make_shared("/tmp/plot/ppsSigSync"); + origSigSync(dumperPrefix + "origSigSync"), + windowdSigSync(dumperPrefix + "windowdSigSync"), + ppsSigSync(dumperPrefix + "ppsSigSync"), #endif - origSigSync = std::make_shared("/tmp/plot/origSigSync"); - phasorPhase = std::make_shared("/tmp/plot/phasorPhase"); - phasorAmplitude = std::make_shared("/tmp/plot/phasorAmplitude"); - phasorFreq = std::make_shared("/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 &ringBuffer, std::vector> &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 }; diff --git a/packaging/docker/Dockerfile.alpine b/packaging/docker/Dockerfile.alpine index 38d4c110d..339aa19b6 100644 --- a/packaging/docker/Dockerfile.alpine +++ b/packaging/docker/Dockerfile.alpine @@ -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