Utilisation de Apple FFT et Accelerate Framework

Quelqu’un a-t-il déjà utilisé l’ Apple FFT pour une application iPhone ou savoir où trouver un exemple d’application pour l’utiliser? Je sais qu’Apple a un exemple de code affiché, mais je ne sais pas trop comment l’implémenter dans un projet réel.

Je viens de recevoir le code FFT pour un projet iPhone:

  • créer un nouveau projet
  • supprimer tous les fichiers à l’exception de main.m et xxx_info.plist
  • aller dans les parameters du projet et rechercher pch et l’empêcher d’essayer de charger un fichier .pch (vu que nous venons de le supprimer)
  • copier coller l’exemple de code sur tout ce que vous avez dans main.m
  • Retirez la ligne qui # inclut le carbone. Le carbone est pour OSX.
  • supprimer tous les frameworks, et append accélérer framework

Vous devrez peut-être également supprimer une entrée de info.plist qui indique au projet de charger un fichier xib, mais je suis sûr à 90% que vous n’avez pas à vous en préoccuper.

REMARQUE: Les sorties du programme sur la console, les résultats sont 0,000, ce n’est pas une erreur – c’est très très rapide

Ce code est vraiment bêtement obscur; il est généreusement commenté, mais les commentaires ne rendent pas la vie plus facile.

Au fond, il y a:

 vDSP_fft_zrip(setupReal, &A, ssortingde, log2n, FFT_FORWARD); vDSP_fft_zrip(setupReal, &A, ssortingde, log2n, FFT_INVERSE); 

FFT sur n flotte, puis retour pour revenir à notre sharepoint départ. ip signifie in-place, ce qui signifie que & A est écrasé C’est la raison de toute cette clé spéciale d’emballage – de sorte que nous pouvons écraser la valeur de retour dans le même espace que la valeur d’envoi.

Pour donner un aperçu (comme dans: pourquoi utiliserions-nous cette fonction en premier lieu?), Supposons que nous souhaitons effectuer une détection de hauteur de son sur l’entrée micro, et nous l’avons configurée pour que des rappels soient déclenchés à chaque fois le microphone reçoit 1024 flotteurs. En supposant que le taux d’échantillonnage du microphone était de 44,1 kHz, soit environ 44 images / sec.

Ainsi, notre fenêtre temporelle est quelle que soit la durée de 1024 échantillons, soit 1/44 s.

Donc, nous devrions emballer A avec 1024 flotteurs du micro, définir log2n = 10 (2 ^ 10 = 1024), précalculer certaines bobines (setupReal) et:

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

A maintenant contiendra n / 2 nombres complexes. Ceux-ci représentent n / 2 intervalles de fréquence:

  • bin [1] .idealFreq = 44Hz – c’est-à-dire que la fréquence la plus basse que nous pouvons détecter de manière fiable est UNE onde complète dans cette fenêtre, soit une onde de 44Hz.

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

  • etc.

  • bin [512] .idealFreq = 512 * 44Hz – La fréquence la plus élevée que nous pouvons détecter (appelée fréquence de Nyquist) est l’endroit où chaque paire de points représente une onde, soit 512 ondes complètes dans la fenêtre, soit 512 * 44Hz, ou: n / 2 * bin [1] .idealFreq

  • En fait, il y a un Bin supplémentaire, Bin [0] qui est souvent appelé “Décalage DC”. Il se trouve que Bin [0] et Bin [n / 2] auront toujours un composant complexe 0, donc A [0] .realp est utilisé pour stocker Bin [0] et A [0] .imagp est utilisé pour stocker Bin [ n / 2]

Et l’ampleur de chaque nombre complexe est la quantité d’énergie qui vibre autour de cette fréquence.

Donc, comme vous pouvez le constater, ce n’est pas un très bon détecteur de hauteur car il n’a pas une granularité suffisante. Il existe une astuce astucieuse Extraction de fréquences précises à partir de bacs FFT en utilisant le changement de phase entre les images pour obtenir la fréquence précise pour un bac donné.

Ok, maintenant sur le code:

Notez le ‘ip’ dans vDSP_fft_zrip, = ‘in place’, c.-à-d. Que la sortie écrase A (‘r’ signifie qu’il faut des entrées réelles)

Regardez la documentation sur vDSP_fft_zrip,

Les données réelles sont stockées sous une forme complexe, avec des réalités impaires stockées du côté imaginaire de la forme complexe et même des réalités stockées sur le côté réel.

C’est probablement la chose la plus difficile à comprendre. Nous utilisons le même conteneur (& A) tout au long du processus. donc au début nous voulons le remplir avec n nombres réels. après la FFT, il va contenir n / 2 nombres complexes. nous lançons ensuite cela dans la transformation inverse, et nous espérons sortir nos nombres réels n originaux.

maintenant la structure de A sa configuration pour les valeurs complexes. Ainsi, vDSP doit standardiser la manière de mettre en mémoire des nombres réels.

Nous commençons donc par générer n nombres réels: 1, 2, …, n

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

Ensuite, nous les emballons dans A comme n / 2 complexes #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, // ssortingde 2, as each complex # is 2 floats &A, 1, // ssortingde 1 in A.realP & .compP nOver2); // n/2 elts 

Vous devriez vraiment regarder comment A est alloué pour obtenir ceci, peut-être regarder COMPLEX_SPLIT dans la documentation.

 A.realp = (float *) malloc(nOver2 * sizeof(float)); A.imagp = (float *) malloc(nOver2 * sizeof(float)); 

