title: Hilbert Transform (Analytical Signal) using Apple's Accelerate Framework?
body: <p>I am having issues with getting a <code>Matlab</code> equivalent Hilbert transform in <code>C++</code> with using Apple's <code>Accelerate Framework</code>. I have been able to get vDSP's FFT algorithm working and, with the help of <a href="https://stackoverflow.com/a/10823152/1802891">Paul R's post</a>, have managed to get the same outcome as Matlab.</p>
  
  <p>I have read both: this <a href="https://stackoverflow.com/q/11942855/1802891">stackoverflow question by Jordan</a> and have read the <a href="http://www.mathworks.co.uk/help/signal/ref/hilbert.html" rel="nofollow noreferrer">Matlab algorithm (under the 'Algorithms' sub-heading)</a>. To sum the algorithm up in 3 stages:</p>
  
  <ol>
  <li>Take forward FFT of input.</li>
  <li>Zero reflection frequencies and double frequencies between DC and Nyquist.</li>
  <li>Take inverse FFT of the modified forward FFT output.</li>
  </ol>
  
  <p>Below are the outputs of both Matlab and C++ for each stage. The examples use the following arrays:</p>
  
  <ul>
  <li>Matlab: <code>m = [1 2 3 4 5 6 7 8];</code></li>
  <li>C++: <code>float m[] = {1,2,3,4,5,6,7,8};</code></li>
  </ul>
  
  <hr>
  
  <h2>Matlab Example</h2>
  
  <hr>
  
  <p><em>Stage 1:</em></p>
  
  <pre><code>  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
  </code></pre>
  
  <p><em>Stage 2:</em></p>
  
  <pre><code>  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
  </code></pre>
  
  <p><em>Stage 3:</em></p>
  
  <pre><code>   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
  </code></pre>
  
  <hr>
  
  <h2>C++ Example (with Apple's Accelerate Framework)</h2>
  
  <hr>
  
  <p><em>Stage 1:</em></p>
  
  <pre><code>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
  </code></pre>
  
  <p><em>Stage 2:</em></p>
  
  <pre><code>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
  </code></pre>
  
  <p><em>Stage 3:</em></p>
  
  <pre><code>Real: -2.0000, Imag: -1.0000
  Real:  2.0000, Imag: 3.0000
  Real:  6.0000, Imag: 7.0000
  Real: 10.0000, Imag: 11.0000
  </code></pre>
  
  <p>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.</p>
  
  <hr>
  
  <h2>Code</h2>
  
  <hr>
  
  <pre><code>void hilbert(std::vector&lt;float&gt; &amp;data, std::vector&lt;float&gt; &amp;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 &lt; dataSize_2; i++)
          workspace.at(i) = 2.0f;
  
      // Multiply forward FFT output by workspace vector.
      for (i = 0; i &lt; 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 &lt; hilbertResult.size()/2; i++)
          printf("Real: %.4f, Imag: %.4f\n", hilbertResult.at(i*2), hilbertResult.at(i*2+1));
  }
  </code></pre>
  
  <p>The code above is what I used to get the 'Stage 3' C++ (with Apple's Accelerate Framework) result. The <code>fft(..)</code> 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).</p>