From ea446914aed8349ad3126ab4dd38941675041a03 Mon Sep 17 00:00:00 2001 From: Carl Philipp Klemm <philipp@uvos.xyz> Date: Tue, 27 Jun 2023 12:26:46 +0200 Subject: [PATCH] add the ability to calulate the impedacne spectra from drt --- drt.cpp | 34 ++++++++++++++++++++++++++++++---- eisdrt/eigendrt.h | 14 +++++++++++++- eisdrt/eisdrt.h | 20 ++++++++++++++++++++ eisdrt/eistorchdrt.h | 2 ++ eisdrt/torchdrt.h | 9 +++++++++ eistoeigen.h | 17 +++++++++++++++++ main.cpp | 10 ++++++++-- 7 files changed, 99 insertions(+), 7 deletions(-) diff --git a/drt.cpp b/drt.cpp index 46f6c46..456458c 100644 --- a/drt.cpp +++ b/drt.cpp @@ -219,13 +219,25 @@ template Eigen::VectorX<float> calcDrt<float>(Eigen::VectorX<std::complex<float> Eigen::VectorX<float>&, FitMetics& fm, const FitParameters& fp, float* rSeries); template<typename fv> -fv testFn() +Eigen::VectorX<std::complex<fv>> calcImpedance(const Eigen::VectorX<fv>& drt, fv rSeries, Eigen::VectorX<fv>& omegaVector) { - fv value = 0.001; - return value; + Eigen::Matrix<fv, Eigen::Dynamic, Eigen::Dynamic> aMatrixImag = aImag<fv>(omegaVector); + Eigen::Matrix<fv, Eigen::Dynamic, Eigen::Dynamic> aMatrixReal = aReal<fv>(omegaVector); + + Eigen::Matrix<fv, Eigen::Dynamic, Eigen::Dynamic> realMul = aMatrixReal*drt; + Eigen::Matrix<fv, Eigen::Dynamic, Eigen::Dynamic> imagMul = aMatrixImag*drt; + Eigen::Matrix<std::complex<fv>, Eigen::Dynamic, Eigen::Dynamic> realMulCplx = realMul.template cast<std::complex<fv>>(); + Eigen::Matrix<std::complex<fv>, Eigen::Dynamic, Eigen::Dynamic> imagMulCplx = imagMul.template cast<std::complex<fv>>()*std::complex<fv>(0,1); + + Eigen::VectorX<std::complex<fv>> z = ((realMulCplx + imagMulCplx).array() + std::complex<fv>(rSeries,0)).matrix(); + return z; } -template double testFn<double>(); +template Eigen::VectorX<std::complex<double>> calcImpedance<double>(const Eigen::VectorX<double>& drt, double rSeries, + Eigen::VectorX<double>& omegaVector); +template Eigen::VectorX<std::complex<float>> calcImpedance<float>(const Eigen::VectorX<float>& drt, float rSeries, + Eigen::VectorX<float>& omegaVector); + #ifdef USE_EISGEN std::vector<fvalue> calcDrt(const std::vector<eis::DataPoint>& data, const std::vector<fvalue>& omegaVector, @@ -247,5 +259,19 @@ std::vector<fvalue> calcDrt(const std::vector<eis::DataPoint>& data, FitMetics& std::vector<fvalue> stdvector(drt.data(), drt.data()+drt.size()); return stdvector; } + +std::vector<eis::DataPoint> calcImpedance(const std::vector<fvalue>& drt, fvalue rSeries, const std::vector<fvalue>& omegaVector) +{ + Eigen::VectorX<fvalue> omega = Eigen::VectorX<fvalue>::Map(omegaVector.data(), omegaVector.size()); + + Eigen::VectorX<fvalue> drtVec = Eigen::VectorX<fvalue>::Map(drt.data(), drt.size()); + Eigen::VectorX<std::complex<fvalue>> spectra = calcImpedance<fvalue>(drtVec, rSeries, omega); + return eigentoeis(spectra, &omega); +} + +std::vector<eis::DataPoint> calcImpedance(const std::vector<fvalue>& drt, fvalue rSeries, const eis::Range& omegaRange) +{ + return calcImpedance(drt, rSeries, omegaRange.getRangeVector()); +} #endif diff --git a/eisdrt/eigendrt.h b/eisdrt/eigendrt.h index 7cc3338..7b4a50a 100644 --- a/eisdrt/eigendrt.h +++ b/eisdrt/eigendrt.h @@ -40,11 +40,23 @@ Api for use with Eigen applications * @param rSeries an optional paramter where the seires resistance is stored * @return a vector with the drt values */ - template<typename fv> Eigen::VectorX<fv> calcDrt(Eigen::VectorX<std::complex<fv>>& impedanceSpectra, Eigen::VectorX<fv>& omegaTensor, FitMetics& fm, const FitParameters& fp, fv* rSeries = nullptr); + +/** + * @brief calculate impedance from drt using eigen datatypes + * + * @tparam fv precision to be used, either double or float + * @param drt the drt to caluclate impedance from + * @param omegaVector vector with the omega values that the impedances where mesured at + * @param rSeries an optional paramter where the seires resistance is stored + * @return a vector with the drt values + */ +template<typename fv> +Eigen::VectorX<std::complex<fv>> calcImpedance(const Eigen::VectorX<fv>& drt, fv rSeries, const Eigen::VectorX<fv>& omegaVector); + /** .... * @} diff --git a/eisdrt/eisdrt.h b/eisdrt/eisdrt.h index 9052df6..b45db5b 100644 --- a/eisdrt/eisdrt.h +++ b/eisdrt/eisdrt.h @@ -53,6 +53,26 @@ std::vector<fvalue> calcDrt(const std::vector<eis::DataPoint>& data, const std:: */ std::vector<fvalue> calcDrt(const std::vector<eis::DataPoint>& data, FitMetics& fm, const FitParameters& fp, fvalue* rSeries = nullptr); +/** + * @brief calculate impedance from drt using eisgenerator datatypes + * + * @param drt the drt to caluclate impedance from + * @param omegaVector vector with the omega values that the impedances where mesured at + * @param rSeries an optional paramter where the seires resistance is stored + * @return a vector with the impedance values + */ +std::vector<eis::DataPoint> calcImpedance(const std::vector<fvalue>& drt, fvalue rSeries, const std::vector<fvalue>& omegaVector); + +/** + * @brief calculate impedance from drt using eisgenerator datatypes + * + * @param drt the drt to caluclate impedance from + * @param omegaRange range that describes the omega values the drt was taken at + * @param rSeries an optional paramter where the seires resistance is stored + * @return a vector with the impedance values + */ +std::vector<eis::DataPoint> calcImpedance(const std::vector<fvalue>& drt, fvalue rSeries, const eis::Range& omegaRange); + /** .... * @} diff --git a/eisdrt/eistorchdrt.h b/eisdrt/eistorchdrt.h index 4f5b9f3..eb8e9b5 100644 --- a/eisdrt/eistorchdrt.h +++ b/eisdrt/eistorchdrt.h @@ -32,6 +32,8 @@ torch::Tensor calcDrtTorch(const std::vector<eis::DataPoint>& data, const std::v torch::Tensor calcDrtTorch(const std::vector<eis::DataPoint>& data, FitMetics& fm, const FitParameters& fp); +torch::Tensor calcImpedance(torch::Tensor& drt, fvalue rSeries, const std::vector<fvalue>& omegaVector); + /** .... * @} diff --git a/eisdrt/torchdrt.h b/eisdrt/torchdrt.h index 6cbbd3f..3888f38 100644 --- a/eisdrt/torchdrt.h +++ b/eisdrt/torchdrt.h @@ -41,6 +41,15 @@ Api for use with libtorch applications template<typename fv> torch::Tensor calcDrtTorch(torch::Tensor& impedanceSpectra, torch::Tensor& omegaTensor, FitMetics& fm, const FitParameters& fp); +/** + * @brief calculates impedance from drt + * + * @tparam fv precision to be used, either double or float + * @param drt the drt to calcualte impedance from + * @return a complex tensor with the impedance spectra + */ +template<typename fv> +torch::Tensor calcImpedance(torch::Tensor& drt, fv rSeries, torch::Tensor& omegaVector); /** .... diff --git a/eistoeigen.h b/eistoeigen.h index df35147..1aca508 100644 --- a/eistoeigen.h +++ b/eistoeigen.h @@ -21,6 +21,7 @@ #pragma once +#include "Eigen/src/Core/Matrix.h" #include <Eigen/Core> #include <eisgenerator/eistype.h> #include <vector> @@ -40,3 +41,19 @@ Eigen::VectorX<std::complex<fvalue>> eistoeigen(const std::vector<eis::DataPoint } return out; } + +std::vector<eis::DataPoint> eigentoeis(const Eigen::VectorX<std::complex<fvalue>>& data, const Eigen::Vector<fvalue, Eigen::Dynamic>* omega = nullptr) +{ + assert(!omega || omega->size() == data.size()); + std::vector<eis::DataPoint> out(data.size()); + + for(ssize_t i = 0 ; i < data.size(); ++i) + { + out[i].im = data[i]; + if(omega) + out[i].omega = (*omega)[i]; + else + out[i].omega = -1; + } + return out; +} diff --git a/main.cpp b/main.cpp index 63cecf8..3c0469c 100644 --- a/main.cpp +++ b/main.cpp @@ -64,7 +64,7 @@ int main(int argc, char** argv) // specify the angular frequency range of the spectra to be // simulated to be 1-1*10^6 Hz with 3 steps and log10 distrobution - eis::Range omega(1, 1e6, 3, true); + eis::Range omega(1, 1e6, 10, true); // specify circut to be simulated eis::Model model("r{10}-r{50}p{0.02, 0.8}"); @@ -84,11 +84,17 @@ int main(int argc, char** argv) assert(x.size() == data.size()); + std::vector<eis::DataPoint> recalculatedSpectra = calcImpedance(x, rSeries, omega); + fvalue dist = eisNyquistDistance(data, recalculatedSpectra); + // print some info on the drt std::cout<<"Iterations: "<<fm.iterations<<'\n'; std::cout<<"fx "<<fm.fx<<"\ndrt: "; printFvalueVector(x); std::cout<<"r series: "<<rSeries<<'\n'; + std::cout<<"dist: "<<dist<<'\n'; + std::cout<<"recalculatedSpectra:\n"; + printImpedance(recalculatedSpectra); - return 0; + return dist < 2 ? 0 : 1; } -- GitLab