group delay option for Y axis

Software/PC_Application/LibreVNA-GUI.pro

@@ -218,6 +218,7 @@ SOURCES += \
     Traces/tracewidget.cpp \
     Traces/tracexyplot.cpp \
     Traces/xyplotaxisdialog.cpp \
+    Util/util.cpp \
     VNA/Deembedding/deembedding.cpp \
     VNA/Deembedding/deembeddingdialog.cpp \
     VNA/Deembedding/deembeddingoption.cpp \

Software/PC_Application/Traces/tracexyplot.cpp

@@ -17,22 +17,6 @@
 
 using namespace std;
 
-const set<TraceXYPlot::YAxisType> TraceXYPlot::YAxisTypes = {TraceXYPlot::YAxisType::Disabled,
-                                           TraceXYPlot::YAxisType::Magnitude,
-                                           TraceXYPlot::YAxisType::Phase,
-                                           TraceXYPlot::YAxisType::VSWR,
-                                           TraceXYPlot::YAxisType::Real,
-                                           TraceXYPlot::YAxisType::Imaginary,
-                                           TraceXYPlot::YAxisType::SeriesR,
-                                           TraceXYPlot::YAxisType::Reactance,
-                                           TraceXYPlot::YAxisType::Capacitance,
-                                           TraceXYPlot::YAxisType::Inductance,
-                                           TraceXYPlot::YAxisType::QualityFactor,
-                                           TraceXYPlot::YAxisType::ImpulseReal,
-                                           TraceXYPlot::YAxisType::ImpulseMag,
-                                           TraceXYPlot::YAxisType::Step,
-                                           TraceXYPlot::YAxisType::Impedance};
-
 TraceXYPlot::TraceXYPlot(TraceModel &model, QWidget *parent)
     : TracePlot(model, parent)
 {
@@ -333,8 +317,11 @@ void TraceXYPlot::updateContextMenu()
 
 bool TraceXYPlot::dropSupported(Trace *t)
 {
-    Q_UNUSED(t)
+    if(domainMatch(t) && !supported(t)) {
-    // all kind of traces can be dropped, the graph will be reconfigured to support the dropped trace if required
+        // correct domain configured but Y axis do not match, prevent drop
+        return false;
+    }
+    // either directly compatible or domain change required
     return true;
 }
 
@@ -823,12 +810,14 @@ QString TraceXYPlot::AxisTypeToName(TraceXYPlot::YAxisType type)
     case YAxisType::Capacitance: return "Capacitance";
     case YAxisType::Inductance: return "Inductance";
     case YAxisType::QualityFactor: return "Quality Factor";
+    case YAxisType::GroupDelay: return "Group delay";
     case YAxisType::ImpulseReal: return "Impulse Response (Real)";
     case YAxisType::ImpulseMag: return "Impulse Response (Magnitude)";
     case YAxisType::Step: return "Step Response";
     case YAxisType::Impedance: return "Impedance";
-    default: return "Unknown";
+    case YAxisType::Last: return "Unknown";
     }
+    return "Missing case";
 }
 
 void TraceXYPlot::enableTraceAxis(Trace *t, int axis, bool enabled)
@@ -861,27 +850,26 @@ void TraceXYPlot::enableTraceAxis(Trace *t, int axis, bool enabled)
     }
 }
 
-bool TraceXYPlot::supported(Trace *t, TraceXYPlot::YAxisType type)
+bool TraceXYPlot::domainMatch(Trace *t)
 {
     switch(XAxis.type) {
     case XAxisType::Frequency:
-        if(t->outputType() != Trace::DataType::Frequency) {
+        return t->outputType() == Trace::DataType::Frequency;
-            return false;
-        }
-        break;
     case XAxisType::Distance:
     case XAxisType::Time:
-        if(t->outputType() != Trace::DataType::Time) {
+        return t->outputType() == Trace::DataType::Time;
-            return false;
-        }
-        break;
     case XAxisType::Power:
-        if(t->outputType() != Trace::DataType::Power) {
+        return t->outputType() == Trace::DataType::Power;
+    case XAxisType::Last:
         return false;
     }
-        break;
+    return false;
-    default:
+}
-        break;
+
+bool TraceXYPlot::supported(Trace *t, TraceXYPlot::YAxisType type)
+{
+    if(!domainMatch(t)) {
+        return false;
     }
 
     switch(type) {
@@ -897,6 +885,11 @@ bool TraceXYPlot::supported(Trace *t, TraceXYPlot::YAxisType type)
             return false;
         }
         break;
+    case YAxisType::GroupDelay:
+        if(t->isReflection()) {
+            return false;
+        }
+        break;
     default:
         break;
     }
@@ -946,6 +939,38 @@ QPointF TraceXYPlot::traceToCoordinate(Trace *t, unsigned int sample, TraceXYPlo
     case YAxisType::QualityFactor:
         ret.setY(Util::SparamToQualityFactor(data.y));
         break;
+    case YAxisType::GroupDelay: {
+        constexpr int requiredSamples = 5;
+        if(t->size() < requiredSamples) {
+            // unable to calculate
+            ret.setY(0.0);
+            break;
+        }
+        // needs at least some samples before/after current sample for calculating the derivative.
+        // For samples too far at either end of the trace, return group delay of "inner" trace sample instead
+        if(sample < requiredSamples / 2) {
+            return traceToCoordinate(t, requiredSamples / 2, type);
+        } else if(sample >= t->size() - requiredSamples / 2) {
+            return traceToCoordinate(t, t->size() - requiredSamples / 2 - 1, type);
+        } else {
+            // got enough samples at either end to calculate derivative.
+            // acquire phases of the required samples
+            std::vector<double> phases;
+            phases.reserve(requiredSamples);
+            for(unsigned int index = sample - requiredSamples / 2;index <= sample + requiredSamples / 2;index++) {
+                phases.push_back(arg(t->sample(index).y));
+            }
+            // make sure there are no phase jumps
+            Util::unwrapPhase(phases);
+            // calculate linearRegression to get derivative
+            double B_0, B_1;
+            Util::linearRegression(phases, B_0, B_1);
+            // B_1 now contains the derived phase vs. the sample. Scale by frequency to get group delay
+            double freq_step = t->sample(sample).x - t->sample(sample - 1).x;
+            ret.setY(-B_1 / (2.0*M_PI * freq_step));
+        }
+    }
+        break;
     case YAxisType::ImpulseReal:
         ret.setY(real(data.y));
         break;
@@ -1103,6 +1128,7 @@ QString TraceXYPlot::AxisUnit(TraceXYPlot::YAxisType type)
     case TraceXYPlot::YAxisType::ImpulseMag: return "db";
     case TraceXYPlot::YAxisType::Step: return "";
     case TraceXYPlot::YAxisType::Impedance: return "Ohm";
+    case TraceXYPlot::YAxisType::GroupDelay: return "s";
     default: return "";
     }
 }

