Korzystanie z platformy Apple FFT i Accelerate Framework

Właśnie dostałem kod FFT pracujący dla projektu iPhone ' a:

• utwórz nowy projekt
• usuń wszystkie pliki z wyjątkiem main.m i xxx_info.plist
• przechodząc do ustawień projektu i poszukując pch i powstrzymując go przed próbą załadowania.pch (ponieważ właśnie go usunęliśmy)
• kopiuj wklej przykład kodu na cokolwiek masz w main.m
• usuń linię, która #include ' s Carbon. Carbon jest dla OSX.
• usuń wszystkie frameworki i dodaj accelerate framework

Może być również konieczne usunięcie wpisu z info.plist, który każe projektowi załadować xib, ale jestem na 90% pewien, że nie musisz się tym przejmować.

Uwaga: Program wychodzi na konsolę, wyniki wychodzą jako 0.000 to nie jest błąd -- to po prostu bardzo szybko

Ten kod jest naprawdę głupio niejasny; jest hojnie komentowany, ale komentarze nie ułatwiają życia.

W zasadzie w samym sercu jest:

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

FFT na N prawdziwych pływakach, a następnie odwrócić, aby wrócić do miejsca, w którym zaczęliśmy. ip oznacza in-place, co oznacza, że &a zostaje nadpisane To jest powód tego całego specjalnego pakowania-tak, że możemy podzielić wartość zwracaną na tę samą przestrzeń, co wartość wysyłaną.

Aby dać jakąś perspektywę (na przykład: dlaczego mielibyśmy używać tej funkcji w pierwszej kolejności?), Powiedzmy, że chcemy wykonać detekcję wysokości na wejściu mikrofonowym i ustawiliśmy ją tak, aby niektóre oddzwanianie jest wyzwalane za każdym razem, gdy mikrofon dostaje się w 1024 pływaków. W 1998 roku, po raz pierwszy w Polsce, w Polsce i za granicą, pojawił się na rynku w 1999 roku.]}

Więc nasze okno czasowe jest bez względu na czas trwania 1024 próbek, tj. 1/44 s. {11]}

Więc spakujemy 1024 pływaki z mikrofonu, ustawimy log2n=10 (2^10=1024), precalculate some bobbins (setupReal) and:

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

Teraz A będzie zawierać n / 2 liczby zespolone. Reprezentują one częstotliwość n / 2

• Kosz [1].idealFreq = 44Hz -- czyli najniższa częstotliwość jaką możemy niezawodnie wykryć to jedna fala kompletna w tym oknie, czyli fala 44Hz.

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

• Itd.

• Kosz [512]idealFreq = 512 * 44Hz -- najwyższa częstotliwość jaką możemy wykryć (znana jako częstotliwość Nyquista) jest tam, gdzie każda para punktów reprezentuje falę, czyli 512 pełnych fal w oknie, czyli 512 * 44Hz, lub: n/2 * Kosz [1].idealFreq

• W rzeczywistości istnieje dodatkowy Bin, Bin[0], który jest często określany jako "Offset DC". Tak się składa, że Bin[0] i Bin [N / 2] zawsze będą miały złożony Składnik 0, a więc[0].realp jest używany do przechowywania Bin[0] I A [0].imagp służy do przechowywania Bin[n/2]

Więc, jak widzisz, nie byłby to zbyt dobry detektor wysokości, ponieważ nie mają prawie wystarczająco drobną ziarnistość. Istnieje sprytna sztuczka wydobywanie precyzyjnych częstotliwości z pojemników FFT za pomocą zmiany faz między klatkami, aby uzyskać dokładną częstotliwość dla danego pojemnika.

Ok, teraz do kodu:

Zwróć uwagę na " ip "w vDSP_fft_zrip, = "na miejscu", czyli wyjście nadpisuje A ("r" oznacza, że pobiera rzeczywiste wejścia)

Spójrz na dokumentację vDSP_fft_zrip,

Prawdziwe dane są przechowywane w kompleksie split forma, z nieparzystymi liczbami przechowywanymi na urojona strona kompleksu split formy, a nawet realiów w przechowywanych na prawdziwa strona.

Teraz struktura konfiguracji its dla złożonych wartości. VDSP musi więc standaryzować sposób zapakowania do niego liczb rzeczywistych.

Więc najpierw generujemy n liczb rzeczywistych: 1, 2,..., n

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

Następnie pakujemy je do a jako 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

