VILLASnode/lib/hooks/pmu_ipdft.cpp

/* ipDFT PMU hook.
 *
 * Author: Manuel Pitz <manuel.pitz@eonerc.rwth-aachen.de>
 * SPDX-FileCopyrightText: 2014-2023 Institute for Automation of Complex Power Systems, RWTH Aachen University
 * SPDX-License-Identifier: Apache-2.0
 */

#include <villas/hooks/pmu.hpp>

namespace villas {
namespace node {

class IpDftPmuHook : public PmuHook {

protected:
	std::complex<double> omega;
	std::vector<std::vector<std::complex<double>>> dftMatrix;
	std::vector<std::complex<double>> dftResult;

	unsigned frequencyCount; // Number of requency bins that are calculated
	double estimationRange; // The range around nominalFreq used for estimation

public:
	IpDftPmuHook(Path *p, Node *n, int fl, int prio, bool en = true):
		PmuHook(p, n, fl, prio, en),
		frequencyCount(0),
		estimationRange(0)

	{ }

	void prepare()
	{
		PmuHook::prepare();

		const double startFrequency = nominalFreq - estimationRange;
		const double endFrequency = nominalFreq + estimationRange;
		const double frequencyResolution = (double)sampleRate / windowSize;

		frequencyCount = ceil((endFrequency - startFrequency) / frequencyResolution);

		// Initialize matrix of dft coeffients
		dftMatrix.clear();
		for (unsigned i = 0; i < frequencyCount; i++)
			dftMatrix.emplace_back(windowSize, 0.0);

		using namespace std::complex_literals;

		omega = exp((-2i * M_PI) / (double)(windowSize));
		unsigned startBin = floor(startFrequency / frequencyResolution);

		for (unsigned i = 0; i <  frequencyCount ; i++) {
			for (unsigned j = 0 ; j < windowSize; j++)
				dftMatrix[i][j] = pow(omega, (i + startBin) * j);
			dftResult.push_back(0.0);
		}
	}

	void parse(json_t *json)
	{
		PmuHook::parse(json);
		int ret;

		json_error_t err;

		assert(state != State::STARTED);

		Hook::parse(json);

		ret = json_unpack_ex(json, &err, 0, "{ s?: F}",
			"estimation_range", &estimationRange
		);

		if (ret)
			throw ConfigError(json, err, "node-config-hook-ip-dft-pmu");

		if (estimationRange <= 0)
			throw ConfigError(json, "node-config-hook-ip-dft-pmu-estimation_range", "Estimation range cannot be less or equal than 0 tried to set {}", estimationRange);

	}

	PmuHook::Phasor estimatePhasor(dsp::CosineWindow<double> *window, PmuHook::Phasor lastPhasor)
	{
		PmuHook::Phasor phasor = {0};

		// Calculate DFT
		for (unsigned i = 0; i < frequencyCount; i++) {
			dftResult[i] = 0;

			const unsigned size = (*window).size();
			for (unsigned j = 0; j < size; j++)
				dftResult[i] += (*window)[j] * dftMatrix[i][j];
		}
		// End calculate DFT

		// Find max bin
		unsigned maxBin = 0;
		double absAmplitude = 0;

		for(unsigned j = 0; j < frequencyCount; j++) {
			if (absAmplitude < abs(dftResult[j])) {
				absAmplitude = abs(dftResult[j]);
				maxBin = j;
			}
		}
		// End find max bin

		if (maxBin == 0 || maxBin == (frequencyCount - 1)) {
			logger->warn("Maximum frequency bin lies on window boundary. Using non-estimated results!");
			//TODO: add handling to not forward this phasor!!
		} else {
			const double startFrequency = nominalFreq - estimationRange;
			const double frequencyResolution = (double)sampleRate / windowSize;

			double a = abs(dftResult[ maxBin - 1 ]);
			double b = abs(dftResult[ maxBin + 0 ]);
			double c = abs(dftResult[ maxBin + 1 ]);
			double bPhase = atan2(dftResult[maxBin].imag(), dftResult[maxBin].real());

			// Estimate phasor
			// Based on https://ieeexplore.ieee.org/document/7980868
			double delta = 0;

			// Paper eq 8
			if (c > a) {
				delta = 1. * (2. * c - b) / (b + c);
			} else {
				delta = -1. * (2. * a - b) / (b + a);
			}

			// Frequency estimation (eq 4)
			phasor.frequency =  startFrequency + ( (double) maxBin + delta) * frequencyResolution;

			// Amplitude estimation (eq 9)
			phasor.amplitude = b * abs( (M_PI * delta) / sin( M_PI * delta) ) * abs( pow(delta, 2) - 1);
			phasor.amplitude *=  2 / (windowSize * window->getCorrectionFactor());

			// Phase estimation (eq 10)
			phasor.phase = bPhase - M_PI * delta;

			// ROCOF estimation
			phasor.rocof = ((phasor.frequency - lastPhasor.frequency) * (double)phasorRate);
			// End estimate phasor
		}

		if (lastPhasor.frequency != 0) // Check if we already calculated a phasor before
			phasor.valid = Status::VALID;

		return phasor;
	}
};

// Register hook
static char n[] = "ip-dft-pmu";
static char d[] = "This hook calculates a phasor based on ipDFT";
static HookPlugin<IpDftPmuHook, n, d, (int) Hook::Flags::NODE_READ | (int) Hook::Flags::NODE_WRITE | (int) Hook::Flags::PATH> p;

} // namespace node
} // namespace villas