caiyu/LibreVNA
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Sliding load support + TRL standard and measurements

Browse Source
This commit is contained in:
Jan Käberich 2022-10-01 16:12:33 +02:00
parent 858e945373
commit fc1d006edd
11 changed files with 621 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

147
Software/PC_Application/Calibration/CalStandardLineEditDialog.ui Normal file
View File

@ -0,0 +1,147 @@
<?xml version="1.0" encoding="UTF-8"?>
<ui version="4.0">
<class>CalStandardLineEditDialog</class>
<widget class="QDialog" name="CalStandardLineEditDialog">
<property name="geometry">
<rect>
<x>0</x>
<y>0</y>
<width>383</width>
<height>202</height>
</rect>
</property>
<property name="windowTitle">
<string>Edit Line Standard</string>
</property>
<property name="modal">
<bool>true</bool>
</property>
<layout class="QVBoxLayout" name="verticalLayout">
<item>
<layout class="QFormLayout" name="formLayout_2">
<item row="0" column="0">
<widget class="QLabel" name="label">
<property name="text">
<string>Name:</string>
</property>
</widget>
</item>
<item row="0" column="1">
<widget class="QLineEdit" name="name"/>
</item>
</layout>
</item>
<item>
<layout class="QFormLayout" name="formLayout">
<item row="0" column="0">
<widget class="QLabel" name="label_35">
<property name="text">
<string>Z0:</string>
</property>
</widget>
</item>
<item row="0" column="1">
<widget class="SIUnitEdit" name="Z0">
<property name="enabled">
<bool>false</bool>
</property>
</widget>
</item>
<item row="1" column="0">
<widget class="QLabel" name="label_21">
<property name="text">
<string>Delay:</string>
</property>
</widget>
</item>
<item row="1" column="1">
<widget class="SIUnitEdit" name="delay"/>
</item>
<item row="2" column="0">
<widget class="QLabel" name="label_2">
<property name="text">
<string>Minimum usable frequency:</string>
</property>
</widget>
</item>
<item row="2" column="1">
<widget class="SIUnitEdit" name="minFreq">
<property name="enabled">
<bool>false</bool>
</property>
</widget>
</item>
<item row="3" column="0">
<widget class="QLabel" name="label_3">
<property name="text">
<string>Maximum usable frequency:</string>
</property>
</widget>
</item>
<item row="3" column="1">
<widget class="SIUnitEdit" name="maxFreq">
<property name="enabled">
<bool>false</bool>
</property>
</widget>
</item>
</layout>
</item>
<item>
<widget class="QDialogButtonBox" name="buttonBox">
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
<property name="standardButtons">
<set>QDialogButtonBox::Cancel|QDialogButtonBox::Ok</set>
</property>
</widget>
</item>
</layout>
</widget>
<customwidgets>
<customwidget>
<class>SIUnitEdit</class>
<extends>QLineEdit</extends>
<header>CustomWidgets/siunitedit.h</header>
</customwidget>
</customwidgets>
<resources/>
<connections>
<connection>
<sender>buttonBox</sender>
<signal>accepted()</signal>
<receiver>CalStandardLineEditDialog</receiver>
<slot>accept()</slot>
<hints>
<hint type="sourcelabel">
<x>248</x>
<y>254</y>
</hint>
<hint type="destinationlabel">
<x>157</x>
<y>274</y>
</hint>
</hints>
</connection>
<connection>
<sender>buttonBox</sender>
<signal>rejected()</signal>
<receiver>CalStandardLineEditDialog</receiver>
<slot>reject()</slot>
<hints>
<hint type="sourcelabel">
<x>316</x>
<y>260</y>
</hint>
<hint type="destinationlabel">
<x>286</x>
<y>274</y>
</hint>
</hints>
</connection>
</connections>
<buttongroups>
<buttongroup name="buttonGroup"/>
</buttongroups>
</ui>

106
Software/PC_Application/Calibration/CalStandardReflectEditDialog.ui Normal file
View File

