Mirror/aicodix___dsp
1
0
Fork 0
mirror of https://github.com/aicodix/dsp.git synced 2026-04-27 14:30:36 +00:00
Code Issues Projects Releases Packages Wiki Activity

added the repeated median estimator of Siegel

Browse source
This commit is contained in:
Ahmet Inan 2021-09-10 00:09:26 +02:00
commit e90a3b601b
2 changed files with 69 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

4
README.md
View file

@ -187,6 +187,10 @@ Implemented [Simple linear regression](https://en.wikipedia.org/wiki/Simple_line
The [Theil-Sen estimator](https://en.wikipedia.org/wiki/Theil%E2%80%93Sen_estimator) is a [robust](https://en.wikipedia.org/wiki/Robust_statistics) [line fitting](https://en.wikipedia.org/wiki/Line_fitting) algorithm.
### [repeated_median.hh](repeated_median.hh)
The [repeated median estimator](https://en.wikipedia.org/wiki/Repeated_median_regression) is a [robust](https://en.wikipedia.org/wiki/Robust_statistics) [line fitting](https://en.wikipedia.org/wiki/Line_fitting) algorithm.
### [complex.hh](complex.hh)
Faster alternative (no Inf/NaN handling) to the std::complex implementation.

65
repeated_median.hh Normal file
View file

@ -0,0 +1,65 @@
/*
Repeated median estimator of Siegel
Copyright 2021 Ahmet Inan <inan@aicodix.de>
*/
#pragma once
#include <algorithm>
namespace DSP {
template <typename TYPE, int SIZE>
class RepeatedMedianEstimator
{
TYPE inner_[SIZE-1], outer_[SIZE];
TYPE xint_, yint_, slope_;
public:
RepeatedMedianEstimator() : xint_(0), yint_(0), slope_(0) {}
void compute(TYPE *x, TYPE *y, int LEN)
{
if (LEN > SIZE)
LEN = SIZE;
for (int i = 0; i < LEN; ++i) {
int count = 0;
for (int j = 0; j < LEN; ++j)
if (x[j] != x[i])
inner_[count++] = (y[j] - y[i]) / (x[j] - x[i]);
std::nth_element(inner_, inner_+count/2, inner_+count);
outer_[i] = inner_[count/2];
}
std::nth_element(outer_, outer_+LEN/2, outer_+LEN);
slope_ = outer_[LEN/2];
for (int i = 0; i < LEN; ++i) {
int count = 0;
for (int j = 0; j < LEN; ++j)
if (x[j] != x[i])
inner_[count++] = (x[j]*y[i] - x[i]*y[j]) / (x[j] - x[i]);
std::nth_element(inner_, inner_+count/2, inner_+count);
outer_[i] = inner_[count/2];
}
std::nth_element(outer_, outer_+LEN/2, outer_+LEN);
yint_ = outer_[LEN/2];
xint_ = - yint_ / slope_;
}
TYPE xint()
{
return xint_;
}
TYPE slope()
{
return slope_;
}
TYPE yint()
{
return yint_;
}
TYPE operator () (TYPE x)
{
return yint_ + slope_ * x;
}
};
}