Utilizzando Apple FFT e Accelerate Framework

Qualcuno ha ancora usato la Apple FFT per un’app per iPhone o sa dove posso trovare un’applicazione di esempio su come usarla? So che Apple ha pubblicato qualche codice di esempio, ma non sono sicuro di come implementarlo in un progetto reale.

Ho appena ricevuto il codice FFT per un progetto iPhone:

  • crea un nuovo progetto
  • elimina tutti i file ad eccezione di main.m e xxx_info.plist
  • andando a proiettare le impostazioni e cercare pch e impedirgli di caricare un file .pch (visto che l’abbiamo appena cancellato)
  • copia incolla l’esempio di codice su qualsiasi cosa tu abbia in main.m
  • rimuovi la riga che # include Carbon. Il carbonio è per OSX.
  • elimina tutti i framework e aggiungi framework accelerato

Potrebbe anche essere necessario rimuovere una voce da info.plist che indica al progetto di caricare uno xib, ma sono sicuro al 90% che non è necessario preoccuparsi di ciò.

NOTA: Programmate le uscite su console, i risultati sono 0.000 e non è un errore – è molto molto veloce

Questo codice è davvero stupidamente oscuro; è generosamente commentato, ma i commenti in realtà non rendono la vita più facile.

Fondamentalmente al centro di esso è:

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

FFT su n float reali, e poi retromarcia per tornare al punto in cui siamo partiti. ip sta per in-place, il che significa che & A viene sovrascritto Questa è la ragione di tutto questo speciale malarkey di impacchettamento – in modo che possiamo schiacciare il valore di ritorno nello stesso spazio del valore di invio.

Per dare un po ‘di prospettiva (come, come in: perché dovremmo usare questa funzione in primo luogo?), Diciamo che vogliamo eseguire il rilevamento del pitch sull’ingresso del microfono, e lo abbiamo impostato in modo che ogni richiamata venga triggersta ogni volta il microfono arriva in 1024 galleggianti. Supponendo che la frequenza di campionamento del microfono fosse 44,1 kHz, quindi ~ 44 fotogrammi / sec.

Quindi, la nostra finestra temporale è qualunque sia la durata di 1024 campioni, ovvero 1/44 s.

Quindi dovremmo impacchettare A con 1024 float dal microfono, impostare log2n = 10 (2 ^ 10 = 1024), precalcolare alcune bobine (setupReal) e:

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

Ora A conterrà n / 2 numeri complessi. Questi rappresentano n / 2 bin di frequenza:

  • bin [1] .idealFreq = 44Hz – ovvero La frequenza più bassa che possiamo rilevare in modo affidabile è UNA onda completa all’interno di quella finestra, ovvero un’onda a 44Hz.

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

  • eccetera.

  • bin [512] .idealFreq = 512 * 44Hz – La frequenza più alta che possiamo rilevare (nota come frequenza di Nyquist) è dove ogni coppia di punti rappresenta un’onda, ovvero 512 onde complete all’interno della finestra, ovvero 512 * 44Hz, o: n / 2 * bin [1] .idealFreq

  • In realtà c’è un Bin, Bin [0] in più, che viene spesso definito come “DC Offset”. Succede che Bin [0] e Bin [n / 2] avranno sempre componenti complessi 0, quindi A [0] .realp è usato per memorizzare Bin [0] e A [0] .imagp è usato per memorizzare Bin [ n / 2]

E l’ quadro di ciascun numero complesso è la quantità di energia che vibra attorno a quella frequenza.

Quindi, come potete vedere, non sarebbe un rilevatore di pitch molto grande in quanto non ha granularità abbastanza fine. C’è un astuto trucco Estrazione di frequenze precise da FFT Bins usando il cambio di fase tra i frame per ottenere la frequenza precisa per un determinato bin.

Ok, ora sul codice:

Nota l”ip’ in vDSP_fft_zrip, = ‘in place’, cioè l’output sovrascrive A (‘r’ significa che prende input reali)

Guarda la documentazione su vDSP_fft_zrip,

I dati reali sono archiviati in forms complesse suddivise, con reali dispari memorizzati sul lato immaginario della forma complessa divisa e anche reali memorizzati nel lato reale.

questa è probabilmente la cosa più difficile da capire. Stiamo utilizzando lo stesso contenitore (& A) per tutto il processo. quindi all’inizio vogliamo riempirlo con n numeri reali. dopo la FFT sarà in possesso di n / 2 numeri complessi. quindi lo gettiamo nella trasformazione inversa, e speriamo di ottenere i nostri numeri originali e reali.

ora la struttura di A è la sua impostazione per valori complessi. Quindi vDSP ha bisogno di standardizzare come comprimere i numeri reali in esso.

quindi prima generiamo n numeri reali: 1, 2, …, n

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

Quindi li impacchettiamo in A come n # 2 complessi:

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

Avresti davvero bisogno di vedere come viene assegnato A per ottenere questo, magari cercare COMPLEX_SPLIT nella documentazione.

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

Quindi eseguiamo un pre-calcolo.


