Chris@39: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@39: /*
Chris@39:     This file is Copyright (c) 2012 Chris Cannam
Chris@39:   
Chris@39:     Permission is hereby granted, free of charge, to any person
Chris@39:     obtaining a copy of this software and associated documentation
Chris@39:     files (the "Software"), to deal in the Software without
Chris@39:     restriction, including without limitation the rights to use, copy,
Chris@39:     modify, merge, publish, distribute, sublicense, and/or sell copies
Chris@39:     of the Software, and to permit persons to whom the Software is
Chris@39:     furnished to do so, subject to the following conditions:
Chris@39: 
Chris@39:     The above copyright notice and this permission notice shall be
Chris@39:     included in all copies or substantial portions of the Software.
Chris@39: 
Chris@39:     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
Chris@39:     EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
Chris@39:     MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
Chris@39:     NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
Chris@39:     ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
Chris@39:     CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
Chris@39:     WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
Chris@39: */
Chris@39: 
Chris@39: #include "PeakInterpolator.h"
Chris@39: 
Chris@40: #include <iostream>
Chris@40: 
Chris@39: double
Chris@42: PeakInterpolator::findPeakLocation(const double *data, int size)
Chris@42: {
Chris@42:     double maxval;
Chris@42:     int maxidx = 0;
Chris@42:     int i;
Chris@42:     for (i = 0; i < size; ++i) {
Chris@42:         if (i == 0 || data[i] > maxval) {
Chris@42:             maxval = data[i];
Chris@42:             maxidx = i;
Chris@42:         }
Chris@42:     }
Chris@42:     return findPeakLocation(data, size, maxidx);
Chris@42: }
Chris@42: 
Chris@42: double
Chris@39: PeakInterpolator::findPeakLocation(const double *data, int size, int peakIndex)
Chris@39: {
Chris@43:     // after jos, 
Chris@43:     // https://ccrma.stanford.edu/~jos/sasp/Quadratic_Interpolation_Spectral_Peaks.html
Chris@40: 
Chris@40:     if (peakIndex < 1 || peakIndex > size - 2) {
Chris@39:         return peakIndex;
Chris@39:     }
Chris@39: 
Chris@43:     double alpha = data[peakIndex-1];
Chris@43:     double beta  = data[peakIndex];
Chris@43:     double gamma = data[peakIndex+1];
Chris@39: 
Chris@43:     double denom = (alpha - 2*beta + gamma);
Chris@39: 
Chris@43:     if (denom == 0) {
Chris@43:         // flat
Chris@43:         return peakIndex;
Chris@39:     }
Chris@39: 
Chris@43:     double p = ((alpha - gamma) / denom) / 2.0;
Chris@39: 
Chris@43:     return double(peakIndex) + p;
Chris@39: }
Chris@39: