Fix calibration measurement interpolation bug

Software/PC_Application/LibreVNA-GUI/Calibration/calibrationmeasurement.cpp

@@ -1,6 +1,7 @@
 #include "calibrationmeasurement.h"
 #include "unit.h"
 #include "calibration.h"
+#include "Util/util.h"
 
 #include <QDateTime>
 #include <QComboBox>
@@ -336,14 +337,10 @@ std::complex<double> CalibrationMeasurement::OnePort::getMeasured(double frequen
         return numeric_limits<complex<double>>::quiet_NaN();
     }
     // frequency within points, interpolate
-    auto lower = lower_bound(points.begin(), points.end(), frequency, [](const Point &lhs, double rhs) -> bool {
-        return lhs.frequency < rhs;
+    Point interp = Util::interpolate<Point, double>(points, frequency, [](const Point &p)->double{
+        return p.frequency;
     });
-    auto highPoint = *lower;
-    auto lowPoint = *prev(lower);
-    double alpha = (frequency - lowPoint.frequency) / (highPoint.frequency - lowPoint.frequency);
-    complex<double> ret;
-    return lowPoint.S * (1.0 - alpha) + highPoint.S * alpha;
+    return interp.S;
 }
 
 std::complex<double> CalibrationMeasurement::OnePort::getActual(double frequency)
@@ -522,19 +519,10 @@ Sparam CalibrationMeasurement::TwoPort::getMeasured(double frequency)
         return Sparam();
     }
     // frequency within points, interpolate
-    auto lower = lower_bound(points.begin(), points.end(), frequency, [](const Point &lhs, double rhs) -> bool {
-        return lhs.frequency < rhs;
+    Point interp = Util::interpolate<Point, double>(points, frequency, [](const Point &p)->double{
+        return p.frequency;
     });
-    auto lowPoint = *lower;
-    advance(lower, 1);
-    auto highPoint = *lower;
-    double alpha = (frequency - lowPoint.frequency) / (highPoint.frequency - lowPoint.frequency);
-    Sparam ret;
-    ret.m11 = lowPoint.S.m11 * (1.0 - alpha) + highPoint.S.m11 * alpha;
-    ret.m12 = lowPoint.S.m12 * (1.0 - alpha) + highPoint.S.m12 * alpha;
-    ret.m21 = lowPoint.S.m21 * (1.0 - alpha) + highPoint.S.m21 * alpha;
-    ret.m22 = lowPoint.S.m22 * (1.0 - alpha) + highPoint.S.m22 * alpha;
-    return ret;
+    return interp.S;
 }
 
 Sparam CalibrationMeasurement::TwoPort::getActual(double frequency)
@@ -699,15 +687,14 @@ std::complex<double> CalibrationMeasurement::Isolation::getMeasured(double frequ
     }
     portRcv--;
     portSrc--;
-    // find correct point, no interpolation yet
-    auto lower = lower_bound(points.begin(), points.end(), frequency, [](const Point &lhs, double rhs) -> bool {
-        return lhs.frequency < rhs;
+    // find correct point, interpolate
+    Point interp = Util::interpolate<Point, double>(points, frequency, [](const Point &p)->double{
+        return p.frequency;
     });
-    Point p = *lower;
-    if(portRcv >= p.S.size() || portSrc >= p.S[portRcv].size()) {
+    if(portRcv >= interp.S.size() || portSrc >= interp.S[portRcv].size()) {
         return numeric_limits<complex<double>>::quiet_NaN();
     } else {
-        return p.S[portRcv][portSrc];
+        return interp.S[portRcv][portSrc];
     }
 }
 

Software/PC_Application/LibreVNA-GUI/Calibration/calibrationmeasurement.h

@@ -94,6 +94,19 @@ public:
     public:
         double frequency;
         std::complex<double> S;
+
+        Point operator*(double scalar) const {
+            Point ret;
+            ret.frequency = frequency * scalar;
+            ret.S = S * scalar;
+            return ret;
+        }
+        Point operator+(const Point &p) const {
+            Point ret;
+            ret.frequency = frequency + p.frequency;
+            ret.S = S + p.S;
+            return ret;
+        }
     };
     std::vector<Point> getPoints() const;
 
@@ -195,6 +208,19 @@ public:
     public:
         double frequency;
         Sparam S;
+
+        Point operator*(double scalar) const {
+            Point ret;
+            ret.frequency = frequency * scalar;
+            ret.S = S * scalar;
+            return ret;
+        }
+        Point operator+(const Point &p) const {
+            Point ret;
+            ret.frequency = frequency + p.frequency;
+            ret.S = S + p.S;
+            return ret;
+        }
     };
     std::vector<Point> getPoints() const;
 
@@ -254,6 +280,37 @@ public:
     public:
         double frequency;
         std::vector<std::vector<std::complex<double>>> S;
+
+        Point operator*(double scalar) const {
+            Point ret;
+            ret.frequency = frequency * scalar;
+            for(const auto &v1 : S) {
+                std::vector<std::complex<double>> v;
+                for(const auto &s : v1) {
+                    v.push_back(s * scalar);
+                }
+                ret.S.push_back(v);
+            }
+            return ret;
+        }
+        Point operator+(const Point &p) const {
+            Point ret;
+            ret.frequency = frequency + p.frequency;
+            if(S.size() != p.S.size()) {
+                throw std::runtime_error("Points to not have the same number of measurements");
+            }
+            for(unsigned int i=0;i<S.size();i++) {
+                std::vector<std::complex<double>> v;
+                if(S[i].size() != p.S[i].size()) {
+                    throw std::runtime_error("Points to not have the same number of measurements");
+                }
+                for(unsigned int j=0;j<S[i].size();j++) {
+                    v.push_back(S[i][j] + p.S[i][j]);
+                }
+                ret.S.push_back(v);
+            }
+            return ret;
+        }
     };
     std::vector<Point> getPoints() const;
 

