Revision 42:45b4401136f6

View differences:

PeakInterpolator.cpp
40 40 
}
41 41 

  		




  42
  42
  
    
double


  		




  
  43
  
    
PeakInterpolator::findPeakLocation(const double *data, int size)


  		




  
  44
  
    
{


  		




  
  45
  
    
    double maxval;


  		




  
  46
  
    
    int maxidx = 0;


  		




  
  47
  
    
    int i;


  		




  
  48
  
    
    for (i = 0; i < size; ++i) {


  		




  
  49
  
    
        if (i == 0 || data[i] > maxval) {


  		




  
  50
  
    
            maxval = data[i];


  		




  
  51
  
    
            maxidx = i;


  		




  
  52
  
    
        }


  		




  
  53
  
    
    }


  		




  
  54
  
    
    return findPeakLocation(data, size, maxidx);


  		




  
  55
  
    
}


  		




  
  56
  
    

  		




  
  57
  
    
double


  		




  43
  58
  
    
PeakInterpolator::findPeakLocation(const double *data, int size, int peakIndex)


  		




  44
  59
  
    
{


  		




  45
  
  
    
    std::cerr << "findPeakLocation: size " << size << ", peakIndex " << peakIndex << std::endl;


  		




  
  60
  
    
//    std::cerr << "findPeakLocation: size " << size << ", peakIndex " << peakIndex << std::endl;


  		




  46
  61
  
    

  		




  47
  62
  
    
    if (peakIndex < 1 || peakIndex > size - 2) {


  		




  48
  
  
    
        std::cerr << "returning " << peakIndex << ", data too short" << std::endl;


  		




  
  63
  
    
//        std::cerr << "returning " << peakIndex << ", data too short" << std::endl;


  		




  49
  64
  
    
        return peakIndex;


  		




  50
  65
  
    
    }


  		




  51
  66
  
    

  		




  ......




  63
  78
  
    
    } else {


  		




  64
  79
  
    
        y[3] = y[2];


  		




  65
  80
  
    
    }


  		




  66
  
  
    
    std::cerr << "a y: " << y[0] << " " << y[1] << " " << y[2] << " " << y[3] << std::endl;


  		




  
  81
  
    
