-
Notifications
You must be signed in to change notification settings - Fork 1
/
thresholdingdemo.cpp
110 lines (89 loc) · 3.53 KB
/
thresholdingdemo.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
//=============================================================================
/**
* Author: Jason Cipriani
* Website: https://github.com/JC3/ZScorePeakDetection
* License: https://github.com/JC3/ZScorePeakDetection/blob/master/LICENSE
*
* This is an implementation of the Robust Peak Detection Algorithm Using
* Z-Scores (Brackel, J.P.G. van) from https://stackoverflow.com/a/22640362
* with some experimental modifications by me:
*
* - Option to process data set in reverse.
* - Absolute minimum level of a peak.
*/
//=============================================================================
#include "thresholdingdemo.h"
#include <QMetaObject>
#include <cmath>
ThresholdingDemo::ThresholdingDemo (QObject *parent)
: QObject(parent)
{
qRegisterMetaType<ThresholdingDemo::Output>("ThresholdingDemo::Output");
}
void ThresholdingDemo::update () {
Output output;
threshold(params_, output, input_);
emit outputChanged(output);
}
template <typename T>
static QVector<T> & reverseif (QVector<T> &vec, bool cond) {
if (cond)
std::reverse(vec.begin(), vec.end());
return vec;
}
void ThresholdingDemo::threshold (const Params ¶ms,
Output &output,
QVector<float> input)
{
// adapted from https://stackoverflow.com/a/46956908/616460.
// there were some bugs in that implementation which have been corrected here,
// as well as porting std::vector to QVector and some API changes.
static const auto mean = [](const QVector<float> &vec, int from, int to) {
from = std::max(0, from);
float sum = 0;
for (int k = from; k <= to; ++ k) sum += vec[k];
return sum / (to - from + 1);
};
static const auto stdDev = [](const QVector<float> &vec, float mean, int from, int to) {
from = std::max(0, from);
float sum = 0;
for (int k = from; k <= to; ++ k) sum += std::pow(vec[k] - mean, 2);
return std::sqrt(sum / (to - from + 1));
};
reverseif(input, params.reverse);
int lag = params.lag;
float threshold = params.threshold;
float influence = params.influence;
output.clear();
if (input.size() <= lag + 2)
return;
QVector<int> signls(input.size(), 0.0);
QVector<float> avgFilter(input.size(), 0.0);
QVector<float> stdFilter(input.size(), 0.0);
QVector<float> filteredY = input;
avgFilter[lag] = mean(filteredY, 0, lag);
stdFilter[lag] = stdDev(filteredY, avgFilter[lag], 0, lag);
for (auto i = lag + 1; i < input.size(); i++) {
if (std::abs(input[i]) >= params.minlevel && std::abs(input[i] - avgFilter[i-1]) > threshold * stdFilter[i-1]) {
signls[i] = (input[i] > avgFilter[i-1]) ? 1 : -1;
filteredY[i] = influence* input[i] + (1 - influence) * filteredY[i-1];
} else {
signls[i] = 0; //# No signal
filteredY[i] = input[i];
}
avgFilter[i] = mean(filteredY, i-lag, i);
stdFilter[i] = stdDev(filteredY, avgFilter[i], i-lag, i);
}
output.input = reverseif(input, params.reverse);
output.params = params;
output.mean = reverseif(avgFilter, params.reverse);
output.stddev = reverseif(stdFilter, params.reverse);
output.signal = reverseif(signls, params.reverse);
output.outfrom = lag;
output.outto = input.size() - 1;
if (params.reverse) {
output.outfrom = input.size() - output.outfrom - 1;
output.outto = input.size() - output.outto - 1;
std::swap(output.outfrom, output.outto);
}
}