Использование Apple FFT и Accelerate Framework
Кто-нибудь использовал Apple FFT для приложения для iPhone еще или знаете, где я могу найти образец приложения о том, как его использовать? Я знаю, что в Apple опубликован пример кода, но я не совсем уверен, как реализовать его в реальном проекте.
4 ответа
Я только что получил код FFT, работающий для проекта iPhone:
- создать новый проект
- удалите все файлы, кроме main.m и xxx_info.plist
- перейти к настройкам проекта и выполнить поиск для pch, чтобы он не пытался загрузить.pch (поскольку мы только что удалили его)
- скопируйте и вставьте пример кода поверх того, что у вас есть в 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 с.
Таким образом, мы должны упаковать 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 = 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 Bins с использованием изменения фазы между кадрами, чтобы получить точную частоту для данного Bin.
Хорошо, теперь на код:
Обратите внимание на "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
Вам действительно нужно посмотреть, как А выделяется, чтобы получить это, возможно, посмотрите 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-точечный сигнал в виде линейной комбинации z0 z1 и 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. Поэтому (если вы не хотите максимальной оптимизации памяти) просто создайте один набор для самого высокого разрешения, которое вам нужно, и используйте его для всего.
Теперь фактические преобразования, использующие то, что мы только что рассчитали:
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., затем мы сравним, просто чтобы убедиться, что мы получили именно то, с чего начали, выпустили наши предварительно рассчитанные катушки и сделали!
Но ждать! перед тем как распаковать, нужно сделать одну последнюю вещь:
// 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);
Вот реальный пример: фрагмент 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);
}
}
В то время как я скажу, что Apple FFT Framework работает быстро... Вам нужно знать, как работает FFT, чтобы получить точное определение шага (т. Е. Рассчитать разность фаз на каждом последующем FFT, чтобы найти точный шаг, а не шаг большинство доминируют в бен).
Я не знаю, поможет ли это, но я загрузил свой объект Pitch Detector из своего приложения для настройки (musicianskit.com/developer.php). Также есть пример проекта xCode 4 для скачивания (чтобы вы могли увидеть, как работает реализация).
Я работаю над загрузкой примера реализации FFT - так что следите за обновлениями, и я обновлю это, как только это произойдет.
Удачного кодирования!
Вот еще один пример из реальной жизни: https://github.com/krafter/DetectingAudioFrequency