//    std::cerr << "a y: " << y[0] << " " << y[1] << " " << y[2] << " " << y[3] << std::endl;


  		




  67
  82
  
    
    for (int i = 0; i < divisions; ++i) {


  		




  68
  83
  
    
        double probe = double(i) / double(divisions);


  		




  69
  84
  
    
        double value = cubicInterpolate(y, probe);


  		




  70
  
  
    
        std::cerr << "probe = " << probe << ", value = " << value << " for location " << peakIndex + probe << std::endl;


  		




  
  85
  
    
//        std::cerr << "probe = " << probe << ", value = " << value << " for location " << peakIndex + probe << std::endl;


  		




  71
  86
  
    
        if (value > maxval) {


  		




  72
  87
  
    
            maxval = value; 


  		




  73
  88
  
    
            location = peakIndex + probe;


  		




  ......




  82
  97
  
    
    } else {


  		




  83
  98
  
    
        y[0] = y[1];


  		




  84
  99
  
    
    }


  		




  85
  
  
    
    std::cerr << "b y: " << y[0] << " " << y[1] << " " << y[2] << " " << y[3] << std::endl;


  		




  
  100
  
    
//    std::cerr << "b y: " << y[0] << " " << y[1] << " " << y[2] << " " << y[3] << std::endl;


  		




  86
  101
  
    
    for (int i = 0; i < divisions; ++i) {


  		




  87
  102
  
    
        double probe = double(i) / double(divisions);


  		




  88
  103
  
    
        double value = cubicInterpolate(y, probe);


  		




  89
  
  
    
        std::cerr << "probe = " << probe << ", value = " << value << " for location " << peakIndex - 1 + probe << std::endl;


  		




  
  104
  
    
//        std::cerr << "probe = " << probe << ", value = " << value << " for location " << peakIndex - 1 + probe << std::endl;


  		




  90
  105
  
    
        if (value > maxval) {


  		




  91
  106
  
    
            maxval = value; 


  		




  92
  107
  
    
            location = peakIndex - 1 + probe;


  		




  93
  108
  
    
        }


  		




  94
  109
  
    
    }


  		




  95
  110
  
    

  		




  96
  
  
    
    std::cerr << "returning " << location << std::endl;


  		




  
  111
  
    
//    std::cerr << "returning " << location << std::endl;


  		




  97
  112
  
    

  		




  98
  113
  
    
    return location;


  		




  99
  114
  
    
}


  		












  

    
      PeakInterpolator.h
    
  





  32
  32
  
    
    ~PeakInterpolator() { }


  		




  33
  33
  
    

  		




  34
  34
  
    
    /**


  		




  
  35
  
    
     * Return the interpolated location (i.e. possibly between sample


  		




  
  36
  
    
     * point indices) of the peak in the given sampled range.


  		




  
  37
  
    
     *


  		




  
  38
  
    
     * "The peak" is defined as the (approximate) location of the


  		




  
  39
  
    
     * maximum of a function interpolating between the points


  		




  
  40
  
    
     * neighbouring the sample index with the maximum value in the


  		




  
  41
  
    
     * range.


  		




  
  42
  
    
     *


  		




  
  43
  
    
     * If multiple local peak samples in the input range are equal,


  		




  
  44
  
    
     * i.e. there is more than one apparent peak in the range, the one


  		




  
  45
  
    
     * with the lowest index will be used. This is the case even if a


  		




  
  46
  
    
     * later peak would be of superior height after interpolation.


  		




  
  47
  
    
     */


  		




  
  48
  
    
    double findPeakLocation(const double *data, int size);


  		




  
  49
  
    

  		




  
  50
  
    
    /**


  		




  35
  51
  
    
     * Return the interpolated location (i.e. between sample point


  		




  36
  
  
    
     * indices) of the peak whose sample is found at peakIndex in a


  		




  37
  
  
    
     * series of size samples.


  		




  
  52
  
    
     * indices) of the peak whose nearest sample is found at peakIndex


  		




  
  53
  
    
     * in the given sampled range. This method allows you to specify


  		




  
  54
  
    
     * which peak to find, if local rather than global peaks are of


  		




  
  55
  
    
     * interest.


  		




  38
  56
  
    
     */


  		




  39
  57
  
    
    double findPeakLocation(const double *data, int size, int peakIndex);


  		




  40
  58
  
    
};


  		












  

    
      test/TestPeakInterpolator.cpp
    
  





  37
  37
  
    
    PeakInterpolator p;


  		




  38
  38
  
    
    double result = p.findPeakLocation(data, 3, 1);


  		




  39
  39
  
    
    BOOST_CHECK_EQUAL(result, 1.0);


  		




  
  40
  
    
    result = p.findPeakLocation(data, 3);


  		




  
  41
  
    
    BOOST_CHECK_EQUAL(result, 1.0);


  		




  40
  42
  
    
}


  		




  41
  43
  
    

  		




  42
  44
  
    
BOOST_AUTO_TEST_CASE(peakAtSample_N5)


  		




  ......




  45
  47
  
    
    PeakInterpolator p;


  		




  46
  48
  
    
    double result = p.findPeakLocation(data, 5, 2);


  		




  47
  49
  
    
    BOOST_CHECK_EQUAL(result, 2.0);


  		




  
  50
  
    
    result = p.findPeakLocation(data, 5);


  		




  
  51
  
    
    BOOST_CHECK_EQUAL(result, 2.0);


  		




  48
  52
  
    
}


  		




  49
  53
  
    

  		




  50
  54
  
    
BOOST_AUTO_TEST_CASE(flat)


  		




  ......




  53
  57
  
    
    PeakInterpolator p;


  		




  54
  58
  
    
    double result = p.findPeakLocation(data, 5, 2);


  		




  55
  59
  
    
    BOOST_CHECK_EQUAL(result, 2.0);


  		




  
  60
  
    
    result = p.findPeakLocation(data, 5);


  		




  
  61
  
    
    BOOST_CHECK_EQUAL(result, 0.0);


  		




  56
  62
  
    
}


  		




  57
  63
  
    

  		




  58
  64
  
    
