iPhone Accelerate Framework FFT vs Matlab FFT

Question

I do not have much math background but part of the project I am working on requires the FFT of a single vector. The matlab function fft(x) works accurately for what I need, but after trying to set up the Accelerate Framework fft functions I get completely inaccurate results. If anyone has more expertise/experience with Accelerate Framework fft I could really use some help trying to figure out what I am doing wrong. I based my fft set-up off an example I found on google, but there were no tutorials or anything that produced different results.

EDIT1: Changed around some stuff based on the answers so far. It seems to be doing calculations but it doesnt output them in any way close to that of matlab

This is the documentation for fft for matlab: http://www.mathworks.com/help/techdoc/ref/fft.html

** NOTE: for example purposes, the x array will be {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} in both examples

Matlab Code:

x = fft(x)

Matlab output:

x =

   1.0e+02 *

  Columns 1 through 4

   1.3600            -0.0800 + 0.4022i  -0.0800 + 0.1931i  -0.0800 + 0.1197i

  Columns 5 through 8

  -0.0800 + 0.0800i  -0.0800 + 0.0535i  -0.0800 + 0.0331i  -0.0800 + 0.0159i

  Columns 9 through 12

  -0.0800            -0.0800 - 0.0159i  -0.0800 - 0.0331i  -0.0800 - 0.0535i

  Columns 13 through 16

  -0.0800 - 0.0800i  -0.0800 - 0.1197i  -0.0800 - 0.1931i  -0.0800 - 0.4022i

Apple Accelerate Framework: http://developer.apple.com/library/mac/#documentation/Accelerate/Reference/vDSPRef/Reference/reference.html#//apple_ref/doc/uid/TP40009464

Objective C code:

int log2n = log2f(16);

FFTSetupD fftSetup = vDSP_create_fftsetupD (log2n, kFFTRadix2);     

DSPDoubleSplitComplex fft_data;
fft_data.realp = (double *)malloc(8 * sizeof(double));
fft_data.imagp = (double *)malloc(8 * sizeof(double));

vDSP_ctoz((COMPLEX *) ffx, 2, &fft_data, 1, nOver2); //split data (1- 16) into odds and evens

vDSP_fft_zrip (fftSetup, &fft_data, 1, log2n, kFFTDirection_Forward); //fft forward

vDSP_fft_zrip (fftSetup, &fft_data, 1, log2n, kFFTDirection_Inverse); //fft inverse

vDSP_ztoc(&fft_data, 2, (COMPLEX *) ffx, 1, nOver2); //combine complex back into real numbers

Objective C output:

ffx now contains:

272.000000
-16.000000
-16.000000
-16.000000
0.000000
0.000000
0.000000
0.000000
0.000000
10.000000
11.000000
12.000000
13.000000
14.000000
15.000000
16.000000

Answer 1

One big problem: C arrays are indexed from 0, unlike MATLAB arrays which are 1-based. So you need to change your loop from

for(int i = 1; i <= 16; i++)

to

for(int i = 0; i < 16; i++)

A second, big problem - you're mixing single precision (float) and double precision (double) routines. Your data is double so you should be using vDSP_ctozD, not vDSP_ctoz, and vDSP_fft_zripD rather than vDSP_fft_zrip, etc.

Another thing to watch out for: different FFT implementations use different definitions of the DFT formula, particularly in regard to scaling factor. It looks like the MATLAB FFT includes a 1/N scaling correction, which most other FFTs do not.

---

Here is a complete working example whose output matches Octave (MATLAB clone):

#include <stdio.h>
#include <stdlib.h>
#include <Accelerate/Accelerate.h>

