使用 Apple FFT 和 Accelerate 框架

Question

有没有人将Apple FFT 用于 iPhone 应用程序，或者知道我在哪里可以找到关于如何使用它的示例应用程序？我知道 Apple 发布了一些示例代码，但我不确定如何将其实施到实际项目中。

Answer 1

我刚刚获得了适用于 iPhone 项目的 FFT 代码：

• 创建一个新项目
• 删除除main.m和xxx_info.plist之外的所有文件
• 进入项目设置并搜索 pch 并阻止它尝试加载 .pch（我们刚刚删除了它）
• 将代码示例复制粘贴到 main.m 中的任何内容上
• 删除#include 的Carbon 行。 Carbon 适用于 OSX。
• 删除所有框架，添加加速框架

您可能还需要从 info.plist 中删除一个告诉项目加载 xib 的条目，但我 90% 确定您不需要为此烦恼。

注意：程序输出到控制台，结果显示为 0.000，这不是错误——它只是非常非常快

这段代码真的很晦涩难懂；评论很慷慨，但 cmets 实际上并没有让生活变得更轻松。

它的核心是：

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

在 n 个真实浮点数上进行 FFT，然后反向返回我们开始的位置。 ip 代表就地，这意味着 &A 被覆盖 这就是所有这些特殊打包错误的原因——这样我们就可以将返回值压缩到与发送值相同的空间中。

为了给出一些观点（例如：为什么我们首先要使用这个函数？），假设我们想要对麦克风输入执行音高检测，并且我们已经设置了它以便获得一些回调每次麦克风进入 1024 浮动时触发。假设麦克风采样率为 44.1kHz，则约为 44 帧/秒。

所以，我们的时间窗是 1024 个样本的持续时间，即 1/44 秒。

所以我们将 A 与来自麦克风的 1024 个浮点数打包在一起，设置 log2n=10 (2^10=1024)，预先计算一些线轴 (setupReal) 并且：

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

现在 A 将包含 n/2 个复数。这些代表 n/2 个频率区间：

• bin[1].idealFreq = 44Hz -- 即我们可以可靠检测到的最低频率是该窗口内的一个完整波，即 44Hz 波。

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

• 等

• bin[512].idealFreq = 512 * 44Hz -- 我们可以检测到的最高频率（称为奈奎斯特频率）是每对点代表一个波，即窗口内的 512 个完整波，即512 * 44Hz，或：n/2 * bin[1].idealFreq

• 实际上还有一个额外的 Bin，Bin[0]，通常称为“DC 偏移”。碰巧 Bin[0] 和 Bin[n/2] 总是有复分量 0，所以 A[0].realp 用于存储 Bin[0] 而 A[0].imagp 用于存储 Bin[ n/2]

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

因此，如您所见，它不会是一个非常出色的音高检测器，因为它没有足够精细的粒度。有一个巧妙的技巧 Extracting precise frequencies from FFT Bins using phase change between frames 可以获取给定 bin 的精确频率。

好的，现在进入代码：

注意 vDSP_fft_zrip 中的 'ip'，='in place' 即输出覆盖 A（'r' 表示它需要实际输入）

查看关于 vDSP_fft_zrip 的文档，

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

这可能是最难理解的事情。我们在整个过程中一直使用相同的容器 (&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 复数 #s 打包到 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 课程： 傅立叶理论需要很长时间才能理解（我已经断断续续地研究了好几年了）

一个cisoid是：

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

即复平面单位圆上的一点。

当您将复数相乘时，角度会相加。所以 z^k 会一直围绕单位圆跳跃； z^k 可以找到角度 k.theta

• 选择 z1 = 0+1i，即从实轴转四分之一圈，注意 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+0i, z0^(anything)=1，这是直流偏移

您可以将任何 4 点信号表示为 z0 z1 和 z2 的线性组合 即您将其投影到这些基向量上

但我听到你问“将信号投射到 cisoid 上意味着什么？”

你可以这样想：针绕着cisoid旋转，所以在样本k处，针指向k.theta方向，长度为signal[k]。与 cisoid 的频率完全匹配的信号将使结果形状在某个方向上凸出。所以如果你把所有的贡献加起来，你会得到一个强大的合成向量。 如果频率几乎匹配，则凸起会更小，并且会绕着圆圈缓慢移动。 对于与频率不匹配的信号，其贡献将相互抵消。

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

但要点是；如果我们选择将 1024 个样本投影到 {z0,...,z512}，我们将预先计算 z0 到 z512，这就是这个预先计算步骤。

---

请注意，如果您在实际代码中执行此操作，您可能希望在应用加载时执行一次，并在应用退出时调用补充释放函数。不要做很多次——它很贵。

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

值得注意的是，如果我们将 log2n 设置为例如 8，您可以将这些预先计算的值放入任何使用分辨率

现在是实际的转换，利用我们刚刚预先计算的东西：

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

此时 A 将包含 n/2 个复数，只有第一个实际上是两个伪装成复数的实数（DC 偏移量，奈奎斯特 #）。文档概述解释了这种包装。它非常简洁——基本上它允许将（复杂的）转换结果打包到与（真实的，但奇怪的打包的）输入相同的内存占用中。

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

然后再返回...我们仍然需要从 A 中解压缩原始数组。然后我们比较只是为了检查我们是否已经完全恢复了我们开始时的状态，释放我们预先计算好的线轴并完成！

但是等等！在您打开包装之前，还有最后一件事需要完成：

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

Answer 2

这是一个真实示例：c++ 的 sn-p，它使用 Accelerate 的 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);
    }
}

Answer 3

虽然我会说 Apple 的 FFT 框架很快……但您需要了解 FFT 的工作原理才能获得准确的音高检测（即计算每个连续 FFT 的相位差以找到准确的音高，而不是最主要的 bin 的间距）。

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

我正在上传一个示例 FFT 实现——敬请期待，一旦发生，我会更新此内容。

编码愉快！

Answer 4

这是另一个真实的例子： https://github.com/krafter/DetectingAudioFrequency