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
");
for (i = 0; i < n; ++i)
{
printf("%d: %8g
", i, input[i]);
}
vDSP_ctozD((DSPDoubleComplex *)input, 2, &fft_data, 1, nOver2);
printf("FFT Input
");
for (i = 0; i < nOver2; ++i)
{
printf("%d: %8g%8g
", i, fft_data.realp[i], fft_data.imagp[i]);
}
vDSP_fft_zripD (fftSetup, &fft_data, 1, log2n, kFFTDirection_Forward);
printf("FFT output
");
for (i = 0; i < nOver2; ++i)
{
printf("%d: %8g%8g
", 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
");
for (i = 0; i < nOver2; ++i)
{
printf("%d: %8g%8g
", i, fft_data.realp[i], fft_data.imagp[i]);
}
printf("Unpacked output
");
printf("%d: %8g%8g
", 0, fft_data.realp[0], 0.0); // DC
for (i = 1; i < nOver2; ++i)
{
printf("%d: %8g%8g
", i, fft_data.realp[i], fft_data.imagp[i]);
}
printf("%d: %8g%8g
", 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.