Nous faisons ensuite un calcul préalable.


Classe DSP rapide pour les maths: la théorie de Fourier prend beaucoup de temps pour se faire une idée (ça fait plusieurs années que je regarde ça)

Un cisoid est:

 z = exp(i.theta) = cos(theta) + i.sin(theta) 

c'est-à-dire un point sur le cercle unitaire dans le plan complexe.

Lorsque vous multipliez des nombres complexes, les angles s’ajoutent. Donc, z ^ k continuera à sauter autour du cercle unitaire; z ^ k peut être trouvé sous un angle k.theta

  • Choisissez z1 = 0 + 1i, c'est-à-dire un quart de tour de l'axe réel, et notez que z1 ^ 2 z1 ^ 3 z1 ^ 4 donne un quart de tour de plus pour que z1 ^ 4 = 1

  • Choisissez z2 = -1, soit un demi-tour. aussi z2 ^ 4 = 1 mais z2 a terminé 2 cycles à ce point (z2 ^ 2 est aussi = 1). On pourrait donc considérer z1 comme la fréquence fondamentale et z2 comme la première harmonique

  • De même, z3 = le point «trois quarts de tour», c'est-à-dire que -i termine exactement 3 cycles, mais aller en avant 3/4 à chaque fois revient à reculer de 1/4 à chaque fois

c'est-à-dire que z3 est juste z1 mais dans la direction opposée

z2 est la fréquence la plus significative, car nous avons choisi 4 échantillons pour tenir une onde complète.

  • z0 = 1 + 0i, z0 ^ (n'importe quoi) = 1, ceci est le décalage DC

Vous pouvez exprimer n'importe quel signal à 4 points sous la forme d'une combinaison linéaire de z0 z1 et z2, c'est-à-dire que vous le projetez sur ces vecteurs de base.

mais je vous entends demander "que signifie projeter un signal sur un cisoid?"

Vous pouvez penser à ceci: L'aiguille tourne autour de la cisoïde, donc à l'échantillon k, l'aiguille pointe dans la direction k.theta et la longueur est un signal [k]. Un signal qui correspond exactement à la fréquence du cisoïde gonfle la forme résultante dans une certaine direction. Donc, si vous additionnez toutes les consortingbutions, vous obtiendrez un vecteur puissant. Si la fréquence correspond presque à celle, le renflement sera plus petit et se déplacera lentement autour du cercle. Pour un signal qui ne correspond pas à la fréquence, les consortingbutions s'annulent mutuellement.

http://complextoreal.com/tutorials/tutorial-4-fourier-analysis-made-easy-part-1/ vous aidera à obtenir une compréhension intuitive.

Mais l'essentiel est; Si nous avons choisi de projeter 1024 échantillons sur {z0, ..., z512}, nous aurions pré-calculé z0 à z512, et c'est ce que cette étape de pré-calcul est .


Notez que si vous le faites en code réel, vous voudrez probablement le faire une fois lorsque l'application se charge et appelez la fonction de publication complémentaire une fois qu'elle se ferme. Ne le faites pas beaucoup de fois - c'est cher.

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

Il convient de noter que si nous définissons log2n à 8, par exemple, vous pouvez lancer ces valeurs précalculées dans toute fonction fft utilisant une résolution <= 2 ^ 8. Donc, à moins d’optimiser l’optimisation de la mémoire, créez simplement un jeu pour la résolution la plus élevée dont vous aurez besoin et utilisez-le pour tout.

Maintenant, les transformations réelles, en utilisant les choses que nous venons juste de calculer:

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

A ce stade, A contiendra n / 2 nombres complexes, seul le premier est en réalité deux nombres réels (décalage DC, Nyquist #) se faisant passer pour un nombre complexe. L'aperçu de la documentation explique cet emballage. C'est plutôt simple - en gros, cela permet de compacter les résultats (complexes) de la transformation dans la même empreinte mémoire que les entrées (réelles mais étrangement empaquetées).

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

et de retour ... nous aurons toujours besoin de déballer notre tableau original de A. puis nous comparons juste pour vérifier que nous avons récupéré exactement ce que nous avons commencé avec, libérez nos bobines précalculées et faites!

Mais attendez! avant de déballer, il y a une dernière chose à faire:

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

Voici un exemple concret: un extrait de c ++ qui utilise les routines vDSP d’Accelerate pour effectuer l’autocorrélation sur l’entrée de l’unité audio IO distante. Utiliser ce framework est assez compliqué, mais la documentation n’est pas trop mauvaise.

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

Alors que je dirai que le Framework FFT d’Apple est rapide … Vous devez savoir comment fonctionne une FFT pour obtenir une détection précise de la hauteur (calcul de la différence de phase sur chaque FFT successive afin de trouver la hauteur exacte, pas la hauteur de la la plupart dominent bin).

Je ne sais pas si cela peut vous aider, mais j’ai téléchargé mon object Pitch Detector depuis mon application tuner (musicianskit.com/developer.php). Il existe un exemple de projet xCode 4 à télécharger également (vous pouvez ainsi voir comment fonctionne l’implémentation).

Je travaille à télécharger un exemple d’implémentation FFT – alors restz à l’écoute et je mettrai à jour cela une fois que cela se produira.

Heureux codage!

Voici un autre exemple concret : https://github.com/krafter/DetectingAudioFrequency