//
// 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 15:27:28 PDT 2001
// Filename:      ...sig/maint/code/Filter/Pitch/Pitch_ML.cpp
// Web Address:   http://sig.sapp.org/src/sigSignal/Pitch_ML.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.
//

#include "Pitch.h"
#include <stdlib.h>
#include "transforms.h"


//////////////////////////////
//
// ML -- Maximum Likelihood
// input parameters are:
//   signal      = 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)
//   tspectrum   = temporary holding space for abs spectrum.
//   refmatrix  = spectrum of test signal
//   csize       = number of columns in refmatrix
//

int Pitch::ML(double* signal, int size, int minIndex, int maxIndex, 
      double* spectrum, double* tspectrum, double** refmatrix, 
      int csize) {
   int i, j;

static int count = 0;
count++;

   // generate the magnitude spectrum
   absfft(tspectrum, signal, 1024); 

   /*  note: normalizing will not affect the outcome
   double norm = 0;
   for (i=0; i<512; i++) {
      if (tspectrum[i] > norm) {
         norm = tspectrum[i];
      }
   }
   for (i=0; i<512; i++) {
      tspectrum[i] = tspectrum[i] / norm;
   }
   */


   // matrix multiply the absolute spectrum with the test matrix
   // and keep track of maximum amplitude value

   int maxLocation = 0;

   for (j=0; j<csize; j++) {
      spectrum[j] = 0.0;
      for (i=0; i<512; i++) {
         spectrum[j] += refmatrix[j][i] * tspectrum[i];
      }
      if (spectrum[j] > spectrum[maxLocation]) {
         maxLocation = j;
      }
   }


   return maxLocation;
}



//////////////////////////////
//
// Pitch::createMLTextMatrix --
// 
// parameters:
//     matrix  =  matrix to store resulting spectra in.
//     rows    =  number of rows in matrix
//     columns =  number of columns in matrix
//     step    =  distance between columns. 1/2 = half step, 1/4 = quarter tone
//     harmonics = number of harmonic to consider
//     srate   =  sampling rate
//     window  =  signal window before transform (length = rows)
//

void Pitch::createMLRefMatrix(double** matrix, int rows, int columns, 
      double minFreq, double step, int harmonics, double srate) {
   int i, j;

   double signal[1024];
   double max = -1.0;
   double freq = 0.0;

   // generate a spectrum for each test signal
   // keeping track of the max amplitude for normalization
   for (j=0; j<columns; j++) {

      freq = minFreq * pow(2.0, j * step / 6.0);
      generateMLSignal(signal, 1024, harmonics, freq, srate);
      for (i=0; i<1024; i++) {
         signal[i] = window[i] * signal[i]; 
      }
      absfft(matrix[j], signal, 1024);
      for (i=0; i<512; i++) {
         if (matrix[j][i] > max) {
            max = matrix[j][i];
         }
      }
   }

   // normalize the spectrum set
   for (j=0; j<columns; j++) {
      for (i=0; i<rows; i++) {
         matrix[j][i] = matrix[j][i]/max;
      }
   }


int z;
for (z=0; z<512; z++) {
   cout << matrix[0][z] << "\taaa1" << endl;
}

for (z=0; z<512; z++) {
   cout << matrix[columns-1][z] << "\taaa2" << endl;
}


}



//////////////////////////////
//
// Pitch::generateMLSignal -- creates a harmonic signal for use in
//      createMLRefMatrix().
//

void Pitch::generateMLSignal(double* signal, int size, int harmonics, 
      double freq, double srate) {
   int i, h;

cout << "Frequency: " << freq << endl;

static int counter = 0;
counter++;

   for (i=0; i<size; i++) {
      signal[i] = 0.0;
      for (h=1; h<=harmonics; h++) {
         signal[i] += cos(freq * i * h * 2.0 * M_PI / srate);
      }
      signal[i] = signal[i] / harmonics;
   }

if (counter==15) {
for(int z=0; z<size; z++) {
cout << signal[z] << "\tgggsignal" << endl;
}
}

}


//////////////////////////////
//
// Pitch::setupML -- create reference matrix for ML algorithm.
//

void Pitch::setupML(void) {
   double minFreq = 146.83;
   double maxFreq = 440.0;
   double step = 1.0/16.0;
   int harmonics = 4;
   double srate = 44100;
   int columns = (int)(6.0 * log10(maxFreq/minFreq)/log10(2.0)/ step);
cout << "COLUMNS = " << columns << endl;

   mlrefmatrix.setSize(columns);
   for (int i=0; i<columns; i++) {
      mlrefmatrix[i] = new double[4096];
   }

   createMLRefMatrix(mlrefmatrix.getBase(), 4096, columns, 
      minFreq, step, harmonics, srate);

   mlinit = 1;
}





// md5sum: 33bf6257384fe8d5125b0419cea7862a Pitch_ML.cpp [20050403]