diff --git a/drt.cpp b/drt.cpp
index 46f6c465c005d5420a5b2f71218528a828405775..456458c4fee2fd29cd96bd044c58ffd10e0c52c2 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 7cc3338f526f2711089b3aeb87b6952cd8f24fbc..7b4a50a11fc502fdb3e36c42974b3a7ae267f656 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 9052df665ac61f4e7d6ebdff3e2331e424eadaed..b45db5bff1c23f457479315f1a8b084c2fd5cc21 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 4f5b9f39e4850ce590bdbe4481a75afca232f7ca..eb8e9b552f71140e6cd37c320302cff26dcc03f6 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 6cbbd3f7456977958705e80dbdba4aee13b3339d..3888f38dfffcf84272eb561351476602660986cf 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 df351474fddbbc198c0c2d1407e7efe4bb624077..1aca508385dcda716f0b64e6ca96c403a545d222 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 63cecf8996711670ffccf4620920ceda78c6e772..3c0469ceae84dfa1ee53df3ca4ea62f48445a1ed 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;
}