Использование Apple FFT и Accelerate Framework - PullRequest
69 голосов
/ 03 августа 2010

Кто-нибудь уже использовал Apple FFT для приложения для iPhone или знаете, где я могу найти пример приложения, как его использовать? Я знаю, что в Apple опубликован пример кода, но я не совсем уверен, как реализовать его в реальном проекте.

Ответы [ 4 ]

133 голосов
/ 04 ноября 2010

Я только что получил код FFT, работающий для проекта iPhone:

  • создать новый проект
  • удалить все файлы, кроме main.m и xxx_info.plist
  • идет к настройкам проекта и ищет pch и не дает возможности загрузить .pch (видя, как мы только что удалили его)
  • copy, вставьте пример кода поверх того, что у вас есть в main.m
  • удалить строку, которая # включает углерод. Углерод для 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 реальных числах с плавающей точкой, а затем обратно, чтобы вернуться туда, откуда мы начали. ip означает «на месте», что означает, что & A перезаписывается Вот причина всей этой специальной упаковочной малярии - чтобы мы могли сдавить возвращаемое значение в то же пространство, что и отправляемое значение.

Чтобы дать некоторую перспективу (как, например, в: почему мы будем использовать эту функцию в первую очередь?), Скажем, мы хотим выполнить определение высоты тона на микрофонном входе, и мы настроили его так, чтобы некоторый обратный вызов получал срабатывает каждый раз, когда микрофон попадает в 1024 поплавка. Предположим, частота дискретизации микрофона составляла 44,1 кГц, то есть ~ 44 кадра / с.

Итак, наше временное окно равно длительности 1024 выборок, т. Е. 1/44 с.

Таким образом, мы должны упаковать А с 1024 поплавками из микрофона, установить log2n = 10 (2 ^ 10 = 1024), предварительно рассчитать несколько катушек (setupReal) и:

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

Теперь A будет содержать n / 2 комплексных чисел. Они представляют n / 2 частотных бинов:

  • bin [1] .idealFreq = 44 Гц - т.е. самая низкая частота, которую мы можем надежно обнаружить, это ОДНА полная волна в этом окне, то есть волна 44 Гц.

  • bin [2] .idealFreq = 2 * 44 Гц

  • и т.д.

  • bin [512] .idealFreq = 512 * 44 Гц - самая высокая частота, которую мы можем обнаружить (известная как частота Найквиста), - это когда каждая пара точек представляет волну, то есть 512 полных волн в окне, т.е. 512 * 44 Гц или: n / 2 * bin [1] .idealFreq

  • На самом деле существует дополнительная корзина [0], которую часто называют «смещение постоянного тока». Бывает, что Bin [0] и Bin [n / 2] всегда будут иметь комплексный компонент 0, поэтому A [0] .realp используется для хранения Bin [0], а A [0] .imagp используется для хранения Bin [ п / 2]

А величина каждого комплексного числа - это количество энергии, колеблющейся вокруг этой частоты.

Итак, как вы можете видеть, это был бы не очень хороший детектор основного тона, так как он не имел почти достаточно тонкой детализации. Существует хитрый трюк Извлечение точных частот из бункеров FFT с использованием изменения фазы между кадрами для получения точной частоты для данного бина.

Хорошо, теперь на код:

Обратите внимание на 'ip' в vDSP_fft_zrip, = 'на месте', т.е. выход перезаписывает A ('r' означает, что он принимает реальные входные данные)

Посмотрите документацию по vDSP_fft_zrip,

Реальные данные хранятся в разделенном комплексе форма со странными реалами, хранящимися на мнимая сторона расколотого комплекса Форма и даже реальная в хранится на реальная сторона.

это, наверное, самая трудная вещь для понимания. Мы используем один и тот же контейнер (& A) на протяжении всего процесса. поэтому вначале мы хотим заполнить его n действительными числами. после БПФ он будет содержать n / 2 комплексных чисел. затем мы бросаем это в обратное преобразование и надеемся получить наши исходные n действительных чисел.

Теперь структура А его настройки для комплексных значений.Поэтому vDSP необходимо стандартизировать, как упаковывать в нее действительные числа.

, поэтому сначала мы генерируем n действительных чисел: 1, 2, ..., n

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

Затем мы упаковываем их вA как n / 2 complex #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 для математических тел: Теория Фурье требует много времени для полученияВаша голова вокруг (я смотрю на это время от времени несколько лет)

Цитата:

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 ^ (что угодно) = 1, это смещение постоянного тока

Вы можете выразить любой 4-точечный сигнал в виде линейной комбинации z0z1 и z2 т.е. вы проецируете его на эти базисные векторы

но я слышу, как вы спрашиваете, "что значит проецировать сигнал на цисоид?"

Вы можете думать об этом так: игла вращается вокруг цисоида, поэтому в образце k стрелка указывает в направлении k.theta, а длина - сигнал [k].Сигнал, который точно соответствует частоте цисоида, будет выпучивать результирующую форму в некотором направлении.Таким образом, если вы сложите все вклады, вы получите сильный результирующий вектор.Если частота почти совпадает, то выпуклость будет меньше и будет медленно двигаться по кругу.Для сигнала, который не соответствует частоте, вклады отменят друг друга.

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, вы можете бросить эти предварительно рассчитанные значения в любую функцию fft, которая использует разрешение <=2 ^ 8.Поэтому (если вы не хотите максимальной оптимизации памяти) просто создайте один набор для самого высокого разрешения, которое вам нужно, и используйте его для всего. </p>

Теперь выполняем фактические преобразования, используя то, что мы только что рассчитали:

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

В этой точке A будет содержать n / 2 комплексных чисел, только первое фактически является двумя действительными числами (DC offset, Nyquist #), маскирующимися под комплексное число.Обзор документации объясняет эту упаковку.Это довольно аккуратно - в основном это позволяет (сложные) результаты преобразования быть упакованы в тот же объем памяти, что и (реальные, но странно упакованные) входы.

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

и обратно ... нам все еще нужно будет распаковать наш исходный массив из A. Затем мы сравним, чтобы убедиться, что мы вернулись именно с того, с чего начали, выпустить наши предварительно рассчитанные катушки и сделать!1182 *

Но подождите!перед тем как распаковать, нужно сделать одну последнюю вещь:

// 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);
26 голосов
/ 21 августа 2010

Вот пример из реальной жизни: фрагмент кода c ++, который использует подпрограммы vDSP Accelerate для выполнения автокорреляции на входе аудиоустройства 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);
    }
}
13 голосов
/ 15 июня 2012

Хотя я скажу, что Apple FFT Framework работает быстро ... Вам нужно знать, как работает FFT, чтобы получить точное определение высоты тона (т. Е. Рассчитать разность фаз для каждого последующего FFT, чтобы найти точный шаг, а не шаг самого доминирующего бен).

Не знаю, поможет ли это, но я загрузил свой объект Pitch Detector из своего приложения для настройки (musicianskit.com/developer.php). Существует также пример проекта xCode 4 для загрузки (так что вы можете увидеть, как работает реализация).

Я работаю над загрузкой примера реализации FFT - так что следите за обновлениями, и я обновлю это, как только это произойдет.

Удачного кодирования!

4 голосов
/ 28 апреля 2014

Вот еще один пример из реальной жизни: https://github.com/krafter/DetectingAudioFrequency

...