使用Apple FFT和加速框架

Question

有没有人用过Apple FFTiPhone应用程序或知道我在哪里可以找到一个如何使用它的示例应用程序？我知道Apple发布了一些示例代码,但我不确定如何将它实现到实际项目中.

Answer 1

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

• 创建一个新项目
• 删除除main.m和xxx_info.plist以外的所有文件
• 转到项目设置并搜索pch并阻止它尝试加载.pch(看到我们刚刚删除它)
• 将代码示例复制粘贴到main.m中的任何内容
• 删除#include's 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被覆盖这就是所有这些特殊包装malarkey的原因 - 这样我们就可以将返回值压缩到与发送值相同的空间.

为了给出一些观点(例如:为什么我们首先要使用这个函数？),假设我们想要对麦克风输入执行音高检测,我们已经设置它以便每次触发一些回调麦克风有1024个浮点数.假设麦克风采样率为44.1kHz,那么~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个频率箱:

• 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 Offset'.碰巧Bin [0]和Bin [n/2]将始终具有复杂的组件0,因此A [0] .realp用于存储Bin [0]和A [0] .imagp用于存储Bin [ N/2]

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

因此,正如您所看到的,它不是一个非常棒的音调检测器,因为它没有足够精细的粒度.有一个狡猾的技巧 使用帧之间的相位变换从FFT区提取精确的频率,以获得给定区间的精确频率.

好的,现在进入代码:

注意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复合体#s:

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

• 选择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 ^(任何)= 1,这是DC偏移

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

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

你可以这样想:针绕着cisoid旋转,所以在样本k,针指向方向k.theta,长度是信号[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复数,只有第一个实际上是两个实数(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 ++片段,它使用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);
    }
}

Answer 3

虽然我会说Apple的FFT框架很快......你需要知道FFT是如何工作的,以便获得准确的音高检测(即计算每个连续FFT的相位差,以便找到精确的音高,而不是音高.最主导的bin).

我不知道它是否有任何帮助,但我从我的调谐器应用程序(musicianskit.com/developer.php)上传了我的Pitch Detector对象.还有一个用于下载的xCode 4项目示例(因此您可以看到实现的工作原理).

我正在努力上传一个示例FFT实现 - 所以请继续关注,一旦发生这种情况我会更新.

快乐的编码!