使用 Apple FFT 和加速框架

下面是一个真实的例子: c + + 的一个片段，它使用 Speed 的 vDSP fft 例程对远程 IO 音频单元的输入进行自相关处理。使用这个框架是相当复杂的，但是文档并不是 也是不好。

OSStatus DSPCore::initialize (double _sampleRate, uint16_t _bufferSize) {
sampleRate = _sampleRate;
bufferSize = _bufferSize;
peakIndex = 0;
frequency = 0.f;
uint32_t maxFrames = getMaxFramesPerSlice();
displayData = (float*)malloc(maxFrames*sizeof(float));
bzero(displayData, maxFrames*sizeof(float));
log2n = log2f(maxFrames);
n = 1 << log2n;
assert(n == maxFrames);
nOver2 = maxFrames/2;
A.realp = (float*)malloc(nOver2 * sizeof(float));
A.imagp = (float*)malloc(nOver2 * sizeof(float));
FFTSetup fftSetup = vDSP_create_fftsetup(log2n, FFT_RADIX2);


return noErr;
}


void DSPCore::Render(uint32_t numFrames, AudioBufferList *ioData) {


bufferSize = numFrames;
float ln = log2f(numFrames);


//vDSP autocorrelation


//convert real input to even-odd
vDSP_ctoz((COMPLEX*)ioData->mBuffers[0].mData, 2, &A, 1, numFrames/2);
memset(ioData->mBuffers[0].mData, 0, ioData->mBuffers[0].mDataByteSize);
//fft
vDSP_fft_zrip(fftSetup, &A, 1, ln, FFT_FORWARD);


// Absolute square (equivalent to mag^2)
vDSP_zvmags(&A, 1, A.realp, 1, numFrames/2);
bzero(A.imagp, (numFrames/2) * sizeof(float));


// Inverse FFT
vDSP_fft_zrip(fftSetup, &A, 1, ln, FFT_INVERSE);


//convert complex split to real
vDSP_ztoc(&A, 1, (COMPLEX*)displayData, 2, numFrames/2);


// Normalize
float scale = 1.f/displayData[0];
vDSP_vsmul(displayData, 1, &scale, displayData, 1, numFrames);


// Naive peak-pick: find the first local maximum
peakIndex = 0;
for (size_t ii=1; ii < numFrames-1; ++ii) {
if ((displayData[ii] > displayData[ii-1]) && (displayData[ii] > displayData[ii+1])) {
peakIndex = ii;
break;
}
}


// Calculate frequency
frequency = sampleRate / peakIndex + quadInterpolate(&displayData[peakIndex-1]);


bufferSize = numFrames;


for (int ii=0; ii<ioData->mNumberBuffers; ++ii) {
bzero(ioData->mBuffers[ii].mData, ioData->mBuffers[ii].mDataByteSize);
}
}

我刚得到了一个 iPhone 项目的 FFT 代码:

• 创建一个新项目
• 删除除 main.m 和 xxx _ info. plist 之外的所有文件
• 进入项目设置并搜索 pch 并阻止它试图加载一个。Pch (因为我们刚刚删除了它)
• 复制粘贴代码示例到 main.m 中的任何内容上
• 删除 # include 的 Carbon。 Carbon 是 OSX 的。
• 删除所有框架，并添加加速框架

您可能还需要从 info.plist 中删除一个条目，该条目告诉项目加载 xib，但是我敢肯定您不需要为此费心。

注意: 程序输出到控制台，结果出来为0.000，这不是一个错误-- 它只是非常非常快

这段代码非常晦涩难懂; 它被大量地注释了，但是注释实际上并没有使生活变得更容易。

基本上它的核心是:

vDSP_fft_zrip(setupReal, &A, stride, log2n, FFT_FORWARD);
vDSP_fft_zrip(setupReal, &A, stride, log2n, FFT_INVERSE);

对 n 个实际的浮点数进行 FFT，然后返回到我们开始的地方。 Ip 代表就位，也就是说 & A 被覆盖了 这就是所有这些特殊的打包废话的原因——这样我们就可以将返回值压缩到与 send 值相同的空间中。

给出一些观点(比如: 我们为什么要首先使用这个函数?)，让我们说，我们想要执行音高检测麦克风输入，我们已经设置，以便一些回调被触发每次麦克风得到1024浮动。假设麦克风采样频率为44.1 kHz，那么就是44帧/秒。

因此，我们的时间窗口是1024个样本的时间长度，即1/44秒。

因此，我们将用麦克风中的1024个浮点数来包装 A，设置 log2n = 10(2 ^ 10 = 1024) ，预先计算一些线轴(setupReal) ，然后:

vDSP_fft_zrip(setupReal, &A, stride, log2n, FFT_FORWARD);

现在 A 包含 n/2个复数，它们代表 n/2个频率箱:

• 垃圾箱[1]。理想频率 = 44赫兹-即我们可以可靠地检测到的最低频率是在该窗口内的一个完整的波，即44赫兹的波。

• Bin [2] . idealFreq = 2 * 44Hz

• 等等。

• 垃圾箱[512]。IdealFreq = 512 * 44Hz ——我们能探测到的最高频率(称为奈奎斯特频率)是每对点代表一个波，即窗口内的512个完整波，即512 * 44Hz，或: n/2 * bin [1]。理想频率

• 实际上还有一个额外的 Bin，Bin [0] ，它通常被称为“ DC 偏移量”。碰巧 Bin [0]和 Bin [ n/2]总是有复杂的组件0，所以 A [0]。Realp 用于存储 Bin [0]和 A [0]。Image 用于存储 Bin [ n/2]

