Why are my frequency values for iPhone FFT incorrect?

Question

I've been trying to get exact frequencies using the FFT in Apple's Accelerate framework, but I'm having trouble working out why my values are off the true frequency.

I have been using this article http://www.dspdimension.com/admin/pitch-shifting-using-the-ft/ as the basis for my implementation, and after really struggling to get to the point I'm at now, I am totally stumped.

So far I've got audio in -> Hanning window -> FFT -> phase calculation -> weird final output. I'd think that there will be a problem with my maths somewhere, but I'm really out of ideas by now.

The outputs are a lot lower what they should be, e.g., I input 440Hz and it prints out 190Hz, or I input 880Hz and it prints out 400Hz. For the most part these results are consistent, but not always, and there doesn't seem to be any common factor between anything either...

Here is my code:

enum
{
    sampleRate  = 44100,
    osamp       = 4,
    samples     = 4096,
    range       = samples * 7 / 16,
    step        = samples / osamp
};

NSMutableArray *fftResults;

static COMPLEX_SPLIT    A;
static FFTSetup         setupReal;
static uint32_t         log2n, n, nOver2;
static int32_t          stride;

static float            expct = 2*M_PI*((double)step/(double)samples);
static float            phase1[range];
static float            phase2[range];
static float            dPhase[range];

- (void)fftSetup
{
    // Declaring integers
    log2n = 12;
    n = 1 << log2n;
    stride = 1;
    nOver2 = n / 2;

    // Allocating memory for complex vectors
    A.realp = (float *) malloc(nOver2 * sizeof(float));
    A.imagp = (float *) malloc(nOver2 * sizeof(float));

    // Allocating memory for FFT
    setupReal = vDSP_create_fftsetup(log2n, FFT_RADIX2);

    // Setting phase
    memset(phase2, 0, range * sizeof(float));

}
        // For each sample in buffer...
        for (int bufferCount = 0; bufferCount < audioBufferList.mNumberBuffers; bufferCount++)
        {
            // Declaring samples from audio buffer list
            SInt16 *samples = (SInt16*)audioBufferList.mBuffers[bufferCount].mData;


            // Creating Hann window function
            for (int i = 0; i < nOver2; i++)
            {
                double hannMultiplier = 0.5 * (1 - cos((2 * M_PI * i) / (nOver2 -  1)));

                // Applying window to each sample
                A.realp[i] = hannMultiplier * samples[i];
                A.imagp[i] = 0;
            }

            // Applying FFT
            vDSP_fft_zrip(setupReal, &A, stride, log2n, FFT_FORWARD);

            // Detecting phase
            vDSP_zvphas(&A, stride, phase1, stride, range);

            // Calculating phase difference
            vDSP_vsub(phase2, stride, phase1, stride, dPhase, stride, range);

            // Saving phase
            memcpy(phase2, phase1, range * sizeof(float));

            // Extracting DSP outputs
            for (size_t j = 0; j < nOver2; j++)
            {
                NSNumber *realNumbers = [NSNumber numberWithFloat:A.realp[j]];
                NSNumber *imagNumbers = [NSNumber numberWithFloat:A.imagp[j]];

                [real addObject:realNumbers];
                [imag addObject:imagNumbers];

            }

            // Combining real and imaginary parts
            [resultsCombined addObject:real];
            [resultsCombined addObject:imag];

            // Filling FFT output array
            [fftResults addObject:resultsCombined];
        }

    }


    int fftCount = [fftResults count];
    NSLog(@"FFT Count: %d",fftCount);

    // For each FFT...
    for (int i = 0; i < fftCount; i++)
    {
        // Declaring integers for peak detection
        float peak = 0;
        float binNumber = 0;

        // Declaring integers for phase detection
        float deltaPhase;

        static float trueFrequency[range];


        for (size_t j = 1; j < range; j++)
        {
            // Calculating bin magnitiude
            float realVal   = [[[[fftResults objectAtIndex:i] objectAtIndex:0] objectAtIndex:j] floatValue];
            float imagVal   = [[[[fftResults objectAtIndex:i] objectAtIndex:1] objectAtIndex:j] floatValue];
            float magnitude = sqrtf(realVal*realVal + imagVal*imagVal);

            // Peak detection
            if (magnitude > peak)
            {
                peak = magnitude;
                binNumber = (float)j;
            }

            // Getting phase difference
            deltaPhase = dPhase[j];

            // Subtract expected difference
            deltaPhase -= (float)j*expct;

            // Map phase difference into +/- pi interval
            int qpd = deltaPhase / M_PI;

            if (qpd >= 0)
                qpd += qpd&1;
            else
                qpd -= qpd&1;

            deltaPhase -= M_PI * (float)qpd;

            // Getting bin deviation from +/i interval
            float deltaFrequency = osamp * deltaPhase / (2 * M_PI);

            // Calculating true frequency at the j-th partial
            trueFrequency[j] = (j * (sampleRate/samples)) + (deltaFrequency * (sampleRate/samples));

        }

        UInt32 mag;
        mag = binNumber;

        // Extracting frequency at bin peak
        float f = trueFrequency[mag];

        NSLog(@"True frequency = %fHz", f);

        float b = roundf(binNumber*(sampleRate/nOver2));

        NSLog(@" Bin frequency = %fHz", b);

    }

Answer 1

Note that the expected phase difference (even for a bin-centered frequency) depends on both the window offset or overlap of the FFT pairs, and the bin number or frequency of the FFT result. e.g. If you offset the windows by very little (1 sample), then the unwrapped phase difference between 2 FFTs will be smaller than with a larger offset. At the same offset, if the frequency is higher, the expected phase difference between the same bin of two FFTs will be greater (or it will wrap more).