BOOST_AUTO_TEST_CASE(multiPeak)


  		




  ......




  61
  67
  
    
    PeakInterpolator p;


  		




  62
  68
  
    
    double result = p.findPeakLocation(data, 5, 3);


  		




  63
  69
  
    
    BOOST_CHECK_EQUAL(result, 3.0);


  		




  
  70
  
    
    result = p.findPeakLocation(data, 5);


  		




  
  71
  
    
    BOOST_CHECK_EQUAL(result, 1.0);


  		




  64
  72
  
    
}


  		




  65
  73
  
    

  		




  66
  74
  
    
BOOST_AUTO_TEST_CASE(start)


  		




  ......




  80
  88
  
    
    PeakInterpolator p;


  		




  81
  89
  
    
    double result = p.findPeakLocation(data, 4, 3);


  		




  82
  90
  
    
    BOOST_CHECK_EQUAL(result, 3.0);


  		




  
  91
  
    
    // But when running without a peak location, we expect idx 2 to be


  		




  
  92
  
    
    // picked as peak, not idx 3, so that will result in interpolation


  		




  
  93
  
    
    result = p.findPeakLocation(data, 4);


  		




  
  94
  
    
    BOOST_CHECK(result > 2.0 && result < 3.0);


  		




  83
  95
  
    
}


  		




  84
  96
  
    

  		




  85
  97
  
    
BOOST_AUTO_TEST_CASE(longHalfway)


  		




  ......




  88
  100
  
    
    PeakInterpolator p;


  		




  89
  101
  
    
    double result = p.findPeakLocation(data, 8, 4);


  		




  90
  102
  
    
    BOOST_CHECK_EQUAL(result, 3.5);


  		




  
  103
  
    
    result = p.findPeakLocation(data, 8);


  		




  
  104
  
    
    BOOST_CHECK_EQUAL(result, 3.5);


  		




  91
  105
  
    
}


  		




  92
  106
  
    

  		




  93
  107
  
    
BOOST_AUTO_TEST_CASE(shortHalfway)


  		




  ......




  96
  110
  
    
    PeakInterpolator p;


  		




  97
  111
  
    
    double result = p.findPeakLocation(data, 4, 1);


  		




  98
  112
  
    
    BOOST_CHECK_EQUAL(result, 1.5);


  		




  
  113
  
    
    result = p.findPeakLocation(data, 4);


  		




  
  114
  
    
    BOOST_CHECK_EQUAL(result, 1.5);


  		




  99
  115
  
    
}


  		




  100
  116
  
    

  		




  101
  117
  
    
BOOST_AUTO_TEST_CASE(aboveHalfway)


  		




  ......




  104
  120
  
    
    PeakInterpolator p;


  		




  105
  121
  
    
    double result = p.findPeakLocation(data, 4, 2);


  		




  106
  122
  
    
    BOOST_CHECK(result > 1.5 && result < 2.0);


  		




  
  123
  
    
    result = p.findPeakLocation(data, 4);


  		




  
  124
  
    
    BOOST_CHECK(result > 1.5 && result < 2.0);


  		




  107
  125
  
    
}


  		




  108
  126
  
    

  		




  109
  127
  
    
BOOST_AUTO_TEST_CASE(belowHalfway)


  		




  ......




  112
  130
  
    
    PeakInterpolator p;


  		




  113
  131
  
    
    double result = p.findPeakLocation(data, 4, 1);


  		




  114
  132
  
    
    BOOST_CHECK(result > 1.0 && result < 1.5);


  		




  
  133
  
    
    result = p.findPeakLocation(data, 4);


  		




  
  134
  
    
    BOOST_CHECK(result > 1.0 && result < 1.5);


  		




  115
  135
  
    
}


  		




  116
  136
  
    

  		




  117
  137
  
    
BOOST_AUTO_TEST_SUITE_END()


  		












Also available in:  Unified diff