@ -0,0 +1,106 @@
<?xml version="1.0" encoding="UTF-8"?>
<ui version="4.0">
<class>CalStandardReflectEditDialog</class>
<widget class="QDialog" name="CalStandardReflectEditDialog">
<property name="geometry">
<rect>
<x>0</x>
<y>0</y>
<width>291</width>
<height>109</height>
</rect>
</property>
<property name="windowTitle">
<string>Edit Reflect Standard</string>
</property>
<property name="modal">
<bool>true</bool>
</property>
<layout class="QVBoxLayout" name="verticalLayout">
<item>
<layout class="QFormLayout" name="formLayout_2">
<item row="0" column="0">
<widget class="QLabel" name="label">
<property name="text">
<string>Name:</string>
</property>
</widget>
</item>
<item row="0" column="1">
<widget class="QLineEdit" name="name"/>
</item>
</layout>
</item>
<item>
<layout class="QFormLayout" name="formLayout">
<item row="0" column="0">
<widget class="QLabel" name="label_2">
<property name="text">
<string>Reflection type:</string>
</property>
</widget>
</item>
<item row="0" column="1">
<widget class="QComboBox" name="type">
<item>
<property name="text">
<string>Open</string>
</property>
</item>
<item>
<property name="text">
<string>Short</string>
</property>
</item>
</widget>
</item>
</layout>
</item>
<item>
<widget class="QDialogButtonBox" name="buttonBox">
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
<property name="standardButtons">
<set>QDialogButtonBox::Cancel|QDialogButtonBox::Ok</set>
</property>
</widget>
</item>
</layout>
</widget>
<resources/>
<connections>
<connection>
<sender>buttonBox</sender>
<signal>accepted()</signal>
<receiver>CalStandardReflectEditDialog</receiver>
<slot>accept()</slot>
<hints>
<hint type="sourcelabel">
<x>248</x>
<y>254</y>
</hint>
<hint type="destinationlabel">
<x>157</x>
<y>274</y>
</hint>
</hints>
</connection>
<connection>
<sender>buttonBox</sender>
<signal>rejected()</signal>
<receiver>CalStandardReflectEditDialog</receiver>
<slot>reject()</slot>
<hints>
<hint type="sourcelabel">
<x>316</x>
<y>260</y>
</hint>
<hint type="destinationlabel">
<x>286</x>
<y>274</y>
</hint>
</hints>
</connection>
</connections>
</ui>

78
Software/PC_Application/Calibration/calibration.cpp
View File

