1
$\begingroup$

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);
}
$\endgroup$
1
  • $\begingroup$ It would be helpful if you could point to some paper which describes the algorithm. $\endgroup$
    – Phonon
    Nov 2, 2013 at 22:59

2 Answers 2

2
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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)
    Real spectrum[] = discardNegativeFrequencies(fullSpectrum);

    // 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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.