Software/PC_Application/Traces/tracexyplot.h

@@ -26,6 +26,7 @@ public:
         Capacitance,
         Inductance,
         QualityFactor,
+        GroupDelay,
         // TDR options
         ImpulseReal,
         ImpulseMag,
@@ -33,7 +34,7 @@ public:
         Impedance,
         Last,
     };
-    static const std::set<YAxisType> YAxisTypes;
+
     enum class XAxisType {
         Frequency,
         Time,
@@ -80,6 +81,7 @@ private:
     YAxisType YAxisTypeFromName(QString name);
     XAxisMode AxisModeFromName(QString name);
     void enableTraceAxis(Trace *t, int axis, bool enabled);
+    bool domainMatch(Trace *t);
     bool supported(Trace *t) override;
     bool supported(Trace *t, YAxisType type);
     QPointF traceToCoordinate(Trace *t, unsigned int sample, YAxisType type);

Software/PC_Application/Traces/xyplotaxisdialog.cpp

@@ -151,7 +151,8 @@ void XYplotAxisDialog::XAxisTypeChanged(int XAxisIndex)
 {
     auto type = (TraceXYPlot::XAxisType) XAxisIndex;
     auto supported = supportedYAxis(type);
-    for(auto t : TraceXYPlot::YAxisTypes) {
+    for(unsigned int i=0;i<(int) TraceXYPlot::YAxisType::Last;i++) {
+        auto t = (TraceXYPlot::YAxisType) i;
         auto enable = supported.count(t) > 0;
         auto index = (int) t;
         enableComboBoxItem(ui->Y1type, index, enable);
@@ -197,6 +198,7 @@ std::set<TraceXYPlot::YAxisType> XYplotAxisDialog::supportedYAxis(TraceXYPlot::X
         ret.insert(TraceXYPlot::YAxisType::Capacitance);
         ret.insert(TraceXYPlot::YAxisType::Inductance);
         ret.insert(TraceXYPlot::YAxisType::QualityFactor);
+        ret.insert(TraceXYPlot::YAxisType::GroupDelay);
         break;
     case TraceXYPlot::XAxisType::Time:
     case TraceXYPlot::XAxisType::Distance:

Software/PC_Application/Util/util.cpp

@@ -0,0 +1,26 @@
+#include "util.h"
+
+void Util::unwrapPhase(std::vector<double> &phase)
+{
+    for (unsigned int i = 1; i < phase.size(); i++) {
+        double d = phase[i] - phase[i-1];
+        d = d > M_PI ? d - 2 * M_PI : (d < -M_PI ? d + 2 * M_PI : d);
+        phase[i] = phase[i-1] + d;
+    }
+}
+
+void Util::linearRegression(const std::vector<double> &input, double &B_0, double &B_1)
+{
+    double x_mean = (input.size() - 1.0) / 2.0;
+    double y_mean = std::accumulate(input.begin(), input.end(), 0.0) / input.size();
+    double ss_xy = 0.0;
+    for(unsigned int i=0;i<input.size();i++) {
+        ss_xy += input[i] * i;
+    }
+    ss_xy -= input.size() * x_mean * y_mean;
+    int n = input.size() - 1;
+    double ss_xx = (1.0/6.0) * n * (n + 1) * (2*n + 1) - input.size() * x_mean * x_mean;
+
+    B_1 = ss_xy / ss_xx;
+    B_0 = y_mean - B_1 * x_mean;
+}

Software/PC_Application/Util/util.h

@@ -4,6 +4,7 @@
 #include <complex>
 #include <math.h>
 #include <limits>
+#include <vector>
 
 #include <QColor>
 
@@ -55,6 +56,11 @@ namespace Util {
         auto brightness = q.redF() * 0.299 + q.greenF() * 0.587 + q.blueF() * 0.114;
         return brightness > 0.6 ? Qt::black : Qt::white;
     }
+
+    void unwrapPhase(std::vector<double> &phase);
+
+    // input values are Y coordinates, assumes evenly spaced linear X values from 0 to input.size() - 1
+    void linearRegression(const std::vector<double> &input, double &B_0, double &B_1);
 }
 
 #endif // UTILH_H