Преобразование Гильберта (аналитический сигнал) с использованием Apple Accelerate Framework?

У меня проблемы с получением Matlab эквивалентное преобразование Гильберта в C++ с использованием Apple Accelerate Framework, Я смог заставить работать алгоритм БПФ vDSP и, с помощью Пост Пола Р., удалось получить тот же результат, что и Matlab.

Я прочитал оба: это StackOverflow вопрос Джордан и прочитал Алгоритм Matlab (в подразделе «Алгоритмы»). Подвести алгоритм в 3 этапа:

1. Возьмите вперед БПФ ввода.
2. Нулевые частоты отражения и двойные частоты между DC и Найквистом.
3. Возьмите обратное БПФ измененного выхода прямого БПФ.

Ниже приведены результаты как Matlab, так и C ++ для каждого этапа. В примерах используются следующие массивы:

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

---

Пример Matlab

---

Этап 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

Этап 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

Этап 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 ++ (с Apple Accelerate Framework)

---

Этап 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

Этап 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

Этап 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

Совершенно ясно, что конечные результаты не совпадают. Я надеюсь, что смогу вычислить результаты Matlab ‘Stage 3’ (или, по крайней мере, мнимые части), но я не уверен, как это сделать, я перепробовал все, что мог придумать, но безуспешно. У меня такое ощущение, что это потому, что я не обнуляю частоты отражения в версии Apple Accelerate, поскольку они не предоставляются (из-за того, что они рассчитываются по частотам между DC и Nyquist) — поэтому я думаю, что БПФ просто принимает сопряженное удвоенных частот в качестве значений отражения, что неправильно. Если бы кто-нибудь мог помочь мне преодолеть эту проблему, я был бы очень признателен.

---

Код

---

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

Приведенный выше код — это то, что я использовал, чтобы получить результат ‘Stage 3’ C ++ (с Apple Accelerate Framework). fft(..) Метод forward fft выполняет подпрограмму vDSP fft, за которой следует масштаб 0,5, а затем распаковывается (согласно записи Пола Р.). Когда выполняется обратное fft, данные упаковываются, масштабируются на 2.0, fft’d с использованием vDSP и, наконец, масштабируются на 1 / (2 * N).

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

  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