//
// Programmer:    Patricio de la Cuadra <pdelac@ccrma.stanford.edu>
// Programmer:    Craig Stuart Sapp <craig@ccrma.stanford.edu>
// Creation Date: Spring 2000
// Last Modified: Mon Jul 16 11:18:09 PDT 2001
// Filename:      ...sig/maint/code/Filter/Pitch/Pitch_HPS.cpp
// Web Address:   http://sig.sapp.org/src/sigSignal/Pitch_HPS.cpp
// Documentation: http://sig.sapp.org/doc/classes/Pitch
// Syntax:        C++; close to C
//
// Reference:     Noll, Michael.  ``Pitch determination of human speech by the
//                harmonic product spectrum, the harmonic sum spectrum, and a 
//                maximum likelihood estimate,'' Proceedings of the Symposium 
//                on Computer Processing in Communications, April 8-10, 1969.
//
// HPS -- Harmonic Product Spectrum
// input parameters are:
//   paddedFrame = array of the windows signal frame already zero padded
//   size        = number of elements in the signal frame
//   minIndex    = minimum frequency index to assign as a pitch
//   maxIndex    = maximum frequency index to assign as a pitch
//   harmonics   = number of harmonics to consider
//   spectrum    = the spectrum values used to generate a pitch value.
//                 (same size as paddedFrame)
//

#include "Pitch.h"
#include "transforms.h"

#include <stdlib.h>


int Pitch::HPS(double* paddedFrame, int size, int minIndex, int maxIndex, 
      int harmonics, double* spectrum) {

   // generate the magnitude spectrum
   absfft(spectrum, paddedFrame, size); 

   // set the maximum index to search for a pitch
   int i, j, maxHIndex;
   maxHIndex = size/harmonics;
   if (maxIndex < maxHIndex) {
      maxHIndex = maxIndex;
   }

   // generate the Harmonic Product Spectrum values and keep track of the
   // maximum amplitude value to assign to a pitch.
   
   int maxLocation = minIndex;
   for (j=minIndex; j<=maxHIndex; j++) {
      for (i=1; i<=harmonics; i++) {  // i=1: double the fundamental
         spectrum[j] *= spectrum[j*i];
      }
      if (spectrum[j] > spectrum[maxLocation]) {
         maxLocation = j;
      }
   }

   // Correct for octave too high errors.  If the maximum subharmonic
   // of the measured maximum is approximately 1/2 of the maximum
   // measured frequency, AND if the ratio of the sub-harmonic
   // to the total maximum is greater than 0.2, THEN the pitch value
   // is assigned to the subharmonic.

   int max2 = minIndex;
   int maxsearch = maxLocation * 3 / 4;
   for (i=minIndex+1; i<maxsearch; i++) {
      if (spectrum[i] > spectrum[max2]) {
         max2 = i;
      }
   }
   if (abs(max2 * 2 - maxLocation) < 4) {
      if (spectrum[max2]/spectrum[maxLocation] > 0.2) {
         maxLocation = max2;
      }
   }

   return maxLocation;
}



// md5sum: 0b10d8d32c70ca28782fbf698c97119e Pitch_HPS.cpp [20050403]