我在使用苹果的Accelerate Framework时,在C++中获得了一个与Matlab等效的Hilbert变换存在问题。我已经成功地让vDSP的FFT算法工作,并且在Paul R的帖子的帮助下,已经成功地得到了与Matlab相同的结果。
我阅读了这个Jordan的stackoverflow问题,并阅读了Matlab算法(在“算法”子标题下)。为了总结这个算法,分为三个阶段:
1. 对输入进行正向FFT。 2. 零反射频率和DC到Nyquist之间的双倍频率。 3. 对修改后的正向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的“第三阶段”结果(或至少是虚部),但我不知道该怎么做,我已经尝试了我能想到的一切,但没有成功。我有一种感觉,这是因为在苹果加速版本中没有将反射频率清零,因为它们没有提供(由于是从DC和Nyquist之间的频率计算出来的)-所以我认为FFT只是将双频的共轭作为反射值,这是错误的。如果有人能帮我解决这个问题,我将非常感激。
我阅读了这个Jordan的stackoverflow问题,并阅读了Matlab算法(在“算法”子标题下)。为了总结这个算法,分为三个阶段:
1. 对输入进行正向FFT。 2. 零反射频率和DC到Nyquist之间的双倍频率。 3. 对修改后的正向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}`。
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
第二阶段:
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++示例(使用苹果的加速框架)
阶段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
第二阶段:
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的“第三阶段”结果(或至少是虚部),但我不知道该怎么做,我已经尝试了我能想到的一切,但没有成功。我有一种感觉,这是因为在苹果加速版本中没有将反射频率清零,因为它们没有提供(由于是从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++(使用苹果的Accelerate Framework)结果的。前向fft的fft(..)
方法执行vDSP fft例程,然后按比例缩放0.5并解压缩(根据Paul R的帖子)。当执行反向fft时,数据被打包,乘以2.0缩放,使用vDSP进行fft,最后乘以1 /(2 * N)缩放。