使用 Apple Accelerate 框架的希尔伯特变换（分析信号）？

Posted 2023-02-25

【中文标题】使用 Apple Accelerate 框架的希尔伯特变换（分析信号）？

【英文标题】：Hilbert Transform (Analytical Signal) using Apple's Accelerate Framework?

【发布时间】：2014-02-19 22:02:01

【问题描述】：

我在使用 Apple 的 Accelerate Framework 在 C++ 中获得 Matlab 等效希尔伯特变换时遇到问题。我已经能够让 vDSP 的 FFT 算法正常工作，并且在 Paul R's post 的帮助下，我设法获得了与 Matlab 相同的结果。

我都读过：这个*** question by Jordan 和Matlab algorithm (under the 'Algorithms' sub-heading)。将算法总结为三个阶段：

对输入进行 FFT。 直流和奈奎斯特之间的零反射频率和双倍频率。 对修改后的正向 FFT 输出进行逆 FFT。

下面是每个阶段的 Matlab 和 C++ 的输出。示例使用以下数组：

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

---

Matlab 示例

---

第一阶段：

  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

第 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

第三阶段：

   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++ 示例（使用 Apple 的 Accelerate 框架）

---

第一阶段：

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

第 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

第三阶段：

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

很明显，最终结果并不相同。我希望能够计算 Matlab 'Stage 3' 结果（或至少是虚构部分），但我不确定如何去做，我已经尝试了我能想到的一切，但没有成功。我有一种感觉，这是因为我没有将 Apple Accelerate 版本中的反射频率归零，因为它们没有提供（由于是根据 DC 和 Nyquist 之间的频率计算的） - 所以我认为 FFT 只是采用共轭将双倍频率作为反射值，这是错误的。如果有人能帮助我解决这个问题，我将不胜感激。

---

代码

---

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));

上面的代码是我用来获得“第 3 阶段”C++（使用 Apple 的 Accelerate Framework）结果的代码。前向 fft 的 fft(..) 方法执行 vDSP fft 例程，后跟 0.5 的比例，然后解包（根据 Paul R 的帖子）。当执行逆 fft 时，将数据打包，缩放 2.0，使用 vDSP 进行 fft，最后缩放 1/(2*N)。

【问题讨论】：

如果您发布您的代码，对您的帮助会更容易一些（只需此处转换的部分就足够了）。

代码已附加到问题中。 :)

tldr 但如果您是 Hilbert 变换信号，您可能希望使用重叠相加。即一次通过信号 1024 个样本，FFT 窗口为 4096。FFT 之前的包络（升余弦或更好的东西）。摆弄你的 FFT 箱，逆 FFT 并合成你的重叠帧。您可以轻松地调整 Stefan's code 并在完成该操作后使用 Accelerate 进行优化。

【参考方案1】：

所以主要问题（据我所知，因为您的代码示例不包括对 vDSP 的实际调用）是您尝试在第三步中使用实数到复数 FFT 例程，其中基本上是复数到复数的逆 FFT（您的 Matlab 结果具有非零虚部这一事实证明了这一点）。

这是一个使用 vDSP 的简单 C 实现，它与您预期的 Matlab 输出相匹配（我使用了更现代的 vDSP_DFT 例程，通常应该优先于旧的 fft 例程，但除此之外这与您正在做的非常相似，并且说明需要实数到复数的正向变换，但需要复数到复数的逆变换）：

#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;

请注意，如果您已经构建了 DFT 设置并分配了存储空间，则计算非常简单； “真正的工作”只是：

  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);

---

样本输出：

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

【讨论】：

对于任何利用这个出色答案的人来说只是一个小提示 - 实数和图像数组就足够了，不需要额外的偶数、奇数和希尔伯特数组。只需用偶数和奇数内容初始化 real 和 imag 的前半部分，然后将后半部分设为零。