So the main problem (as far as I can tell, as your code sample doesn't include the actual calls into vDSP) is that you’re attempting to use the real-to-complex FFT routines for step three, which is fundamentally a complex-to-complex inverse FFT (as evidenced by the fact that your Matlab results have non-zero imaginary parts).
Here’s a simple C implementation using vDSP that matches your expected Matlab output (I used the more modern vDSP_DFT routines, which should generally be preferred to the older fft routines, but otherwise this is very similar to what you’re doing, and illustrates the need for a real-to-complex forward transform, but complex-to-complex inverse transform):
#include <Accelerate/Accelerate.h>
#include <stdio.h>
int main(int argc, char *argv[]) {
const vDSP_Length n = 8;
vDSP_DFT_SetupD forward = vDSP_DFT_zrop_CreateSetupD(NULL, n, vDSP_DFT_FORWARD);
vDSP_DFT_SetupD inverse = vDSP_DFT_zop_CreateSetupD(forward, n, vDSP_DFT_INVERSE);
// Look like a typo? The real-to-complex DFT takes its input separated into
// the even- and odd-indexed elements. Since the real signal is [ 1, 2, 3, ... ],
// signal[0] is 1, signal[2] is 3, and so on for the even indices.
double even[n/2] = { 1, 3, 5, 7 };
double odd[n/2] = { 2, 4, 6, 8 };
double real[n] = { 0 };
double imag[n] = { 0 };
vDSP_DFT_ExecuteD(forward, even, odd, real, imag);
// At this point, we have the forward real-to-complex DFT, which agrees with
// MATLAB up to a factor of two. Since we want to double all but DC and NY
// as part of the Hilbert transform anyway, I'm not going to bother to
// unscale the rest of the frequencies -- they're already the values that
// we really want. So we just need to move NY into the "right place",
// and scale DC and NY by 0.5. The reflection frequencies are already
// zeroed out because the real-to-complex DFT only writes to the first n/2
// elements of real and imag.
real[0] *= 0.5; real[n/2] = 0.5*imag[0]; imag[0] = 0.0;
printf("Stage 2:\n");
for (int i=0; i<n; ++i) printf("%f%+fi\n", real[i], imag[i]);
double hilbert[2*n];
double *hilbertreal = &hilbert[0];
double *hilbertimag = &hilbert[n];
vDSP_DFT_ExecuteD(inverse, real, imag, hilbertreal, hilbertimag);
// Now we have the completed hilbert transform up to a scale factor of n.
// We can unscale using vDSP_vsmulD.
double scale = 1.0/n; vDSP_vsmulD(hilbert, 1, &scale, hilbert, 1, 2*n);
printf("Stage 3:\n");
for (int i=0; i<n; ++i) printf("%f%+fi\n", hilbertreal[i], hilbertimag[i]);
vDSP_DFT_DestroySetupD(inverse);
vDSP_DFT_DestroySetupD(forward);
return 0;
}
Note that if you already have your DFT setups built and your storage allocated, the computation is extremely straightforward; the “real work” is just:
vDSP_DFT_ExecuteD(forward, even, odd, real, imag);
real[0] *= 0.5; real[n/2] = 0.5*imag[0]; imag[0] = 0.0;
vDSP_DFT_ExecuteD(inverse, real, imag, hilbertreal, hilbertimag);
double scale = 1.0/n; vDSP_vsmulD(hilbert, 1, &scale, hilbert, 1, 2*n);
Sample output:
Stage 2:
36.000000+0.000000i
-8.000000+19.313708i
-8.000000+8.000000i
-8.000000+3.313708i
-4.000000+0.000000i
0.000000+0.000000i
0.000000+0.000000i
0.000000+0.000000i
Stage 3:
1.000000+3.828427i
2.000000-1.000000i
3.000000-1.000000i
4.000000-1.828427i
5.000000-1.828427i
6.000000-1.000000i
7.000000-1.000000i
8.000000+3.828427i