# How does one know how many times to downsample in Harmonic Product Spectrum?

I have been working on a Harmonic Product Spectrum Algorithm. Now, all the literature I've read about the subject tells me to downsample an N amount of times. How does one determine what this N value should be? Here is my implementation of Harmonic Product Spectrum so far. Feel free to tell me if I have gone wrong somewhere.

private int HarmonicProductSpectrum(Complex[] fftData, int n){
Complex[][] data = new Complex[n][fftData.length/n];
for(int i = 0; i<n; i++){
for(int j = 0; j<data[0].length; j++){
data[i][j] = fftData[j*(i+1)];
}
}
Complex[] result = new Complex[fftData.length/n];//Combines the arrays
for(int i = 0; i<result.length; i++){
Complex tmp = new Complex(1,0);
for(int j = 0; j<n; j++){//multiplies arrays together
tmp = tmp.times(data[j][i]);
}
result[i] = tmp;
}
//Calculates Maximum Magnitude of the array
double max = Double.MIN_VALUE;
int index = -1;
for(int i = 0; i<result.length; i++){
Complex c = result[i];
double tmp = c.getMagnitude();
if(tmp>max){
max = tmp;;
index = i;
}
}
return index*getFFTBinSize(fftData.length);
}

• It would be helpful if you could point to some paper which describes the algorithm. – Phonon Nov 2 '13 at 22:59

Depends on the pitch source spectrum, the lowest pitch possible, and the FFT length.

If N is too small, the algorithm might miss some higher harmonics that contain a significant fraction of a pitch spectrums energy. So you need to know how many overtones might be important in your particular pitch source.

However, if N is too large, multiple overtones of a single low pitch could end up in the same bin after downsampling the spectrum, confusing the results.

For very low pitches with extremely rich higher harmonics, these 2 constraints may overlap and thus indicate the need for a longer FFT window for HPS, or even the need for a completely different pitch estimation method.

Take a look at this paper on page 4. It has a great graphical representation of the algorithm. I'll try to write up some pseudocode here that you can adapt to your language (looks like Java to me) as necessary.

/* Harmonic Product spectrum PSEUDOCODE
* @param fftData Discrete Fourier transform of data
* @param N Number of times we downsample the spectrum to get HPS
*/
int HPS( Complex[] fftData, int N )
{
// Find magnitude of the FFT
Real fullSpectrum[] = absOfComplex(fftData);

// Keep only the positive frequencies (DC to Nyquist)

// Make a new array to store HPS
Real hps[] = copyOf(spectrum);

// Perfrom HPS:
// Go through each downsampling factor
for (int downsamplingFactor = 1; downsamplingFactor <= N; downsamplingFactor++)
{
// Go through samples of the downsampled signal and compute HPS at this iteration
for(int idx = 0; idx < spectrum.length()/downsamplingFactor; idx++)
{
hps[idx] *= spectrum[idx * downsamplingFactor];
}
}

return findIndexOfMax(hps);
}


This isn't bug-proof but it should get you started.