@ -3,6 +3,7 @@
#include "CustomWidgets/informationbox.h"
#include "Util/app_common.h"
#include "unit.h"
#include "Util/util.h"
#include "LibreCAL/librecaldialog.h"
#include "Eigen/Dense"
@ -296,6 +297,7 @@ QString Calibration::TypeToString(Calibration::Type type)
case Type::None: return "None";
case Type::SOLT: return "SOLT";
case Type::ThroughNormalization: return "ThroughNormalization";
case Type::TRL: return "TRL";
case Type::Last: return "Invalid";
}
}
@ -560,6 +562,12 @@ void Calibration::edit()
}
}
}
// update calibration (may have changed due to deleted measurement)
if(canCompute(caltype)) {
compute(caltype);
} else {
deactivate();
}
updateMeasurementTable();
updateCalibrationList();
});
@ -605,6 +613,12 @@ void Calibration::edit()
for(auto s : selected) {
measurements[s.row()]->clearPoints();
}
// update calibration (may have changed due to deleted measurement)
if(canCompute(caltype)) {
compute(caltype);
} else {
deactivate();
}
updateMeasurementTable();
updateCalibrationList();
});
@ -647,8 +661,11 @@ CalibrationMeasurement::Base *Calibration::newMeasurement(CalibrationMeasurement
case CalibrationMeasurement::Base::Type::Open: m = new CalibrationMeasurement::Open(this); break;
case CalibrationMeasurement::Base::Type::Short: m = new CalibrationMeasurement::Short(this); break;
case CalibrationMeasurement::Base::Type::Load: m = new CalibrationMeasurement::Load(this); break;
case CalibrationMeasurement::Base::Type::SlidingLoad: m = new CalibrationMeasurement::SlidingLoad(this); break;
case CalibrationMeasurement::Base::Type::Reflect: m = new CalibrationMeasurement::Reflect(this); break;
case CalibrationMeasurement::Base::Type::Through: m = new CalibrationMeasurement::Through(this); break;
case CalibrationMeasurement::Base::Type::Isolation: m = new CalibrationMeasurement::Isolation(this); break;
case CalibrationMeasurement::Base::Type::Line: m = new CalibrationMeasurement::Line(this); break;
}
return m;
}
@ -679,13 +696,33 @@ Calibration::Point Calibration::computeSOLT(double f)
auto p = caltype.usedPorts[i];
auto _short = static_cast<CalibrationMeasurement::Short*>(findMeasurement(CalibrationMeasurement::Base::Type::Short, p));
auto open = static_cast<CalibrationMeasurement::Open*>(findMeasurement(CalibrationMeasurement::Base::Type::Open, p));
auto load = static_cast<CalibrationMeasurement::Load*>(findMeasurement(CalibrationMeasurement::Base::Type::Load, p));
auto s_m = _short->getMeasured(f);
auto o_m = open->getMeasured(f);
auto l_m = load->getMeasured(f);
auto s_c = _short->getActual(f);
auto o_c = open->getActual(f);
auto l_c = load->getActual(f);
complex<double> l_c, l_m;
auto slidingMeasurements = findMeasurements(CalibrationMeasurement::Base::Type::SlidingLoad, p);
if(slidingMeasurements.size() >= 3) {
// use sliding load
vector<complex<double>> slidingMeasured;
for(auto m : slidingMeasurements) {
auto slidingLoad = static_cast<CalibrationMeasurement::SlidingLoad*>(m);
auto value = slidingLoad->getMeasured(f);
if(isnan(abs(value))) {
throw runtime_error("missing sliding load measurement");
}
slidingMeasured.push_back(value);
}
// use center of measurement for ideal load measurement
l_m = Util::findCenterOfCircle(slidingMeasured);
// assume perfect sliding load
l_c = 0.0;
} else {
// use normal load standard
auto load = static_cast<CalibrationMeasurement::Load*>(findMeasurement(CalibrationMeasurement::Base::Type::Load, p));
l_c = load->getActual(f);
l_m = load->getMeasured(f);
}
auto denom = l_c * o_c * (o_m - l_m) + l_c * s_c * (l_m - s_m) + o_c * s_c * (s_m - o_m);
point.D[i] = (l_c * o_m * (s_m * (o_c - s_c) + l_m * s_c) - l_c * o_c * l_m * s_m + o_c * l_m * s_c * (s_m - o_m)) / denom;
point.S[i] = (l_c * (o_m - s_m) + o_c * (s_m - l_m) + s_c * (l_m - o_m)) / denom;
@ -1203,8 +1240,14 @@ bool Calibration::canCompute(Calibration::CalType type, double *startFreq, doubl
for(auto p : type.usedPorts) {
required.push_back({.type = CalibrationMeasurement::Base::Type::Short, .port1 = p});
required.push_back({.type = CalibrationMeasurement::Base::Type::Open, .port1 = p});
if(findMeasurements(CalibrationMeasurement::Base::Type::SlidingLoad, p).size() >= 3) {
// got enough sliding load measurements, use these
required.push_back({.type = CalibrationMeasurement::Base::Type::SlidingLoad, .port1 = p});
} else {
// not enough sliding load measurement, use normal load
required.push_back({.type = CalibrationMeasurement::Base::Type::Load, .port1 = p});
}
}
// through measurements between all ports
for(int i=1;i<=type.usedPorts.size();i++) {
for(int j=i+1;j<=type.usedPorts.size();j++) {
@ -1220,6 +1263,19 @@ bool Calibration::canCompute(Calibration::CalType type, double *startFreq, doubl
}
}
break;
case Type::TRL:
// Reflect measurement for every port
for(auto p : type.usedPorts) {
required.push_back({.type = CalibrationMeasurement::Base::Type::Reflect, .port1 = p});
}
// through and line measurements between all ports
for(int i=1;i<=type.usedPorts.size();i++) {
for(int j=i+1;j<=type.usedPorts.size();j++) {
required.push_back({.type = CalibrationMeasurement::Base::Type::Through, .port1 = i, .port2 = j});
required.push_back({.type = CalibrationMeasurement::Base::Type::Line, .port1 = i, .port2 = j});
}
}
break;
}
if(required.size() > 0) {
vector<CalibrationMeasurement::Base*> foundMeasurements;
@ -1278,6 +1334,7 @@ int Calibration::minimumPorts(Calibration::Type type)
switch(type) {
case Type::SOLT: return 1;
case Type::ThroughNormalization: return 2;
case Type::TRL: return 2;
}
return -1;
}
@ -1415,8 +1472,9 @@ bool Calibration::hasFrequencyOverlap(std::vector<CalibrationMeasurement::Base *
}
}
CalibrationMeasurement::Base *Calibration::findMeasurement(CalibrationMeasurement::Base::Type type, int port1, int port2)
std::vector<CalibrationMeasurement::Base *> Calibration::findMeasurements(CalibrationMeasurement::Base::Type type, int port1, int port2)
{
vector<CalibrationMeasurement::Base*> ret;
for(auto m : measurements) {
if(m->getType() != type) {
continue;
@ -1434,10 +1492,20 @@ CalibrationMeasurement::Base *Calibration::findMeasurement(CalibrationMeasuremen
}
}
// if we get here, we have a match
return m;
ret.push_back(m);
}
return ret;
}
CalibrationMeasurement::Base *Calibration::findMeasurement(CalibrationMeasurement::Base::Type type, int port1, int port2)
{
auto meas = findMeasurements(type, port1, port2);
if(meas.size() > 0) {
return meas[0];
} else {
return nullptr;
}
}
QString Calibration::CalType::getReadableDescription()
{

5
Software/PC_Application/Calibration/calibration.h
View File

@ -19,6 +19,7 @@ public:
None,
SOLT,
ThroughNormalization,
TRL,
Last,
};
class CalType {
@ -121,8 +122,10 @@ private:
void deleteMeasurements();
bool hasFrequencyOverlap(std::vector<CalibrationMeasurement::Base*> m, double *startFreq = nullptr, double *stopFreq = nullptr, int *points = nullptr);
// returns all measurements that match the paramaters
std::vector<CalibrationMeasurement::Base*> findMeasurements(CalibrationMeasurement::Base::Type type, int port1 = 0, int port2 = 0);
// returns the first measurement in the list that matches the parameters
CalibrationMeasurement::Base* findMeasurement(CalibrationMeasurement::Base::Type type, int port1 = 0, int port2 = 0);
CalibrationMeasurement::Base *newMeasurement(CalibrationMeasurement::Base::Type type);
class Point {

8
Software/PC_Application/Calibration/calibrationmeasurement.cpp
View File

@ -86,8 +86,11 @@ QString CalibrationMeasurement::Base::TypeToString(CalibrationMeasurement::Base:
case Type::Open: return "Open";
case Type::Short: return "Short";
case Type::Load: return "Load";
case Type::SlidingLoad: return "SlidingLoad";
case Type::Reflect: return "Reflect";
case Type::Through: return "Through";
case Type::Isolation: return "Isolation";
case Type::Line: return "Line";
case Type::Last: return "Invalid";
}
}
@ -678,3 +681,8 @@ std::vector<CalibrationMeasurement::Isolation::Point> CalibrationMeasurement::Is
{
return points;
}
QWidget *CalibrationMeasurement::SlidingLoad::createStandardWidget()
{
return new QLabel("Connect sliding load");
}

38
Software/PC_Application/Calibration/calibrationmeasurement.h
View File

@ -21,8 +21,11 @@ public:
Open,
Short,
Load,
SlidingLoad,
Reflect,
Through,
Isolation,
Line,
Last,
};
@ -138,6 +141,31 @@ public:
virtual Type getType() override {return Type::Load;}
};
class SlidingLoad : public OnePort
{
Q_OBJECT
public:
SlidingLoad(Calibration *cal) :
OnePort(cal){setFirstSupportedStandard();}
virtual QWidget* createStandardWidget() override;
virtual std::set<CalStandard::Virtual::Type> supportedStandardTypes() override {return {};}
virtual Type getType() override {return Type::SlidingLoad;}
};
class Reflect : public OnePort
{
Q_OBJECT
public:
Reflect(Calibration *cal) :
OnePort(cal){setFirstSupportedStandard();}
virtual std::set<CalStandard::Virtual::Type> supportedStandardTypes() override {return {CalStandard::Virtual::Type::Open,
CalStandard::Virtual::Type::Short, CalStandard::Virtual::Type::Reflect};}
virtual Type getType() override {return Type::Reflect;}
};
class TwoPort : public Base
{
Q_OBJECT
@ -234,5 +262,15 @@ protected:
std::vector<Point> points;
};
class Line : public TwoPort
{
Q_OBJECT
public:
Line(Calibration *cal) :
TwoPort(cal){setFirstSupportedStandard();}
virtual std::set<CalStandard::Virtual::Type> supportedStandardTypes() override {return {CalStandard::Virtual::Type::Line};}
virtual Type getType() override {return Type::Line;}
};
}
#endif // CALIBRATIONMEASUREMENT_H

