Использование Apple FFT и Accelerate Framework

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

Я только что получил код FFT для проекта iPhone:

  • создать новый проект
  • удалите все файлы, кроме main.m и xxx_info.plist
  • перейти к настройкам проекта и найти pch и остановить его от попытки загрузить .pch (видя, что мы только что удалили его)
  • скопируйте пример кода над всем, что у вас есть в main.m.
  • удалите строку, в которой # include 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); 

FFT на n реальных поплавках, а затем назад, чтобы вернуться туда, где мы начали. ip означает «на месте», что означает, что «A» перезаписывается. Это причина для всей этой специальной упаковки malarkey – чтобы мы могли выдавить возвращаемое значение в то же пространство, что и значение отправки.

Чтобы дать некоторую перспективу (например, как в: почему мы должны использовать эту функцию в первую очередь?), Предположим, мы хотим выполнить определение высоты тона на микрофонном входе, и мы настроили его так, чтобы каждый обратный вызов запускался каждый раз микрофон получает 1024 поплавков. Предположим, что частота дискретизации микрофона составляет 44,1 кГц, поэтому это ~ 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 – т.е. самая низкая частота, которую мы можем надежно обнаружить, – это ОДНА полная волна внутри этого windows, то есть волна 44 Гц.

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

  • и т.п.

  • bin [512] .idealFreq = 512 * 44Hz. Самая высокая частота, которую мы можем обнаружить (так называемая частота Найквиста), – это где каждая пара точек представляет собой волну, т.е. 512 полных волн в окне, т.е. 512 * 44 Гц или: 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 [ п / 2]

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

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

Итак, Теперь на код:

Обратите внимание на «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)); 

Затем мы делаем предварительный расчет.


Быстрый class DSP для maths bods: теория Фурье занимает много времени, чтобы окунуться в голову (я смотрел на нее в течение нескольких лет)

Цисоид:

 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-точечный сигнал как линейную комбинацию z0 z1 и z2, т. Е. Вы проецируете его на эти базисные векторы

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

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

http://complextoreal.com/tutorials/tutorial-4-fourier-analysis-made-easy-part-1/ поможет вам получить интуитивное понимание.

Но суть такова; если мы выбрали проект 1024 образцов на {z0, ..., z512}, мы бы предварительно вычислили z0 через z512, и это то, что это шаг предварительного расчета.


Обратите внимание, что если вы делаете это в реальном коде, вы, вероятно, захотите сделать это один раз, когда приложение загрузится и вызовет дополнительную функцию release один раз, когда он уйдет. НЕ делайте это много раз - это дорого.

 // 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. Таким образом (если вы не хотите, чтобы оптимизация памяти была максимальной), просто создайте один набор для максимального разрешения, который вам понадобится, и используйте его для всего.

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

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

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

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

Вот пример реального мира: fragment c ++, который использует подпрограммы fc для ускорения vdSP Accelerate, чтобы выполнить автокорреляцию на входе аудиоустройства удаленного ввода-вывода. Использование этой структуры довольно сложно, но документация не так уж плоха.

 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; iimNumberBuffers; ++ii) { bzero(ioData->mBuffers[ii].mData, ioData->mBuffers[ii].mDataByteSize); } } 

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

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

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

Счастливое кодирование!

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

Давайте будем гением компьютера.