mirror of
https://github.com/ultimatepp/ultimatepp.git
synced 2026-05-16 06:05:58 -06:00
352 lines
No EOL
8.2 KiB
C++
352 lines
No EOL
8.2 KiB
C++
#include <plugin/Eigen/Eigen.h>
|
|
using namespace Eigen;
|
|
|
|
#include "ScatterDraw.h"
|
|
|
|
NAMESPACE_UPP
|
|
|
|
|
|
#define Membercall(fun) (this->*fun)
|
|
|
|
double DataSource::Min(Getdatafun getdata, int64& id) {
|
|
double minVal = -DOUBLE_NULL;
|
|
for (int64 i = 0; i < GetCount(); ++i) {
|
|
double d = Membercall(getdata)(i);
|
|
if (!IsNull(d) && minVal > d) {
|
|
minVal = d;
|
|
id = i;
|
|
}
|
|
}
|
|
if (minVal == -DOUBLE_NULL)
|
|
return Null;
|
|
return minVal;
|
|
}
|
|
|
|
double DataSource::Max(Getdatafun getdata, int64& id) {
|
|
double maxVal = DOUBLE_NULL;
|
|
for (int64 i = 0; i < GetCount(); ++i) {
|
|
double d = Membercall(getdata)(i);
|
|
if (!IsNull(d) && maxVal < d) {
|
|
maxVal = d;
|
|
id = i;
|
|
}
|
|
}
|
|
if (maxVal == DOUBLE_NULL)
|
|
return Null;
|
|
return maxVal;
|
|
}
|
|
|
|
double DataSource::Avg(Getdatafun getdata) {
|
|
double ret = 0;
|
|
int count = 0;
|
|
for (int64 i = 0; i < GetCount(); ++i) {
|
|
double d = Membercall(getdata)(i);
|
|
if (!IsNull(d)) {
|
|
ret += d;
|
|
count++;
|
|
}
|
|
}
|
|
if (count == 0)
|
|
return Null;
|
|
return ret/count;
|
|
}
|
|
|
|
double DataSource::Variance(Getdatafun getdata, double avg) {
|
|
if (IsNull(avg))
|
|
avg = Avg(getdata);
|
|
if (IsNull(avg))
|
|
return Null;
|
|
double ret = 0;
|
|
int count = 0;
|
|
for (int64 i = 0; i < GetCount(); ++i) {
|
|
double d = Membercall(getdata)(i);
|
|
if (!IsNull(d)) {
|
|
d -= avg;
|
|
ret += d*d;
|
|
count++;
|
|
}
|
|
}
|
|
if (count <= 0)
|
|
return Null;
|
|
return ret/(count - 1);
|
|
}
|
|
|
|
Vector<int64> DataSource::Envelope(Getdatafun getdataY, Getdatafun getdataX, double width, bool (*fun)(double a, double b)) {
|
|
Vector<int64> ret;
|
|
double width_2 = width/2.;
|
|
|
|
for (int i = 0; i < GetCount(); ++i) {
|
|
double y = Membercall(getdataY)(i);
|
|
double x = Membercall(getdataX)(i);
|
|
if (IsNull(x) || IsNull(y))
|
|
continue;
|
|
int numComparisons = 0;
|
|
for (int j = i-1; j >= 0; --j) {
|
|
double ynext = Membercall(getdataY)(j);
|
|
double xnext = Membercall(getdataX)(j);
|
|
if (IsNull(xnext) || IsNull(ynext))
|
|
continue;
|
|
if ((x - xnext) > width_2)
|
|
break;
|
|
if (!fun(y, ynext)) {
|
|
numComparisons = Null;
|
|
break;
|
|
}
|
|
numComparisons++;
|
|
}
|
|
if (IsNull(numComparisons))
|
|
continue;
|
|
for (int j = i+1; j < GetCount(); ++j) {
|
|
double ynext = Membercall(getdataY)(j);
|
|
double xnext = Membercall(getdataX)(j);
|
|
if (IsNull(xnext) || IsNull(ynext))
|
|
continue;
|
|
if ((xnext - x) > width_2)
|
|
break;
|
|
if (!fun(y, ynext)) {
|
|
numComparisons = Null;
|
|
break;
|
|
}
|
|
numComparisons++;
|
|
}
|
|
if (IsNull(numComparisons))
|
|
continue;
|
|
if (numComparisons > 2)
|
|
ret << i;
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
bool GreaterThan(double a, double b) {return a > b;}
|
|
bool LowerThan(double a, double b) {return a < b;}
|
|
|
|
Vector<int64> DataSource::UpperEnvelope(Getdatafun getdataY, Getdatafun getdataX, double width) {return Envelope(getdataY, getdataX, width, GreaterThan);}
|
|
Vector<int64> DataSource::LowerEnvelope(Getdatafun getdataY, Getdatafun getdataX, double width) {return Envelope(getdataY, getdataX, width, LowerThan);}
|
|
|
|
Vector<Pointf> DataSource::MovingAverage(Getdatafun getdataY, Getdatafun getdataX, double width) {
|
|
Vector<Pointf> ret;
|
|
double width_2 = width/2.;
|
|
|
|
for (int i = 0; i < GetCount(); ++i) {
|
|
double y = Membercall(getdataY)(i);
|
|
double x = Membercall(getdataX)(i);
|
|
if (IsNull(x) || IsNull(y))
|
|
continue;
|
|
int numAvg = 1;
|
|
double sum = y;
|
|
int j;
|
|
for (j = i-1; j >= 0; --j) {
|
|
double ynext = Membercall(getdataY)(j);
|
|
double xnext = Membercall(getdataX)(j);
|
|
if (IsNull(xnext) || IsNull(ynext))
|
|
continue;
|
|
if ((x - xnext) > width_2)
|
|
break;
|
|
sum += ynext;
|
|
numAvg++;
|
|
}
|
|
if (j < 0)
|
|
continue;
|
|
for (j = i+1; j < GetCount(); ++j) {
|
|
double ynext = Membercall(getdataY)(j);
|
|
double xnext = Membercall(getdataX)(j);
|
|
if (IsNull(xnext) || IsNull(ynext))
|
|
continue;
|
|
if ((xnext - x) > width_2)
|
|
break;
|
|
sum += ynext;
|
|
numAvg++;
|
|
}
|
|
if (j == GetCount())
|
|
continue;
|
|
ret << Pointf(x, sum/numAvg);
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
double DataSource::StdDev(Getdatafun getdata, double avg) {
|
|
double var = Variance(getdata, avg);
|
|
return IsNull(var) ? Null : sqrt(var);
|
|
}
|
|
|
|
double DataSource::SinEstim_Amplitude(double avg) {
|
|
return sqrt(2.*VarianceY(avg));
|
|
}
|
|
|
|
bool DataSource::SinEstim_FreqPhase(double &frequency, double &phase, double avg) {
|
|
if (GetCount() < 4)
|
|
return false;
|
|
if (IsNull(avg))
|
|
avg = AvgY();
|
|
int64 firstId;
|
|
for (firstId = 0; firstId < GetCount(); ++firstId)
|
|
if (!IsNull(x(firstId)) && !IsNull(y(firstId)))
|
|
break;
|
|
bool firstIsToPositive = (y(firstId) - avg) < 0;
|
|
bool isPossitive = !firstIsToPositive;
|
|
double T = 0;
|
|
int numT = 0;
|
|
double lastZero = Null;
|
|
double firstZero;
|
|
firstId++;
|
|
for (int64 id = firstId; id < GetCount(); ++id) {
|
|
if (IsNull(x(id)) || IsNull(y(id)))
|
|
continue;
|
|
if (((y(id) - avg) > 0) != isPossitive) {
|
|
isPossitive = !isPossitive;
|
|
double zero = x(id-1) - (y(id-1) - avg)*(x(id) - x(id-1))/(y(id) - y(id-1));
|
|
if (IsNull(lastZero))
|
|
firstZero = zero;
|
|
else {
|
|
T += zero - lastZero;
|
|
numT++;
|
|
}
|
|
lastZero = zero;
|
|
}
|
|
}
|
|
T = 2*T/numT;
|
|
frequency = 2*M_PI/T;
|
|
phase = -frequency*firstZero;
|
|
if (!firstIsToPositive)
|
|
phase += M_PI;
|
|
return true;
|
|
}
|
|
|
|
|
|
double CArray::znFixed(int n, int64 id) {
|
|
if (n == 0)
|
|
return zData[id];
|
|
NEVER();
|
|
return Null;
|
|
}
|
|
|
|
|
|
Vector<Pointf> DataSource::FFT(Getdatafun getdata, double tSample) {
|
|
int numData = int(GetCount());
|
|
VectorXd timebuf(numData);
|
|
int num = 0;
|
|
for (int i = 0; i < numData; ++i) {
|
|
double data = Membercall(getdata)(i);
|
|
if (!IsNull(data)) {
|
|
timebuf[i] = Membercall(getdata)(i);
|
|
num++;
|
|
}
|
|
}
|
|
Vector<Pointf> res;
|
|
if (num < 3)
|
|
return res;
|
|
timebuf.resize(num);
|
|
|
|
VectorXcd freqbuf;
|
|
try {
|
|
Eigen::FFT<double> fft;
|
|
fft.SetFlag(fft.HalfSpectrum);
|
|
fft.fwd(freqbuf, timebuf);
|
|
} catch(...) {
|
|
return res;
|
|
}
|
|
for (int i = 1; i < freqbuf.size(); ++i) {
|
|
double T = (tSample*numData)/i;
|
|
res << Pointf(T, 2*std::abs(freqbuf[i])/numData);
|
|
}
|
|
return res;
|
|
}
|
|
|
|
|
|
bool Spline::Load(Vector<double>& xdata, Vector<double>& ydata)
|
|
{
|
|
if(!IsInputSane(xdata, ydata))
|
|
return false;
|
|
|
|
int n = xdata.GetCount() - 1;
|
|
|
|
Vector<double> h;
|
|
for(int i = 0; i < n; ++i)
|
|
h << (xdata[i+1] - xdata[i]);
|
|
|
|
Vector<double> alpha;
|
|
for(int i = 0; i < n; ++i)
|
|
alpha << (3.*(ydata[i+1] - ydata[i])/h[i] - 3*(ydata[i]-ydata[i-1])/h[i-1]);
|
|
|
|
Vector<double> l, mu, z;
|
|
l.SetCount(n+1);
|
|
mu.SetCount(n+1);
|
|
z.SetCount(n+1);
|
|
l[0] = 1;
|
|
mu[0] = 0;
|
|
z[0] = 0;
|
|
|
|
for(int i = 1; i < n; ++i) {
|
|
l[i] = 2.*(xdata[i+1] - xdata[i-1]) - h[i-1]*mu[i-1];
|
|
mu[i] = h[i]/l[i];
|
|
z[i] = (alpha[i] - h[i-1]*z[i-1])/l[i];
|
|
}
|
|
|
|
l[n] = 1.;
|
|
z[n] = 0;
|
|
|
|
Vector<double> c;
|
|
c.SetCount(n+1);
|
|
c[n] = 0;
|
|
|
|
Vector<double> b, d;
|
|
b.SetCount(n);
|
|
d.SetCount(n);
|
|
for(int j = n-1; j >= 0; --j) {
|
|
c[j] = z[j] - mu[j] * c[j+1];
|
|
b[j] = (ydata[j+1] - ydata[j])/h[j] - h[j]*(c[j+1] + 2*c[j])/3.;
|
|
d[j] = (c[j+1] - c[j])/3./h[j];
|
|
}
|
|
|
|
coeff.SetCount(n);
|
|
for(int i = 0; i < n; ++i) {
|
|
coeff[i].a = ydata[i];
|
|
coeff[i].b = b[i];
|
|
coeff[i].c = c[i];
|
|
coeff[i].d = d[i];
|
|
coeff[i].x = xdata[i];
|
|
}
|
|
return true;
|
|
}
|
|
|
|
double Spline::GetY(double x)
|
|
{
|
|
int j;
|
|
for (j = 0; j < coeff.size(); j++) {
|
|
if(coeff[j].x > x) {
|
|
if(j == 0)
|
|
j++;
|
|
break;
|
|
}
|
|
}
|
|
j--;
|
|
|
|
double dx = x - coeff[j].x;
|
|
return coeff[j].a + coeff[j].b*dx + coeff[j].c*dx*dx + coeff[j].d*dx*dx*dx;
|
|
}
|
|
|
|
bool Spline::IsInputSane(Vector<double>& xdata, Vector<double>& ydata)
|
|
{
|
|
if(xdata.IsEmpty() || ydata.IsEmpty())
|
|
return false;
|
|
|
|
if(!IsSorted(xdata))
|
|
return false;
|
|
|
|
bool first = true;
|
|
double xold;
|
|
for(int i = 0; i < xdata.GetCount(); ++i) {
|
|
double x = xdata[i];
|
|
if(first)
|
|
xold = x;
|
|
else {
|
|
first = false;
|
|
if( fabs(x - xold) < 1)
|
|
return false;
|
|
xold = x;
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
|
|
END_UPP_NAMESPACE |