int main(void)
{
    const int log2n = 4;
    const int n = 1 << log2n;
    const int nOver2 = n / 2;

    FFTSetupD fftSetup = vDSP_create_fftsetupD (log2n, kFFTRadix2);

    double *input;

    DSPDoubleSplitComplex fft_data;

    int i;

    input = malloc(n * sizeof(double));
    fft_data.realp = malloc(nOver2 * sizeof(double));
    fft_data.imagp = malloc(nOver2 * sizeof(double));

    for (i = 0; i < n; ++i)
    {
       input[i] = (double)(i + 1);
    }

    printf("Input\n");

    for (i = 0; i < n; ++i)
    {
        printf("%d: %8g\n", i, input[i]);
    }

    vDSP_ctozD((DSPDoubleComplex *)input, 2, &fft_data, 1, nOver2);

    printf("FFT Input\n");

    for (i = 0; i < nOver2; ++i)
    {
        printf("%d: %8g%8g\n", i, fft_data.realp[i], fft_data.imagp[i]);
    }

    vDSP_fft_zripD (fftSetup, &fft_data, 1, log2n, kFFTDirection_Forward);

    printf("FFT output\n");

    for (i = 0; i < nOver2; ++i)
    {
        printf("%d: %8g%8g\n", i, fft_data.realp[i], fft_data.imagp[i]);
    }

    for (i = 0; i < nOver2; ++i)
    {
        fft_data.realp[i] *= 0.5;
        fft_data.imagp[i] *= 0.5;
    }

    printf("Scaled FFT output\n");

    for (i = 0; i < nOver2; ++i)
    {
        printf("%d: %8g%8g\n", i, fft_data.realp[i], fft_data.imagp[i]);
    }

    printf("Unpacked output\n");

    printf("%d: %8g%8g\n", 0, fft_data.realp[0], 0.0); // DC
    for (i = 1; i < nOver2; ++i)
    {
        printf("%d: %8g%8g\n", i, fft_data.realp[i], fft_data.imagp[i]);
    }
    printf("%d: %8g%8g\n", nOver2, fft_data.imagp[0], 0.0); // Nyquist

    return 0;
}

Output is:

Input
0:        1
1:        2
2:        3
3:        4
4:        5
5:        6
6:        7
7:        8
8:        9
9:       10
10:       11
11:       12
12:       13
13:       14
14:       15
15:       16
FFT Input
0:        1       2
1:        3       4
2:        5       6
3:        7       8
4:        9      10
5:       11      12
6:       13      14
7:       15      16
FFT output
0:      272     -16
1:      -16 80.4374
2:      -16 38.6274
3:      -16 23.9457
4:      -16      16
5:      -16 10.6909
6:      -16 6.62742
7:      -16  3.1826
Scaled FFT output
0:      136      -8
1:       -8 40.2187
2:       -8 19.3137
3:       -8 11.9728
4:       -8       8
5:       -8 5.34543
6:       -8 3.31371
7:       -8  1.5913
Unpacked output
0:      136       0
1:       -8 40.2187
2:       -8 19.3137
3:       -8 11.9728
4:       -8       8
5:       -8 5.34543
6:       -8 3.31371
7:       -8  1.5913
8:       -8       0

Comparing with Octave we get:

octave-3.4.0:15> x = [ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 ]
x =

    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16

octave-3.4.0:16> fft(x)
ans =

 Columns 1 through 7:

   136.0000 +   0.0000i    -8.0000 +  40.2187i    -8.0000 +  19.3137i    -8.0000 +  11.9728i    -8.0000 +   8.0000i    -8.0000 +   5.3454i    -8.0000 +   3.3137i

 Columns 8 through 14:

    -8.0000 +   1.5913i    -8.0000 +   0.0000i    -8.0000 -   1.5913i    -8.0000 -   3.3137i    -8.0000 -   5.3454i    -8.0000 -   8.0000i    -8.0000 -  11.9728i

 Columns 15 and 16:

    -8.0000 -  19.3137i    -8.0000 -  40.2187i

octave-3.4.0:17>

Note that the outputs from 9 to 16 are just a complex conjugate mirror image or the bottom 8 terms, as is the expected case with a real-input FFT.

Note also that we needed to scale the vDSP FFT by a factor of 2 - this is due to the fact that it's a real-to-complex FFT, which is based on an N/2 point complex-to-complex FFT, hence the outputs are scaled by N/2, whereas a normal FFT would be scaled by N.

Answer 2

I think it could also be an array packing issue. I have just been looking at their sample code, and I see they keep calling conversion routines like

vDSP_ctoz

Here is the code sample link from apple: http://developer.apple.com/library/ios/#documentation/Performance/Conceptual/vDSP_Programming_Guide/SampleCode/SampleCode.html#//apple_ref/doc/uid/TP40005147-CH205-CIAEJIGF

I dont think it is the full answer yet, but I also agree with Paul R.

By the way, just curiously if you go to Wolfram Alpha, they give a completely different answer for FFT{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16}