Revision 40:8f56ef28b0b1
| PeakInterpolator.cpp | ||
|---|---|---|
| 24 | 24 |
|
| 25 | 25 |
#include "PeakInterpolator.h" |
| 26 | 26 |
|
| 27 |
#include <iostream> |
|
| 28 |
|
|
| 27 | 29 |
static double cubicInterpolate(const double y[4], double x) |
| 28 | 30 |
{
|
| 29 | 31 |
double a0 = y[3] - y[2] - y[0] + y[1]; |
| ... | ... | |
| 40 | 42 |
double |
| 41 | 43 |
PeakInterpolator::findPeakLocation(const double *data, int size, int peakIndex) |
| 42 | 44 |
{
|
| 43 |
if (peakIndex < 2 || peakIndex > size - 3) {
|
|
| 45 |
std::cerr << "findPeakLocation: size " << size << ", peakIndex " << peakIndex << std::endl; |
|
| 46 |
|
|
| 47 |
if (peakIndex < 1 || peakIndex > size - 2) {
|
|
| 48 |
std::cerr << "returning " << peakIndex << ", data too short" << std::endl; |
|
| 44 | 49 |
return peakIndex; |
| 45 | 50 |
} |
| 46 | 51 |
|
| ... | ... | |
| 53 | 58 |
y[0] = data[peakIndex-1]; |
| 54 | 59 |
y[1] = data[peakIndex]; |
| 55 | 60 |
y[2] = data[peakIndex+1]; |
| 56 |
y[3] = data[peakIndex+2]; |
|
| 61 |
if (peakIndex < size - 2) {
|
|
| 62 |
y[3] = data[peakIndex+2]; |
|
| 63 |
} else {
|
|
| 64 |
y[3] = y[2]; |
|
| 65 |
} |
|
| 66 |
std::cerr << "a y: " << y[0] << " " << y[1] << " " << y[2] << " " << y[3] << std::endl; |
|
| 57 | 67 |
for (int i = 0; i < divisions; ++i) {
|
| 58 | 68 |
double probe = double(i) / double(divisions); |
| 59 | 69 |
double value = cubicInterpolate(y, probe); |
| 70 |
std::cerr << "probe = " << probe << ", value = " << value << " for location " << peakIndex + probe << std::endl; |
|
| 60 | 71 |
if (value > maxval) {
|
| 61 | 72 |
maxval = value; |
| 62 | 73 |
location = peakIndex + probe; |
| ... | ... | |
| 66 | 77 |
y[3] = y[2]; |
| 67 | 78 |
y[2] = y[1]; |
| 68 | 79 |
y[1] = y[0]; |
| 69 |
y[0] = data[peakIndex-2]; |
|
| 80 |
if (peakIndex > 1) {
|
|
| 81 |
y[0] = data[peakIndex-2]; |
|
| 82 |
} else {
|
|
| 83 |
y[0] = y[1]; |
|
| 84 |
} |
|
| 85 |
std::cerr << "b y: " << y[0] << " " << y[1] << " " << y[2] << " " << y[3] << std::endl; |
|
| 70 | 86 |
for (int i = 0; i < divisions; ++i) {
|
| 71 | 87 |
double probe = double(i) / double(divisions); |
| 72 | 88 |
double value = cubicInterpolate(y, probe); |
| 89 |
std::cerr << "probe = " << probe << ", value = " << value << " for location " << peakIndex - 1 + probe << std::endl; |
|
| 73 | 90 |
if (value > maxval) {
|
| 74 | 91 |
maxval = value; |
| 75 | 92 |
location = peakIndex - 1 + probe; |
| 76 | 93 |
} |
| 77 | 94 |
} |
| 78 | 95 |
|
| 96 |
std::cerr << "returning " << location << std::endl; |
|
| 97 |
|
|
| 79 | 98 |
return location; |
| 80 | 99 |
} |
| 81 | 100 |
|
| test/TestPeakInterpolator.cpp | ||
|---|---|---|
| 39 | 39 |
BOOST_CHECK_EQUAL(result, 1.0); |
| 40 | 40 |
} |
| 41 | 41 |
|
| 42 |
BOOST_AUTO_TEST_CASE(peakAtSample_N4) |
|
| 43 |
{
|
|
| 44 |
double data[] = { 0.0, 1.0, 2.0, 0.0 };
|
|
| 45 |
PeakInterpolator p; |
|
| 46 |
double result = p.findPeakLocation(data, 4, 2); |
|
| 47 |
//!!! Actually, this isn't certain at all I think? It could quite |
|
| 48 |
//!!! reasonably be smaller than 2.0 |
|
| 49 |
BOOST_CHECK_EQUAL(result, 2.0); |
|
| 50 |
} |
|
| 51 |
|
|
| 52 | 42 |
BOOST_AUTO_TEST_CASE(peakAtSample_N5) |
| 53 | 43 |
{
|
| 54 | 44 |
double data[] = { 0.0, 1.0, 2.0, 1.0, 0.0 };
|
| ... | ... | |
| 73 | 63 |
BOOST_CHECK_EQUAL(result, 3.0); |
| 74 | 64 |
} |
| 75 | 65 |
|
| 66 |
BOOST_AUTO_TEST_CASE(start) |
|
| 67 |
{
|
|
| 68 |
// Can't meaningfully interpolate if we're identifying element 0 |
|
| 69 |
// as the peak |
|
| 70 |
double data[] = { 1.0, 1.0, 0.0, 0.0 };
|
|
| 71 |
PeakInterpolator p; |
|
| 72 |
double result = p.findPeakLocation(data, 4, 0); |
|
| 73 |
BOOST_CHECK_EQUAL(result, 0.0); |
|
| 74 |
} |
|
| 75 |
|
|
| 76 |
BOOST_AUTO_TEST_CASE(end) |
|
| 77 |
{
|
|
| 78 |
// Likewise for the final element |
|
| 79 |
double data[] = { 0.0, 0.0, 1.0, 1.0 };
|
|
| 80 |
PeakInterpolator p; |
|
| 81 |
double result = p.findPeakLocation(data, 4, 3); |
|
| 82 |
BOOST_CHECK_EQUAL(result, 3.0); |
|
| 83 |
} |
|
| 84 |
|
|
| 85 |
BOOST_AUTO_TEST_CASE(longHalfway) |
|
| 86 |
{
|
|
| 87 |
double data[] = { 1.0, 1.0, 1.0, 2.0, 2.0, 1.0, 1.0, 1.0 };
|
|
| 88 |
PeakInterpolator p; |
|
| 89 |
double result = p.findPeakLocation(data, 8, 4); |
|
| 90 |
BOOST_CHECK_EQUAL(result, 3.5); |
|
| 91 |
} |
|
| 92 |
|
|
| 93 |
BOOST_AUTO_TEST_CASE(shortHalfway) |
|
| 94 |
{
|
|
| 95 |
double data[] = { 1.0, 2.0, 2.0, 1.0 };
|
|
| 96 |
PeakInterpolator p; |
|
| 97 |
double result = p.findPeakLocation(data, 4, 1); |
|
| 98 |
BOOST_CHECK_EQUAL(result, 1.5); |
|
| 99 |
} |
|
| 100 |
|
|
| 76 | 101 |
BOOST_AUTO_TEST_SUITE_END() |
| 77 | 102 |
|
Also available in: Unified diff