How to get fundamental frequency of a signal using autocorrelation?

Question

I am trying to get the fundamental frequency of a signal that only has a single pitch. I coded out the autocorrelation function using FFT and already got the autocorrelation result. Unfortunately, I don't know how to get the fundamental frequency from the autocorrelation result. Can someone help me? My code is below:

public double getPitch(double[] buffer, int firstSample, int lastSample, double sampleRate)
{
    int lengthOfFFTWindow=lastSample-firstSample;
    double[] input_buffer=new double[lengthOfFFTWindow];
    DoubleFFT_1D fft = new DoubleFFT_1D(lengthOfFFTWindow);
    double[] autocorrelation_values=new double[lengthOfFFTWindow];
    double[] fftData = new double[lengthOfFFTWindow * 2];
    double max=-1;
    double max_i=-1;
    //FFT on each sample in each window
    for (int i = 0; i < lengthOfFFTWindow; i++) {
        // copying audio data to the fft data buffer, imaginary part is 0
        fftData[2 * i] = buffer[i+firstSample];
        fftData[2 * i + 1] = 0;
    }
    fft.complexForward(fftData);
    for (int i = 0,j=0; i < fftData.length; i += 2,j++) {
        // complex numbers -> vectors, so we compute the length of the vector, which is sqrt(realpart^2+imaginarypart^2)
        autocorrelation_values[j] = Math.sqrt((fftData[i] * fftData[i]) + (fftData[i + 1] * fftData[i + 1]));
    }
    fft.complexInverse(fftData, false);
    for(int i=0;i<fftData.length;i++)
    {
        if(max<fftData[i])
        {
            max=fftData[i];
            max_i=i;
        }
    }
    return (max_i * 44100 )/ lengthOfFFTWindow;
}

Is it correct to return the maximum autocorrelation value divide by 2 as the fundamental frequency? I keep getting wrong answers when I do that.

EDIT: An example of the single pitch test file: http://dl.dropbox.com/u/24274475/testfile_A4.wav

Answer 1

Autocorrellation is one way to help find the dominant frequency of a signal, but I don't see what a FFT has to do with that. The autocorrellation will produce peaks with the period of any strong frequency components. If you then take the FFT of that to find the frequency of those peaks, you might as well take the FFT of the original signal in the first place.

Instead of showing us code, show us the data at various stages of your process. The details of the code are your issue, and are separate from the conceptual processes of going thru the various convolutions, filters, or whatever.

You say your signal only has a single pitch, meaning it's a pure sine wave. In that case, I really don't see the advantage of a autocorrellation pass. You can find the period directly by looking at the time between zero crossings.

In the past I've had to find the fundamental frequency of mostly repeating signals with significant noise on them. What I usually did was apply several stages of low pass filtering. One big advantage of digital filtering is that smaller signals out of a filter don't mean less signal to noise ratio as long as you keep adding the necessary bits at the low end. Using floating point, for example, does this automatically. You can then aggressively low pass filter a signal such that it would be only µV in analog, but still have the same meaningful bits left at the end.

Each LPF pass attenuates the harmonics relative to the fundamental. After enough passes, you are left with mostly the fundamental. Once you have attenuated the harmonics enough to guarantee only two zero crossings per cycle, you look at the zero crossing period, perhaps apply a little low pass filtering to successive ones, and infer the frequency from there.

Added:

Now that you have provided some data we can see what is really going on:

It looks like you have about a 440 Hz signal, but this is clearly far from a "single pitch" as the shape is far from a sine. Just from inspection we can see that the second harmonic is particularly strong. It may be so strong that this "note" is perceived to be 880 Hz instead of the fundamental of 440 Hz.

In this case, what is it you want the answer to be, 440 Hz or 880 Hz? With enough low pass filtering, eventually you get mostly the fundamental and measuring 440 Hz shouldn't be that hard. If you want the answer to be the possibly perceptual tone of 880 Hz, then things get a lot more complicated. One possibility would be to identify the fundamental in all cases. Once you have that, it's easy to find the relative amplitude of the first few harmonics. Then you can decide based on the strength of those harmonics whether you want to report one of them or the fundamental.

Answer 2

I see a couple of things that might be problems in your code.

First, to use the FFT to calculate an autocorrelation, there are three steps:

1. transform the input signal, \$F_R(f) = \mathrm{FFT}(x[t])\$
2. compute the power spectrum, \$S(f) = F_R(f) F_R^*(f)\$
3. inverse transform to get the autocorrelation, \$R[\tau]=\mathrm{IFFT}(S(f))\$

I see you doing step 1 and step 2, but then you do something completely different in step 3.

Note that if you did take the autocorrelation, the peaks would indicate the period of your signal, not the frequency.

What your code actually seems to be trying to do is take the FFT, then search the bins for the highest peak, and take that as the frequency of the signal. You also calculate the IFFT, but then throw away the result of that calculation.

A peak-search in the power spectrum can be used to estimate the signal frequency, but there you have another, more basic bug:

for(int i=0;i<autocorrelation_values.length;i++)
{
    if(max<autocorrelation_values[i])
    {
        max=autocorrelation_values[i];
        max_i=i;
    }
}
return max/2;

Your max is storing the maximum amplitude of the power spectrum values (which you confusingly named autocorrelation_values). max_i is storing the bin number where the maximum amplitude was found. You should be returning something based on max_i rather than max. You need to use the sampling rate and number of samples used in the FFT to scale the bin number into an actual frequency.

Edit I'd also recommend only searching the first half of the power spectrum for peaks.

for(int i=0; i < 0.5 * autocorrelation_values.length; i++)
   ...

The higher bins correspond to an alias of the spectrum in the negative frequencies, so don't contain any new information.

Answer 3

I want to find fundametal frequency with autocorrelation method. I read your code but I think that the last part of code could be like this:

int max=0; 

for (int i=0; i < autocorrelation_value.length; i++) {
  if (autocorrelation_value[i]>autocorrelation_value[max]) {
    max = i;
  }
  return 1/autocorrelation_value[max]
}

because the maximum of autocorrelation is the fundamental period and so in the "return" I calculate the inverse of autocorrelation's max.

And after I don't understand what means "max_i".

Answer 4

Why don't you do it the old tried and tested way by using a time-delayed version of the signal. Multiply the signal by delayed version of itself and store the result. For a given window in time containing (say) 1000 samples, you'll need to store the biggest result of the individual multiplications.

Then re-test with a range of different delays. The result that is numerically the greatest corresponds with the period of the signal. Be aware that once the delay has found the "period" of the signal, it will also find an "image" period at twice the delay (because the signal peaks coincide once more).

Alternatively use a zero-cross detector and measure the time. You said the signal has only a single pitch so this should be pefectly ok - watch out for zero-crossing noise by averaging a few cycles - the more the better.