126
Software/PC_Application/Calibration/calstandard.cpp
View File

@ -2,7 +2,9 @@
#include "ui_CalStandardOpenEditDialog.h"
#include "ui_CalStandardShortEditDialog.h"
#include "ui_CalStandardLoadEditDialog.h"
#include "ui_CalStandardReflectEditDialog.h"
#include "ui_CalStandardThroughEditDialog.h"
#include "ui_CalStandardLineEditDialog.h"
#include "unit.h"
#include "Util/util.h"
@ -23,8 +25,9 @@ Virtual *Virtual::create(Virtual::Type type)
case Type::Open: return new Open;
case Type::Short: return new Short;
case Type::Load: return new Load;
case Type::Reflect: return new Reflect;
case Type::Through: return new Through;
case Type::Line: // TODO
case Type::Line: return new Line;
case Type::Last:
break;
}
@ -46,6 +49,7 @@ QString Virtual::TypeToString(Virtual::Type type)
case Type::Open: return "Open";
case Type::Short: return "Short";
case Type::Load: return "Load";
case Type::Reflect: return "Reflect";
case Type::Through: return "Through";
case Type::Line: return "Line";
case Type::Last: return "Invalid";
@ -677,3 +681,123 @@ void Through::fromJSON(nlohmann::json j)
delay = j.value("delay", 0.0);
loss = j.value("loss", 0.0);
}
Reflect::Reflect()
{
isShort = true;
}
std::complex<double> Reflect::toS11(double freq)
{
Q_UNUSED(freq)
return std::numeric_limits<complex<double>>::quiet_NaN();
}
void Reflect::edit(std::function<void ()> finishedCallback)
{
auto d = new QDialog;
auto ui = new Ui::CalStandardReflectEditDialog;
ui->setupUi(d);
ui->name->setText(name);
ui->type->setCurrentIndex(isShort ? 1 : 0);
QObject::connect(d, &QDialog::accepted, [=](){
name = ui->name->text();
isShort = ui->type->currentIndex() == 1;
if(finishedCallback) {
finishedCallback();
}
});
d->show();
}
nlohmann::json Reflect::toJSON()
{
auto j = OnePort::toJSON();
j["isShort"] = isShort;
return j;
}
void Reflect::fromJSON(nlohmann::json j)
{
OnePort::fromJSON(j);
isShort = j.value("isShort", true);
}
Line::Line()
{
Z0 = 50.0;
setDelay(0.0);
}
Sparam Line::toSparam(double freq)
{
Sparam ret;
ret.m11 = numeric_limits<complex<double>>::quiet_NaN();
ret.m12 = numeric_limits<complex<double>>::quiet_NaN();
ret.m21 = numeric_limits<complex<double>>::quiet_NaN();
ret.m22 = numeric_limits<complex<double>>::quiet_NaN();
return ret;
}
void Line::edit(std::function<void ()> finishedCallback)
{
auto d = new QDialog;
auto ui = new Ui::CalStandardLineEditDialog;
ui->setupUi(d);
ui->name->setText(name);
ui->Z0->setUnit("Ω");
ui->Z0->setPrecision(2);
ui->Z0->setValue(Z0);
ui->delay->setValue(delay);
ui->delay->setUnit("s");
ui->delay->setPrefixes("pnum ");
ui->delay->setPrecision(4);
ui->minFreq->setUnit("Hz");
ui->minFreq->setPrecision(3);
ui->minFreq->setPrefixes(" kMG");
ui->minFreq->setValue(minFreq);
ui->maxFreq->setUnit("Hz");
ui->maxFreq->setPrecision(3);
ui->maxFreq->setPrefixes(" kMG");
ui->maxFreq->setValue(maxFreq);
QObject::connect(ui->delay, &SIUnitEdit::valueChanged, [=](double val){
ui->minFreq->setValue(1.0 / val * 20 / 360);
ui->maxFreq->setValue(1.0 / val * 160 / 360);
});
QObject::connect(d, &QDialog::accepted, [=](){
name = ui->name->text();
Z0 = ui->Z0->value();
setDelay(ui->delay->value());
if(finishedCallback) {
finishedCallback();
}
});
d->show();
}
nlohmann::json Line::toJSON()
{
auto j = TwoPort::toJSON();
j["delay"] = delay;
return j;
}
void Line::fromJSON(nlohmann::json j)
{
TwoPort::fromJSON(j);
setDelay(j.value("delay", 0.0));
}
void Line::setDelay(double delay)
{
this->delay = delay;
minFreq = 1.0 / delay * 20 / 360;
maxFreq = 1.0 / delay * 160 / 360;
}

