Hilbert Transform (Analytical Signal) using Apple&

I am having issues with getting a Matlab equivalent Hilbert transform in C++ with using Apple's Accelerate Framework. I have been able to get vDSP's FFT algorithm working and, with the help of Paul R's post, have managed to get the same outcome as Matlab.

I have read both: this stackoverflow question by Jordan and have read the Matlab algorithm (under the 'Algorithms' sub-heading). To sum the algorithm up in 3 stages:

Take forward FFT of input.
Zero reflection frequencies and double frequencies between DC and Nyquist.
Take inverse FFT of the modified forward FFT output.

Below are the outputs of both Matlab and C++ for each stage. The examples use the following arrays:

Matlab: m = [1 2 3 4 5 6 7 8];
C++: float m[] = {1,2,3,4,5,6,7,8};

Matlab Example

Stage 1:

  36.0000 + 0.0000i
  -4.0000 + 9.6569i
  -4.0000 + 4.0000i
  -4.0000 + 1.6569i
  -4.0000 + 0.0000i
  -4.0000 - 1.6569i
  -4.0000 - 4.0000i
  -4.0000 - 9.6569i

Stage 2:

  36.0000 + 0.0000i
  -8.0000 + 19.3137i
  -8.0000 + 8.0000i
  -8.0000 + 3.3137i
  -4.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i

Stage 3:

   1.0000 + 3.8284i
   2.0000 - 1.0000i
   3.0000 - 1.0000i
   4.0000 - 1.8284i
   5.0000 - 1.8284i
   6.0000 - 1.0000i
   7.0000 - 1.0000i
   8.0000 + 3.8284i

C++ Example (with Apple's Accelerate Framework)

Stage 1:

Real: 36.0000, Imag: 0.0000
Real: -4.0000, Imag: 9.6569
Real: -4.0000, Imag: 4.0000
Real: -4.0000, Imag: 1.6569
Real: -4.0000, Imag: 0.0000

Stage 2:

Real: 36.0000, Imag: 0.0000
Real: -8.0000, Imag: 19.3137
Real: -8.0000, Imag: 8.0000
Real: -8.0000, Imag: 3.3137
Real: -4.0000, Imag: 0.0000

Stage 3:

Real: -2.0000, Imag: -1.0000
Real:  2.0000, Imag: 3.0000
Real:  6.0000, Imag: 7.0000
Real: 10.0000, Imag: 11.0000

It's quite clear that the end results are not the same. I am looking to be able to compute the Matlab 'Stage 3' results (or at least the imaginary parts) but I am unsure how to go about it, I've tried everything I can think of with no success. I have a feeling that it's because I'm not zeroing out the reflection frequencies in the Apple Accelerate version as they are not provided (due to being calculated from the frequencies between DC and Nyquist) - so I think the FFT is just taking the conjugate of the doubled frequencies as the reflection values, which is wrong. If anyone could help me overcome this issue I would greatly appreciate it.

Code

void hilbert(std::vector<float> &data, std::vector<float> &hilbertResult){

    // Set variable.
    dataSize_2 = data.size() * 0.5;

    // Clear and resize vectors.
    workspace.clear();
    hilbertResult.clear();

    workspace.resize(data.size()/2+1, 0.0f);
    hilbertResult.resize(data.size(), 0.0f);

    // Copy data into the hilbertResult vector.
    std::copy(data.begin(), data.end(), hilbertResult.begin());

    // Perform forward FFT.
    fft(hilbertResult, hilbertResult.size(), FFT_FORWARD);

    // Fill workspace with 1s and 2s (to double frequencies between DC and Nyquist).
    workspace.at(0) = workspace.at(dataSize_2) = 1.0f;

    for (i = 1; i < dataSize_2; i++)
        workspace.at(i) = 2.0f;

    // Multiply forward FFT output by workspace vector.
    for (i = 0; i < workspace.size(); i++) {
        hilbertResult.at(i*2)   *= workspace.at(i);
        hilbertResult.at(i*2+1) *= workspace.at(i);
    }

    // Perform inverse FFT.
    fft(hilbertResult, hilbertResult.size(), FFT_INVERSE);

    for (i = 0; i < hilbertResult.size()/2; i++)
        printf("Real: %.4f, Imag: %.4f\n", hilbertResult.at(i*2), hilbertResult.at(i*2+1));
}

The code above is what I used to get the 'Stage 3' C++ (with Apple's Accelerate Framework) result. The fft(..) method for the forward fft performs the vDSP fft routine followed by a scale of 0.5 and then unpacked (as per Paul R's post). When the inverse fft is performed, the data is packed, scaled by 2.0, fft'd using vDSP and finally scaled by 1/(2*N).

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