每个复数的大小就是围绕这个频率振动的能量。

所以，正如你所看到的，它不会是一个非常好的音高检测器，因为它没有足够细的粒度。有一个狡猾的技巧 利用帧间相位变化从 FFT 帧中提取精确频率获得精确的频率给定的垃圾桶。

好了，现在说说密码:

注意 vDSP _ fft _ zlip 中的‘ ip’，= ‘ in place’即输出覆盖 A (‘ r’表示它接受真正的输入)

查看 vDSP _ fft _ zlip 上的文档,

实际数据存储在拆分复合体中 形式，与奇数雷亚尔存储在 分裂复合体的虚边 形式，甚至真实存储在 真实的一面。

这可能是最难理解的事情。我们在整个过程中使用相同的容器(& A)。所以一开始我们要用 n 个实数来填充。FFT 之后是 n/2复数。然后我们把它转换成逆变换，希望得到原来的 n 个实数。

现在是 A 的结构它的复数值设置。所以 vDSP 需要标准化如何将真实数字包装进去。

首先我们生成 n 个实数: 1,2，... ，n

for (i = 0; i < n; i++)
originalReal[i] = (float) (i + 1);

接下来我们把它们打包成 n/2复数形式的 A # :

// 1. masquerades n real #s as n/2 complex #s = {1+2i, 3+4i, ...}
// 2. splits to
//   A.realP = {1,3,...} (n/2 elts)
//   A.compP = {2,4,...} (n/2 elts)
//
vDSP_ctoz(
(COMPLEX *) originalReal,
2,                            // stride 2, as each complex # is 2 floats
&A,
1,                            // stride 1 in A.realP & .compP
nOver2);                      // n/2 elts

您确实需要了解 A 是如何分配的，也许可以在文档中查找 COMPLEX _ SPLIT。

A.realp = (float *) malloc(nOver2 * sizeof(float));
A.imagp = (float *) malloc(nOver2 * sizeof(float));

接下来我们做一个预先计算。

数学天才的快速 DSP 课程: 傅立叶理论需要很长时间才能让你理解(我已经断断续续看了好几年了)

一个坐标系是:

z = exp(i.theta) = cos(theta) + i.sin(theta)

也就是复平面上单位圆上的一点。

当你把复数相乘时，角相加。所以 z ^ k 会一直绕着单位圆跳动，z ^ k 可以在角 k.θ 处找到

• 选择 z1 = 0 + 1 i，也就是离实轴四分之一圈，注意 z1 ^ 2 z1 ^ 3 z1 ^ 4每个都给出另一个四分之一圈，这样 z1 ^ 4 = 1

• 选择 z2 = -1，即半圈。也是 z2 ^ 4 = 1，但是 z2在这一点上已经完成了2个循环(z2 ^ 2也是 = 1)。所以你可以认为 z1是基频，而 z2是第一谐波

• 类似地，z3 = “四分之三的旋转”点，也就是-i 刚好完成3个周期，但实际上每次向前3/4和向后1/4是一样的

也就是说，z3只是 z1，但是在相反的方向上，它被称为混叠

Z2是最高的有意义的频率，因为我们选择了4个样本持有一个完整的波。

• Z0 = 1 + 0 i，z0 ^ (anything) = 1，这是直流偏移量

你可以用 z0z1和 z2的线性组合表示任意四点信号 也就是说，你把它投影到这些基本向量上

但我听到你在问“将信号投射到一个西索上意味着什么?”

你们可以这样想: 指针绕着纵坐标旋转，所以在样本 k 上，指针指向 k.theta 方向，长度是信号[ k ]。如果信号的频率与坐标系的频率完全匹配，那么就会在某个方向上形成凸起的形状。所以如果你把所有的贡献加起来，你会得到一个强合成矢量。 如果频率几乎匹配，则凸起将更小，并且在圆周上缓慢移动。 对于一个与频率不匹配的信号，贡献会相互抵消

Http://complextoreal.com/tutorials/tutorial-4-fourier-analysis-made-easy-part-1/将帮助你获得直观的理解。

// let's say log2n = 8, so n=2^8=256 samples, or 'harmonics' or 'terms'
// if we pre-calculate the 256th roots of unity (of which there are 256)
// that will save us time later.
//
// Note that this call creates an array which will need to be released
// later to avoid leaking
setupReal = vDSP_create_fftsetup(log2n, FFT_RADIX2);

vDSP_fft_zrip(setupReal, &A, stride, log2n, FFT_FORWARD);

vDSP_fft_zrip(setupReal, &A, stride, log2n, FFT_INVERSE);

// Need to see the documentation for this one...
// in order to optimise, different routines return values
// that need to be scaled by different amounts in order to
// be correct as per the math
// In this case...
scale = (float) 1.0 / (2 * n);


vDSP_vsmul(A.realp, 1, &scale, A.realp, 1, nOver2);
vDSP_vsmul(A.imagp, 1, &scale, A.imagp, 1, nOver2);

虽然我会说苹果的 FFT 框架是快速的... 你需要知道一个 FFT 是如何工作的，以便获得准确的基音检测(即计算每个连续的 FFT 的相位差，以便找到准确的基音，而不是最占主导地位的垃圾桶)。

我不知道这是否有帮助，但我上传了我的音调检测器对象从我的调谐器应用程序( musicianskit.com/developer.php )。还有一个示例 xCode4项目可供下载(这样您就可以看到实现是如何工作的)。

我正在上传一个 FFT 实现示例——所以请继续关注，一旦发生这种情况，我将更新它。

编程愉快！

下面是另一个真实世界的例子: Https://github.com/krafter/detectingaudiofrequency