35
Software/PC_Application/Calibration/calstandard.h
View File

@ -20,6 +20,7 @@ public:
Open,
Short,
Load,
Reflect,
Through,
Line,
Last
@ -126,6 +127,23 @@ private:
bool Cfirst;
};
class Reflect : public OnePort
{
public:
Reflect();
Reflect(QString name, bool isShort)
: OnePort(name), isShort(isShort){}
virtual std::complex<double> toS11(double freq) override;
virtual void edit(std::function<void(void)> finishedCallback = nullptr) override;
virtual Type getType() override {return Type::Reflect;}
virtual nlohmann::json toJSON() override;
virtual void fromJSON(nlohmann::json j) override;
private:
bool isShort;
};
class TwoPort : public Virtual
{
public:
@ -161,6 +179,23 @@ private:
double Z0, delay, loss;
};
class Line : public TwoPort
{
public:
Line();
Line(QString name, double Z0, double delay)
: TwoPort(name), Z0(Z0), delay(delay){}
virtual Sparam toSparam(double freq) override;
virtual void edit(std::function<void(void)> finishedCallback = nullptr) override;
virtual Type getType() override {return Type::Line;}
virtual nlohmann::json toJSON() override;
virtual void fromJSON(nlohmann::json j) override;
private:
void setDelay(double delay);
double Z0, delay;
};
}
#endif // CALSTANDARD_H