Naprawdę trzeba by spojrzeć na to, jak A jest przydzielane, może poszukać COMPLEX_SPLIT w dokumentacji.

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

Następnie wykonujemy wstępne obliczenia.

Szybka lekcja DSP z matematyki bods: teoria Fouriera zajmuje dużo czasu, aby ogarnąć twoją głowę (patrzyłem na to od kilku lat)

cisoid to:

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

tj. punkt na okręgu jednostkowym w płaszczyźnie zespolonej.

gdy mnożysz liczby zespolone, kąty dodają. Tak więc z^k będzie skakać po okręgu jednostki; z^k można znaleźć pod kątem K. theta

• Wybierz z1 = 0 + 1i, czyli ćwierć obrotu od rzeczywistej osi i zauważ, że z1^2 Z1^3 Z1^4 każdy daje kolejne ćwierć obrotu tak, że z1^4 = 1

• Wybierz z2 = -1, czyli pół obrotu. również z2^4 = 1, ale z2 zakończył 2 cykle w tym momencie(z2^2 to również = 1). Można więc myśleć o Z1 jako o częstotliwości podstawowej, a o z2 jako o pierwszej harmonicznej

• podobnie z3 = punkt 'trzy czwarte obrotu' tzn.-ja kończę dokładnie 3 cykle, ale faktycznie Idę Do przodu po 3/4 każdy czas jest taki sam jak cofanie się 1/4 za każdym razem

tzn. z3 jest po prostu z1, ale w przeciwnym kierunku -- nazywa się aliasing

z2 to najwyższa znacząca częstotliwość, ponieważ wybraliśmy próbki 4, aby utrzymać pełną falę.

• z0 = 1 + 0i, z0^(cokolwiek)=1, to jest DC offset

można wyrazić dowolny sygnał 4-punktowy jako kombinację liniową z0 z1 i z2 tzn. rzutujesz go na te wektory bazowe

ale słyszę, jak pytasz: "co to znaczy wysyłać sygnał na cisoid?"

można myśleć o tym w ten sposób: igła obraca się wokół cisoidu, więc w próbce k igła jest skierowana w kierunku K. theta, a długość jest sygnałem [k]. Sygnał, który pasuje do częstotliwości cisoidu dokładnie wybrzuszi powstały kształt w pewnym kierunku. Więc jeśli zsumujesz wszystkie składki, otrzymasz silny wektor wypadkowy. Jeśli częstotliwość prawie pasuje, niż wybrzuszenie będzie mniejsze i będzie poruszać się powoli wokół okręgu. Dla sygnału, który nie pasuje do częstotliwości, wkłady anulują się nawzajem.

Http://complextoreal.com/tutorials/tutorial-4-fourier-analysis-made-easy-part-1/ pomoże ci uzyskać intuicyjne zrozumienie.

ale gist jest; jeśli zdecydowaliśmy się na projekt 1024 próbek na {z0,...przeliczalibyśmy Z0 do z512, i tym właśnie jest ten etap prekalkulacyjny .

Zauważ, że jeśli robisz to w prawdziwym kodzie, prawdopodobnie chcesz to zrobić raz, gdy aplikacja ładuje się i wywołać uzupełniającą funkcję release raz, gdy kończy pracę. Nie rób tego wiele razy - to jest drogie.

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

Warto zauważyć, że jeśli ustawimy log2n na np. 8, można te wstępnie obliczone wartości wrzucić do dowolnej funkcji fft, która używa rozdzielczości

Teraz rzeczywiste transformacje, korzystając z rzeczy, które właśnie obliczyliśmy:

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

W tym punkcie A będzie zawierać n / 2 liczby zespolone, tylko pierwsza z nich to w rzeczywistości dwie liczby rzeczywiste (przesunięcie DC, Nyquist #) maskujące się jako liczba zespolona. Omówienie dokumentacji wyjaśnia to opakowanie. Jest dość schludny - w zasadzie pozwala na (złożone) wyniki przekształć, aby zostać spakowanym do tego samego śladu pamięci, co (rzeczywiste, ale dziwnie zapakowane) wejścia.

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

I z powrotem... nadal będziemy musieli rozpakować naszą oryginalną tablicę z A. Następnie porównujemy tylko po to, aby sprawdzić, czy mamy z powrotem dokładnie to, co zaczęliśmy, zwolnić nasze wstępnie obliczone szpulki i gotowe!

Ale czekaj! przed rozpakowaniem jest jeszcze jedna rzecz do zrobienia: {]}

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