Software/PC_Application/LibreVNA-GUI/Tools/parameters.h

@@ -30,7 +30,7 @@ public:
     Sparam(const Tparam &t);
     Sparam(const ABCDparam &a, Type Z01, Type Z02);
     Sparam(const ABCDparam &a, Type Z0);
-    Sparam operator+(const Sparam &r) {
+    Sparam operator+(const Sparam &r) const {
         Sparam p;
         p.m11 = this->m11+r.m11;
         p.m12 = this->m12+r.m12;
@@ -38,7 +38,7 @@ public:
         p.m22 = this->m22+r.m22;
         return p;
     }
-    Sparam operator*(const Type &r) {
+    Sparam operator*(const Type &r) const {
         Sparam p(m11*r, m12*r, m21*r, m22*r);
         return p;
     }

Software/PC_Application/LibreVNA-GUI/Traces/tracemodel.cpp

@@ -11,6 +11,7 @@ TraceModel::TraceModel(QObject *parent)
 {
     traces.clear();
     source = DataSource::Unknown;
+    lastSweepPosition = 0.0;
 }
 
 TraceModel::~TraceModel()

Software/PC_Application/LibreVNA-GUI/Util/util.h

@@ -69,6 +69,42 @@ namespace Util {
         return brightness > 0.6 ? Qt::black : Qt::white;
     }
 
+    /*
+     * Performs interpolation of a list of sorted values.
+     * T: type of the elements in the list. Must contain a value by which these elements are sorted in the list.
+     * C: Target value of the interpolation (e.g. frequency).
+     * extract: A function returning the "sorted by"-value of the element
+     *
+     * The following operators must be defined:
+     * T*double
+     * T+T
+     *
+     * Will return an interpolated element T. If no interpolation is possible because target is outside
+     * of the list, either the first or last element in the list is returned instead.
+     */
+    template<typename T, typename C>
+    T interpolate(const std::vector<T> &list, C target, std::function<C(const T &lhs)> extract) {
+        T ret = {};
+        if(list.size() > 0) {
+            auto it = std::lower_bound(list.begin(), list.end(), target, [=](const T &lhs, C rhs)->bool{
+                return extract(lhs) < rhs;
+            });
+            if(it == list.begin()) {
+                // just return the first element
+                ret = list.front();
+            } else if(it == list.end()) {
+                // outside of the list of provided values, just use the last one
+                ret = list.back();
+            } else {
+                // needs to interpolate
+                auto highPoint = *it;
+                auto lowPoint = *std::prev(it);
+                double alpha = (target - extract(lowPoint)) / (extract(highPoint) - extract(lowPoint));
+                ret = lowPoint * (1.0 - alpha) + highPoint * alpha;
+            }
+        }
+        return ret;
+    }
     void unwrapPhase(std::vector<double> &phase, unsigned int start_index = 0);
 
     // input values are Y coordinates, assumes evenly spaced linear X values from 0 to input.size() - 1