logo标签列表

使用苹果FFT和Accelerate框架

iphoneaudiosignal-processingfftaccelerate-framework
73

有人在iPhone应用程序中使用过Apple FFT吗？或者知道我可以在哪里找到如何使用它的示例应用程序吗？我知道苹果公司已经发布了一些示例代码，但我不太确定如何将其实现到实际项目中。

26
好的，文件说明太糟糕了。你说得对。 - P i
1
@Pi 特别是关于特殊数据排序的部分 - 实际上在许多情况下并不适用。 - marko
4个回答
139

我刚刚为一个iPhone项目编写的FFT代码已经可以使用：

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

您可能还需要从info.plist中删除一个条目，该条目告诉项目加载xib，但我90％确定您不需要担心这个问题。

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

这段代码真的非常晦涩难懂；它有大量注释，但这些注释实际上并没有使生活变得更轻松。

基本上，它的核心是：

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），预先计算一些bobbins（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]

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

因此，正如您所看到的，它不会成为非常好的音高检测器，因为它没有足够精细的细分度。有一个巧妙的技巧使用帧之间的相位变化从FFT Bin中提取精确频率来获取给定Bin的精确频率。

好了，现在进入代码：

请注意vDSP_fft_zrip中的'ip'，表示“就地”即输出覆盖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);

接下来我们将它们打包成 A，作为 n/2 个复数：

// 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.theta 处找到

  • 选择 z1 = 0+1i，即从实轴旋转1/4圈，并注意 z1^2 z1^3 z1^4 每个都会再旋转1/4圈，所以 z1^4 = 1

  • 选择 z2 = -1，即旋转1/2圈。也有 z2^4 = 1，但此时 z2 已完成了两个周期（z2^2 也等于 1）。因此，您可以将 z1 视为基本频率，将 z2 视为第一谐波

  • 同样，z3 = “三分之四的旋转”点，即 -i 正好完成了 3 个周期，但每次向前走 3/4 实际上等于每次向后走 1/4

即 z3 只是 z1 的相反方向 -- 这被称为别名效应

z2 是最高有效频率，因为我们选择了 4 个样本来表示一个完整的波。

  • z0 = 1+0i，z0^(任何数)=1，这是直流偏置

你可以用z0 z1和z2的线性组合来表示任何4点信号 也就是说，你将其投影到这些基向量上

但我听到你在问“将信号投影到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，您可以将这些预先计算的值投入任何使用分辨率<=2 ^ 8的fft函数中。因此（除非您想要最终的内存优化），只需为您需要的最高分辨率创建一个集合，并将其用于所有内容。
现在，实际的转换，利用我们刚刚预先计算的东西：
vDSP_fft_zrip(setupReal, &A, stride, log2n, FFT_FORWARD);

此时 A 将包含 n/2 个复数，只有第一个实际上是两个伪装成复数的实数（直流偏移量、奈奎斯特数字）。文档概述解释了这种打包方式。它非常巧妙 -- 基本上允许将变换的（复数）结果打包到与（实数但奇怪地打包的）输入相同的内存占用中。
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);
这不是44，而是43！在更高的二进制中，这非常重要！22050/512 = 43！ - Curnelious
1
深入解释。你能否发布一下这个苹果链接呢？我搜索了但是它引导我到多个样例，我真的想通过你的解释来理解它。谢谢！ - Nirav Bhatt
1
这是一篇很棒的帖子。是否有一个可供逐步调试代码的 GitHub 项目？ - Michael
1
你好。我们能在哪里看到完整的代码吗？我找不到这里提到的Apple示例。谢谢。 - Andrei Filip
26

以下是一个真实案例：使用Accelerate的vDSP fft例程在Remote 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);
    }
}
1
很好的例子，但你能告诉我这两个函数的实现方向吗：getMaxFramesPerSlice() 和 quadInterpolate()？ - CJ Hanson
1
@CJ：看起来 getMaxFramesPerSlice() 是检索每次回调触发时发送的帧数。我认为这同样可以是一个 #define。 - P i
3
这是一个简单的音高检测算法，使用输入信号的自相关。在这种情况下，“getMaxFramesPerSlice()”不能被定义为“#define”，因为它可能会随着每次运行而变化。这个方法实际上是对应音频单元属性访问器的包装器。代码将输入清零，因为相同的缓冲区传递给设备的输出，清零可以防止反馈环路。 - Art Gillespie
1
我认为不应该对元素0应用vDSP_zvmags，因为它的虚部实际上是奈奎斯特桶的实部。难道你不只是要平方A.realp[0]和A.imagp[0]，而不是清零A.imagp[0]吗？ - pat
根据这篇文章：http://www.dspguru.com/dsp/howtos/how-to-interpolate-fft-peak，以及我在该网站上看到的另一个自动校正示例中使用的插值函数，我认为您的quadinterpolate()函数名称是不准确的，因为您将结果添加到索引中，而插值函数实际上只计算三个项目，如下所示：`float quadinterpolate(float *p) { float y1 = *p++, y2 = *p++, y3 = *p++; return (y3 - y1) / (2 * (2 * y2 - y1 - y3)); }` - clearlight
显示剩余4条评论
14

虽然我会说苹果的FFT框架很快...但你需要知道FFT的工作原理才能进行准确的音高检测(即在每个连续的FFT上计算相位差以找到精确的音高，而不是最强信号的音高)。

我不知道这是否有帮助，但我上传了我的音高检测器对象（来自我的调音应用程序musicianskit.com/developer.php）。还提供了一个示例xCode 4项目进行下载（这样您就可以看到实现的工作方式）。

我正在努力上传一个FFT实现的示例 - 因此请继续关注，我会更新这篇文章。

Happy coding!

感谢您的分享，但是您的示例代码存在以下错误无法编译通过：1）错误：'interp'类型冲突[3]。2）文件Auto Correllation/Auto Correllation/AudioController.m:92:32，错误：未声明的标识符 'recordingCallback'[3]。 - Meir
2
抱歉，我有点懒......但这是项目链接：https://github.com/kevmdev/PitchDetectorExample。它应该可以正确编译（至少在几周前我上次尝试时是这样的），但今晚我会再次检查！ - Kpmurphy91
4
回答您的问题 - 是的，您可以。在输出数组中，您有索引作为频率，值作为幅度。因此，第一个元素是最低频率，最后一个元素是最高频率（或反之亦然）。谢谢！ - krafter
好的，太棒了。我仍然只想看看是否可以在出现其他更显著噪音的同时检测到高频率，例如18000hz。不确定是否可能？在ViewController.mm中这个函数内，maxIndex代表频谱中发现的最高频率吗？ static Float32 strongestFrequencyHZ(Float32 *buffer, FFTHelperRef *fftHelper, UInt32 frameSize, Float32 *freqValue) - Alan Scarpa
要在存在其他强频率的情况下检测18000Hz需要进行一些修改。您可以仅限于听取16000-18000Hz范围内的频率范围，并查看其内部的高峰和低谷。 - krafter
嘿，Krafter，你的例子很好。你能帮我调一下我的吉他调音器吗？ - tryKuldeepTanwar
感谢@krafter分享代码。它确实为一个困难的话题提供了很好的见解。 - 02fentym
显示剩余5条评论
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接