iPhone FFT with Accelerate framework vDSP

Question

I'm having difficulty implementing an FFT using vDSP. I understand the theory but am looking for a specific code example please.

I have data from a wav file as below:

Question 1. How do I put the audio data into the FFT?

Question 2. How do I get the output data out of the FFT?

Question 3. The ultimate goal is to check for low frequency sounds. How would I do this?

-(OSStatus)open:(CFURLRef)inputURL{
OSStatus result = -1;

result = AudioFileOpenURL (inputURL, kAudioFileReadPermission, 0, &mAudioFile);
if (result == noErr) {
    //get  format info
    UInt32 size = sizeof(mASBD);

    result = AudioFileGetProperty(mAudioFile, kAudioFilePropertyDataFormat, &size, &mASBD);

    UInt32 dataSize = sizeof packetCount;
    result = AudioFileGetProperty(mAudioFile, kAudioFilePropertyAudioDataPacketCount, &dataSize, &packetCount);
    NSLog([NSString stringWithFormat:@"File Opened, packet Count: %d", packetCount]);

    UInt32 packetsRead = packetCount;
    UInt32 numBytesRead = -1;
    if (packetCount > 0) { 
        //allocate  buffer
        audioData = (SInt16*)malloc( 2 *packetCount);
        //read the packets
        result = AudioFileReadPackets (mAudioFile, false, &numBytesRead, NULL, 0, &packetsRead,  audioData); 
        NSLog([NSString stringWithFormat:@"Read %d  bytes,  %d packets", numBytesRead, packetsRead]);
    }
}
return result;
}

FFT code below:

log2n = N;
n = 1 << log2n;
stride = 1;
nOver2 = n / 2;

printf("1D real FFT of length log2 ( %d ) = %d\n\n", n, log2n);

/* Allocate memory for the input operands and check its availability,
 * use the vector version to get 16-byte alignment. */

A.realp = (float *) malloc(nOver2 * sizeof(float));
A.imagp = (float *) malloc(nOver2 * sizeof(float));
originalReal = (float *) malloc(n * sizeof(float));
obtainedReal = (float *) malloc(n * sizeof(float));

if (originalReal == NULL || A.realp == NULL || A.imagp == NULL) {
printf("\nmalloc failed to allocate memory for  the real FFT"
"section of the sample.\n");
exit(0);
}

/* Generate an input signal in the real domain. */
for (i = 0; i < n; i++)

    originalReal[i] = (float) (i + 1);

/* Look at the real signal as an interleaved complex vector  by
 * casting it.  Then call the transformation function vDSP_ctoz to
 * get a split complex vector, which for a real signal, divides into
 * an even-odd configuration. */

vDSP_ctoz((COMPLEX *) originalReal, 2, &A, 1, nOver2);

/* Set up the required memory for the FFT routines and check  its
 * availability. */

setupReal = vDSP_create_fftsetup(log2n, FFT_RADIX2);

if (setupReal == NULL) {

printf("\nFFT_Setup failed to allocate enough memory  for"
"the real FFT.\n");

exit(0);
}

/* Carry out a Forward and Inverse FFT transform. */
vDSP_fft_zrip(setupReal, &A, stride, log2n, FFT_FORWARD);
vDSP_fft_zrip(setupReal, &A, stride, log2n, FFT_INVERSE);

/* Verify correctness of the results, but first scale it by  2n. */
scale = (float) 1.0 / (2 * n);
vDSP_vsmul(A.realp, 1, &scale, A.realp, 1, nOver2);
vDSP_vsmul(A.imagp, 1, &scale, A.imagp, 1, nOver2);

/* The output signal is now in a split real form.  Use the  function
 * vDSP_ztoc to get a split real vector. */
vDSP_ztoc(&A, 1, (COMPLEX *) obtainedReal, 2, nOver2);

/* Check for accuracy by looking at the inverse transform  results. */
Compare(originalReal, obtainedReal, n);

Thanks

Answer 1

One thing you need to be careful to is the DC component of the calculated FFT. I compared my results with the fftw library FFT and the imaginary part of the transform calculated with the vDSP library always had a different value at index 0 (which means 0 frequency, so DC). Another measure I applied was to divide both real and imaginary parts by a factor of 2. I guess this is due to the algorithm used in the function. Also, both these problems occurred in the FFT process but not in the IFFT process.

I used vDSP_fft_zrip.

Answer 2

1. You put your audio sample data into the real part of the input, and zero the imaginary part.

2. If you are just interested in the magnitude of each bin in the frequency domain then you calculate sqrt(re*re + im*im) for each output bin. If you're only interested in relative magnitude then you can drop the sqrt and just calculate the squared magnitude, (re*re + im*im).

3. You would look at the magnitudes of the bin or bins (see (2)) that correspond to your frequency or frequencies of interest. If your sample rate is Fs, and your FFT size is N, then the corresponding frequency for output bin i is given by f = i * Fs / N. Conversely if you are interested in a specific frequency f then the bin of interest, i, is given by i = N * f / Fs.

Additional note: you will need to apply a suitable window function (e.g. Hann aka Hanning) to your FFT input data, prior to calculating the FFT itself.