Classe DSP rapida per i bod di matematica: la teoria di Fourier richiede molto tempo per capirlo (l'ho guardata avanti e indietro per diversi anni)

Un cisoid è:

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

cioè un punto sul cerchio unitario nel piano complesso.

Quando moltiplichi i numeri complessi, gli angoli aggiungono. Quindi z ^ k continuerà a saltare intorno al cerchio dell'unità; z ^ k può essere trovato ad un angolo k.theta

  • Scegli z1 = 0 + 1i, ovvero un quarto di giro dall'asse reale, e nota che z1 ^ 2 z1 ^ 3 z1 ^ 4 danno ciascuno un altro quarto di giro in modo che z1 ^ 4 = 1

  • Scegli z2 = -1, ovvero un mezzo giro. anche z2 ^ 4 = 1 ma z2 ha completato 2 cicli a questo punto (z2 ^ 2 è anche = 1). Quindi puoi pensare a z1 come frequenza fondamentale e z2 come prima armonica

  • Allo stesso modo z3 = il punto "tre quarti di una rivoluzione" vale a dire -i completa esattamente 3 cicli, ma in realtà andando avanti di 3/4 ogni volta equivale a retrocedere di 1/4 ogni volta

cioè z3 è solo z1 ma nella direzione opposta - Si chiama aliasing

z2 è la più alta frequenza significativa, in quanto abbiamo scelto 4 campioni per tenere un'onda completa.

  • z0 = 1 + 0i, z0 ^ (nulla) = 1, questo è l'offset DC

Puoi esprimere qualsiasi segnale a 4 punti come una combinazione lineare di z0 z1 e z2, ovvero lo stai proiettando su questi vettori di base

ma ti sento chiedere "cosa significa proiettare un segnale su un cisoide?"

Puoi pensare in questo modo: l'ago gira intorno al cisoide, quindi al campione k, l'ago punta in direzione k.theta e la lunghezza è segnale [k]. Un segnale che corrisponda esattamente alla frequenza del cisoide farà rigonfiare la forma risultante in qualche direzione. Quindi se sommi tutti i contributi, otterrai un vettore risultante forte. Se la frequenza è quasi uguale, il rigonfiamento sarà più piccolo e si muoverà lentamente attorno al cerchio. Per un segnale che non corrisponde alla frequenza, i contributi si annulleranno a vicenda.

http://complextoreal.com/tutorials/tutorial-4-fourier-analysis-made-easy-part-1/ ti aiuterà ad ottenere una comprensione intuitiva.

Ma l'essenza è; se abbiamo scelto di proiettare 1024 campioni su {z0, ..., z512} avremmo calcificato il precalcolo da z0 a z512, ed è proprio questo il passo di precalcolo .


Nota che se stai facendo questo in codice reale probabilmente lo vuoi fare una volta quando l'app carica e chiama la funzione di rilascio complementare una volta quando si chiude. NON farlo molte volte: è costoso.

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

Vale la pena notare che se impostiamo log2n ad es. 8, potete lanciare questi valori precalcolati in qualsiasi funzione fft che usi la risoluzione <= 2 ^ 8. Quindi (a meno che tu non voglia la massima ottimizzazione della memoria) crea solo un set per la risoluzione più alta di cui avrai bisogno e usalo per tutto.

Ora le trasformazioni effettive, facendo uso delle cose che abbiamo appena precalcolato:

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

A questo punto A contiene n / 2 numeri complessi, solo il primo è in realtà due numeri reali (DC offset, Nyquist #) mascherati da un numero complesso. La panoramica della documentazione spiega questo imballaggio. È abbastanza pulito - in pratica consente ai risultati (complessi) della trasformazione di essere impacchettati nello stesso spazio di memoria degli input (reali, ma stranamente impacchettati).

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

e di nuovo ... avremo ancora bisogno di decomprimere il nostro array originale da A. quindi confrontiamo solo per verificare che abbiamo recuperato esattamente ciò che abbiamo iniziato, rilasciamo le nostre bobine precalcolate e fatto!

Ma aspetta! prima di disimballare, c'è un'ultima cosa che deve essere fatta:

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

Ecco un esempio del mondo reale: uno snippet di c ++ che utilizza le routine di fft di vDSP di Accelerate per eseguire l’autocorrelazione sull’input dell’unità audio remoto IO. L’utilizzo di questo framework è piuttosto complicato, ma la documentazione non è male.

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

Mentre dirò che il FFT Framework di Apple è veloce … Devi sapere come funziona un FFT per ottenere un’accurata rilevazione del pitch (cioè calcolare la differenza di fase su ogni FFT successiva per trovare l’intonazione esatta, non il pitch del la maggior parte dei domini bin).

Non so se è di alcun aiuto, ma ho caricato il mio object Pitch Detector dalla mia app di sintonizzazione (musicianskit.com/developer.php). C’è anche un esempio di progetto xCode 4 per il download (così puoi vedere come funziona l’implementazione).

Sto lavorando per caricare un’implementazione FFT di esempio – quindi rimanete sintonizzati e aggiornerò ciò una volta che ciò accadrà.

Buona programmazione!

Ecco un altro esempio del mondo reale: https://github.com/krafter/DetectingAudioFrequency