2
Software/PC_Application/LibreVNA-GUI.pro
View File

@ -278,8 +278,10 @@ osx:LIBS += $(shell pkg-config --libs libusb-1.0)
QT += widgets network
FORMS += \
Calibration/CalStandardLineEditDialog.ui \
Calibration/CalStandardLoadEditDialog.ui \
Calibration/CalStandardOpenEditDialog.ui \
Calibration/CalStandardReflectEditDialog.ui \
Calibration/CalStandardShortEditDialog.ui \
Calibration/CalStandardThroughEditDialog.ui \
Calibration/LibreCAL/librecaldialog.ui \

79
Software/PC_Application/Util/util.cpp
View File

@ -90,3 +90,82 @@ unsigned long long Util::random(unsigned long long max)
std::uniform_int_distribution<unsigned long long> distribute(0, max);
return distribute(generator);
}
std::complex<double> Util::findCenterOfCircle(const std::vector<std::complex<double> > &points)
{
int i,iter,IterMAX=99;
double Xi,Yi,Zi;
double Mz,Mxy,Mxx,Myy,Mxz,Myz,Mzz,Cov_xy,Var_z;
double A0,A1,A2,A22;
double Dy,xnew,x,ynew,y;
double DET,Xcenter,Ycenter;
// find means
double meanX = 0.0, meanY = 0.0;
for(auto p : points) {
meanX += p.real();
meanY += p.imag();
}
meanX /= points.size();
meanY /= points.size();
// computing moments
Mxx=Myy=Mxy=Mxz=Myz=Mzz=0.;
for (i=0; i<(int) points.size(); i++)
{
Xi = points[i].real() - meanX; // centered x-coordinates
Yi = points[i].imag() - meanY; // centered y-coordinates
Zi = Xi*Xi + Yi*Yi;
Mxy += Xi*Yi;
Mxx += Xi*Xi;
Myy += Yi*Yi;
Mxz += Xi*Zi;
Myz += Yi*Zi;
Mzz += Zi*Zi;
}
Mxx /= points.size();
Myy /= points.size();
Mxy /= points.size();
Mxz /= points.size();
Myz /= points.size();
Mzz /= points.size();
// computing the coefficients of the characteristic polynomial
Mz = Mxx + Myy;
Cov_xy = Mxx*Myy - Mxy*Mxy;
Var_z = Mzz - Mz*Mz;
A2 = 4.0*Cov_xy - 3.0*Mz*Mz - Mzz;
A1 = Var_z*Mz + 4.0*Cov_xy*Mz - Mxz*Mxz - Myz*Myz;
A0 = Mxz*(Mxz*Myy - Myz*Mxy) + Myz*(Myz*Mxx - Mxz*Mxy) - Var_z*Cov_xy;
A22 = A2 + A2;
// finding the root of the characteristic polynomial
// using Newton's method starting at x=0
// (it is guaranteed to converge to the right root)
for (x=0.,y=A0,iter=0; iter<IterMAX; iter++) // usually, 4-6 iterations are enough
{
Dy = A1 + x*(A22 + 16.*x*x);
xnew = x - y/Dy;
if ((xnew == x)||(!isfinite(xnew))) break;
ynew = A0 + xnew*(A1 + xnew*(A2 + 4.0*xnew*xnew));
if (abs(ynew)>=abs(y)) break;
x = xnew; y = ynew;
}
// computing paramters of the fitting circle
DET = x*x - x*Mz + Cov_xy;
Xcenter = (Mxz*(Myy - x) - Myz*Mxy)/DET/2.0;
Ycenter = (Myz*(Mxx - x) - Mxz*Mxy)/DET/2.0;
// assembling the output
return std::complex<double>(Xcenter + meanX, Ycenter + meanY);
}

2
Software/PC_Application/Util/util.h
View File

@ -77,6 +77,8 @@ namespace Util {
double distanceToLine(QPointF point, QPointF l1, QPointF l2, QPointF *closestLinePoint = nullptr, double *pointRatio = nullptr);
unsigned long long random(unsigned long long max);
std::complex<double> findCenterOfCircle(const std::vector<std::complex<double>> &points);
}
#endif // UTILH_H