DSP dla kaûdego (cz. 6.) ------------------------ FFT, BASIC i ASEMBLER William Mobius >ZACZYNAMY UTRUDNIAÊ< Tak narysowany wykres moûna teraz zinterpretowaê. Otóû mieliômy 32 próbki. W FFT otrzymamy z nich 16 lub 17 sîupków (17 liczâc tzw. czëstotliwoôê zerowâ, czyli skîadowâ staîâ), a nie 15 harmonicznych, jak to byîo z szeregami. Tu wyjaônië przy okazji, co to jest ta skîadowa staîa, nazywana takûe czëstotliwoôciâ zerowâ (f0). Otóû przy symetrycznym wykresie (rys. 1., czëôê A) jest ona zerowa lub bliska zeru. Natomiast przy wiëkszej liczbie górek niû doîków (rys. 1., czëôê B) skîadowa ta zaczyna sië objawiaê i na wykresie sygnaîu (czasowym), i w widmie amplitudowym. Jest to zrozumiaîe, poniewaû dodajâc do sygnaîu napiëcie staîe powodujemy, ûe przesuwa sië on w kierunku wartoôci dodatnich. Mamy wiëc 16 harmonicznych. W zasadzie ostatnia skîadowa (16) zawsze bëdzie ze wzglëdu na zjawisko nakîadania sië widm (aliasing) nieco znieksztaîcona i powinno sië braê pod uwagë tylko pierwsze 15, niemniej niektórzy to dopuszczajâ (w zaleûnoôci od literatury). Skoro sâ to harmoniczne, to kolejne czëstotliwoôci powinny przedstawiaê >wielokrotnoôci< podstawowej pierwszej skîadowej f1. I tak teû jest. Na przykîad sîupek nr 5 (jeôli nie liczyê skîadowej staîej) bëdzie pokazywaî czëstotliwoôê 5 razy wiëkszâ od podstawowej czëstotliwoôci sygnaîu. I tu dochodzimy do dwóch waûnych zaleûnoôci w metodie fft, okreôlajâcych liczbë mierzalnych pasm i rozdzielczoôê czëstotliwoôciowâ. Ta pierwsza wîasnoôê zostaîa juû wîaôciwie omówiona, poniewaû im wiëcej weúmiemy próbek, tym wiëkszâ liczbë pasm bëdziemy mogli zanalizowaê. Dla 32 próbek wiëc -- 16 pasm, dla 64 próbek -- 32 pasma itd. Druga zaleûnoôê jest rozdzielczoôciâ samej metody, to znaczy okreôla dokîadnoôê w hercach, inaczej >raster< czëstotliwoôci, czyli skok, z jakim bëdâ sië zwiëkszaê kolejne skîadowe. Dochodzimy wiëc do wzoru, pozwalajâcego okreôliê liczbowo rzeczywistâ czëstotliwoôê pierwszej skîadowej i kolejnych. A wyglâda on: f1=1/Tpr Oznacza to po prostu, ûe czëstotliwoôê pierwszej skîadowej i zarazem rozdzielczoôê caîego FFT to odwrotnoôê okresu, przez jaki robiîo sië pomiar. Jeûeli np. samplowaliômy 32 próbki przez sekundë, to rozdzielczoôê metody wynosi dokîadnie 1 Hz i kolejne czëstotliwoôci bëdâ wynosiîy: f1=1Hz, f2=2 Hz, f3=3 Hz i tak aû do 16 Hz (rys. 2 A). Jest to bardzo duûa rozdzielczoôê, na ogóî nie stosowana. Takâ samâ rozdzielczoôê uzyskamy wszakûe, jeûeli przeprowadzimy np. 1024 pomiary w ciâgu sekundy. Bëdzie ona nadal wynosiîa 1 Hz, zmieni sië tylko liczba analizowanych czëstotliwoôci. Jeûeli natomiast przeprowadzimy pomiar tych 1024 próbek w ciâgu 1/10 sekundy, to zgodnie ze wzorem rozdzielczoôê czëstotliwoôciowa zmaleje do 10 Hz na pasmo, co przedstawia rys. 2 B. Teraz wiemy, ile czego zawiera sygnaî. Pzy czym wystëpuje tu ciekawa zaleûnoôê, zwana >zasadâ nieoznaczonoôci<, co zresztâ wynika z samego wzoru, poniewaû dla maîej liczby harmonicznych (np. 16) pobieramy bardzo maîe wycinki sygnaîu (32 próbki). Zatem mamy wtedy dobrâ rozdzielczoôê w czasie, lecz kiepskâ w dziedzinie czëstotliwoôci i na odwrót. Chcâc zobaczyê bardzo dokîadnie, z czego skîada sië widmo dúwiëku, pobieramy do analizy znaczny fragment muzyki przez jakiô czas, np. 4096 próbek, co niestety uôrednia w tym czasie zmiany widma, jakie by ewentualnie nastâpiîy, lecz za to moûemy zobaczyê takie widmo aû dla 2048 harmonicznych! Tutaj rozdzielczoôê czasowa jest kiepska, ale za to rozdzielczoôê czëstotliwoôciowa bardzo dobra. I jest to druga charakterystyczna cecha metody analizy widma dúwiëku. >MACHINE CODE< Juû sam algorytm FFT daî znaczne poprawienie szybkoôci obliczeï, jednakûe czîowiek wciâû coô ulepsza. Kiedyô potrzebowaîem i ja coô sobie ulepszyê, jako ûe przetwarzaîem za pomocâ FFT olbrzymie ciâgi danych, dlatego w koïcu napisaîem caîoôê w assemblerze. Otóû wystëpuje tu zasadniczy problem, polegajâcy na tym, ûe wszelkie algorytmy przetwarzajâce sygnaî, a szególnie FFT, sâ bardzo czuîe na bîëdy zaokrâgleï. Zresztâ moûna zrobiê próbë i przerobiê caîy listing 1. w poprzednim odcinku tak, aby komputer uûywaî tylko liczb caîkowitych. Wychodzâ wtedy gîupoty. Problem z przenoszeniem obliczeï zmiennoprzecinkowych do asemblera na pewno kaûdy koder dobrze zna. Sâ dwa rozwiâzania tych niedogodnoôci. Pierwsze (mniej dokîadne, ale szybsze) polega na wstëpnym mnoûeniu danych przez jakieô duûe liczby, przez co przesuwa sië pozornie przecinek i zyskujemy tym samâ wiëkszâ dokîadnoôê, a po zakoïczeniu obliczeï dzieli sië z kolei wynik przez të samâ liczbë. Naleûy tylko uwaûaê, aby liczby te nie byîy zbyt duûe, poniewaû moûe wtedy dojôê do przesterowania, czyli mówiâc cyfrowo, do przepeînienia bitowego, tzn. przekroczenia zekresu np. liczby 32-bitowej, ato jest juû znacznie powaûniejszym bîëdem. Drugie (lepsze) rozwiâzanie to zastosowanie do obliczeï koprocesora matematycznego, uûywajâcego szybkich maszynowych dziaîaï zmiennoprzecinkowych. Poniewaû jednak sam nie mam tego cudownego wynalazku zamontowanego w swojej maszynce, to nie chcâc ani go emulowaê, ani wywoîywaê bibiotek matematycznych Amigi zastosowaîem metodë pierwszâ, choê na poczâtku wychodziîy dziwne rzeczy. W koïcu znalazîem bîâd. Pierwszym stadium do przeróbki na asembler jest takie przerobienie programu, aby pracowaî tylko na liczbach caîkowitych. Gdy caîoôê dziaîa, to po prostu zamieniamy wszystko na odpowiednie komendy asemblera. Cierpliwym polecam teû metodë, polegajâcâ na skorzystaniu z bibliotek matematycznych Amigi, które zawierajâ gotowe procedury do dziaîaï i przeliczeï zmiennoprzecinkowych, choê wtedy mimo duûej dokîadnoôci i zastosowanego jëzyka maszynowego nie bëdzie to juû takie szybkie. Teraz sprawa najwaûniejsza. Program (listing 2.) ma wyliczone pewne wspóîczynniki oraz stablicowane metody sortowania, przez co jest ograniczony do obliczania 512 harmonicznych. Przy czym w asemblerze zrobiony jest tylko gîówny, najbardziej czasochîonny, rdzeï programu. Inne sprawy, takie jak np. rysowanie sîupków czy obliczenia moduîu wg podanego wczeôniej wzoru, naleûy obliczyê juû z poziomu jëzyka wysokiego poziomu, GFA. Dopiero jeûeli ktoô chce sobie zrobiê np. ôwietlny equalizer, dziaîajâcy w czasie rzeczywistym (lub prawie rzeczywistym, bo 512 harmonicznych jednak kawaîek czasu procesora zajmuje -- na A 500 jakieô 0,5 sek.), to musi dodatkowo rysowaê sîupki z poziomu asemblera oraz stablicowaê pierwiastkowanie i podnoszenie do kwadratu, co zostawiam jako pracë domowâ dla wytrwaîych maniaków. A oto program w asemblerze (listing 2.). -----------tu listing 2.-------------------- A jego wykorzystanie jest nastëpujâce. Najpierw îaduje sië do pamiëci skompilowany program, wykonywalny w asemblerze (obiekt), pod jakiô wolny adres, zapisany na przykîad w zmiennej: proc% GFA musi mieê takûe zadeklarowane dwie tablice caîkowitoliczbowe 32-bitowe ze znakiem, kaûda o dîugoôci 1024 pól. Do pierwszej x%() wpisuje sië wartoôci sygnaîu z zakresu -128..127, a drugâ jx%() wypeînia sië zerami. Nastëpnie przez poczwórne poke, czyli LPOKE dla GFA (lub LOKE dla AMOS), wpisuje sië wskaúniki tablic dla programu maszynowego, tzn. wbija sië adresy tablic x%() i jx%(), które umieszcza sië w dwóch kolejnych danych typu dîugie sîowo. Jak widaê, w listingu 2. sâ zarezerwowane dwa miejsca po 32 bity na poczâtku procedury, a wîaôciwy program maszynowy zaczyna sië pod adresem: proc%+8 Stâd program maszynowy wie, gdzie sâ te tablice, nie musi korzystaê ze stosu i moûe sobie z danymi tam umieszczonymi dowolnie poczynaê. Teraz wywoîuje sië program maszynowy komendâ: CALL proc%+8 i za moment z poziomu GFA odczytujemy z tych samych tablic widmo zespolone. Proste i skuteczne. Jeszcze tylko zamieniê na widmo amplitudowe za pomocâ podanego wczeôniej wzoru i caîoôê jest widoczna jak na patelni. A wiëc proszë teraz wpisaê listingi do komputera. Ûyczë powodzenia.