7 obrazków DSP555_?.GIF w tym katalogu, OSTATNI MA MIEÊ 1,5x wys i szer (3 szpalty) DSP dla kaûdego (cz. 5.) ------------------------ FFT, BASIC i ASEMBLER W tym miesiâcu nieco na luzie otwieram piâtâ czëôê cyfrowego przetwarzania sygnaîów dla kaûdego. Poniewaû mamy juû wakacje, a moûe wîaônie dlatego, zajmiemy sië baaaardzo trudnym zagadnieniem, a mianowicie praktycznym zastosowaniem algorytmu FFT dla potrzeb wizualizacji widma sygnaîu muzycznego. Koledzy poczâtkujâcy nieco sië tu zawiodâ, gdyû zagadnienia tu poruszane istotnie do îatwych nie naleûâ, ale za to odcinek ten ucieszy na pewno obeznanych nieco z tematem. No cóû, wszystkich zadowoliê sië nie da. William Mobius Pokazywanie widma sygnaîu muzycznego ma bardzo duûo zastosowaï: w oglâdaniu zakîóceï sygnaîu, wizualizacji mowy ludzkiej, analizie dúwiëku itd. Podam przykîad. Otóû podczas rekonstrukcji nagrania, pochodzâcego ze starych taôm w Polskim radiu, reûyser dúwiëku usîyszaî buczenie na tle muzyki. Poniewaû buczenie to zostaîo juû wmiksowane razem z dostarczonym nagraniem, przeto naleûaîo je, mówiâc fachowo, "wyczyôciê". Posîuûono sië zatem analizatorem widma, to jest przyrzâdem, pokazujâcym na osi czëstotliwoôci poszczególne skîadowe sygnaîu. Oglâdajâc takie widmo w studiu dla fragmentu muzyki, reûyser zauwaûyî tam dwie duûe skîadowe. Patrz rys. 2. Byîy to: 50 Hz i 100 Hz, czyli podstawowa skîadowa i jej wielokrotnoôê, charakterystyczne dla prâdu zmiennego w sieci. Innymi sîowy, kiepska aparatura elektroakustyczna, uûyta w poprzednim nagraniu, daîa o sobie znaê, poniewaû skîadowe te, czyli tzw. tëtnienie, przeniknëîo z sieci zasilajâcej do nagrania. Dziëki wyodrëbnieniu zakîócenia wystarczyîo wîâczyê w tor miksera dwa pasmowo-zaporowe filtry szpilkowe, czyli tzw. filtry klasy notch, dziëki czemu kopia nagrania pozbawiona zostaîa uciâûliwego buczenia, a samo nagranie minimalnie znieksztaîcone. Dlaczego? Wydawaê by sië mogîo, ûe sama muzyka takûe moûe zawieraê w pewnych momentach te skîadowe (instrumenty basowe), lecz braku tych dwóch czëstotliwoôci, tj. 50 Hz i 100 Hz, tak naprawdë sîuchacz nie zauwaûy (patrz rys. 4.). W rzeczywistoôci stosuje sië bardziej wymyôlne techniki usuwania zakîóceï (np. dynamiczne filtrowania, przez co w zasadzie nie ma sîyszalnej przerwy w paômie), ale pominiemy tu te sprawy, poniewaû chodziîo mi tylko o pokazanie jednego z zastosowaï zobrazowania widma sygnaîu. A przy okazji: podana tu metoda przyda sië i amigowym muzykom. Swego czasu dostaîem na gieîdzie sporo sampli, które miaîy w tle zakîócenie takim wîaônie przydúwiëkiem sieciowym (prawdopodobnie samplowane za pomocâ jakiegoô samplera firmy "no name"). Normalnie tego nie sîychaê (wiëkszoôê ludzi "puszcza sobie" moduîy na monitorze 1084 S, który ma fatalnâ fonië). Ale po zaîoûeniu dobrych sîuchawek wychodzi wszystko jak na dîoni, wiëkszoôê sampli jest "brudna" (zaszumiona, z przydúwiëkami). Dopiero dziëki odpowiedniej analizie i filtracji moûna sië byîo pozbyê tego paskudnego brzëczenia. Podana wczeôniej metoda analizy dúwiëku przez zastosowanie szeregów trygonometrycznych Fouriera, miaîa pewne ograniczenia, miëdzy innymi byîa bardzo powolna i miaîa w zasadzie zastosowanie tylko do sygnaîów okresowych, takich jak np. piîa, sinus i inne... Natomiast do analizy dúwiëków nieokresowych, a zwîaszcza takich o widmie zîoûonym, z duûâ zawartoôciâ skîadowych losowych (szumu) i innych skomplikowanych dúwiëków lepiej nadaje sië tzw. >przeksztaîcenie caîkowe Fouriera<, w wersji dyskretnej, czyli DFT (Discrete Fourier Transform), które jest podobne do szeregów Fouriera, choê otrzymany wynik jest nieco inny. W tym artykule przedstawië szybkâ odmianë tej metody, czyli tzw. metodë FFT (ang. Fast Fourier Transform -- pol. SPF Szybkie Przeksztaîcenie Fouriera), które dziëki specjalnemu algorytmowi jest przy duûej liczbie próbek ôrednio kilkadziesiât razy szybsze, i to tym wiëcej, im wiëcej jest próbek sygnaîu. Wystarczy tu wiedzieê, ûe klasyczne DFT wymaga: n^2 mnoûeï i dodawaï podczas gdy zastosowana tu wersja algorytmu FFT wymaga tylko: N*Log2(n) ************************************** DO SKÎADU: N razy logarytm przy podstawie 2 z n. dwójkë skîad powinien umieôciê nieco niûej tekstu i mniejszâ czcionkâ, tak jak sië oznacza logarytmy ************************************** *** Ma ona tylko të wadë, ûe wymaga liczby próbek, bëdâcej wynikiem szeregu N=2^p, czyli N=2,4,8,16,32,64 itd. Ale jak o tym byîo wczeôniej, w praktyce nie jest to zbyt dokuczliwe, bo moûna sygnaî analizowaê porcjami, np. po 256 bajtów, a ostatniâ koïcówkë uzupeîniê zerami do 256. Moûna takûe, w celu dokîadniejszego zobrazowania, zastosowaê tzw. Krótkoterminowâ Analizë Fouriera, co jest znacznie lepszym rozwiâzaniem niû klasyczne FFT, niemniej znowu z powodów oczywistych (artykuî zamieniîby sië w wykîad akademicki) pominë tu te sprawy. Ale co by tu nie mówiê, zysk przy obliczeniach dla FFT jest znaczny. Np. przy 256 próbkach dla klasycznego DFT naleûaîoby wykonaê: 256^2=65536 dziaîaï, a dla metody FFT tylko: 256*Log2(256)=2048 czyli róûnica jest 32-krotna. Przy wiëkszej liczbie próbek jest jeszcze wiëksza. Mówiâc obrazowo, jeûeli obliczenie widma trwaîoby dla tradycyjnej metody, powiedzmy, 10 minut, to dla FFT trwaîoby mniej niû 20 sekund. W artykule tym pokaûë metodë obliczeï widma dla 1024 próbek, czyli 512 skîadowych, a wiëc bardzo dokîadnej analizy widmowej. Co jednoczeônie, jak îatwo policzyê, da zysk przy obliczeniach, a wiëc takûe zysk czasowy: (1024^2)/(1024*Log2(1024))= ponad 102 krotny! Ûeby byîo jeszcze szybciej, zastosujemy Asembler, co sprawi, ûe moûliwe bëdzie pokazywanie widma sygnaîu nawet w czasie rzeczywistym! Ale po kolei. Kod U2 W trzecim odcinku DSP przyrównaîem tablicë liczb, uûywanâ w jëzykach wysokiego poziomu i zapisywanâ w sposób nastëpujâcy: x() do ciâgu liczb, uzyskanych z wtîoczenia jakichô danych (muzycznych czy graficznych) do wnëtrza maszyny. Wtîoczenie, czy moûe raczej przepisanie, danych do komputera moûna w wypadku danych muzycznych uzyskaê za pomocâ dowolnego samplera i jakiegoô programu muzycznego. Otrzymany tak ciâg liczb, odpowiadajâcych w przybliûeniu zapisywanemu dúwiëkowi, moûna byîo po zsamplowaniu zapisaê na dysk i dalej zrobiê z nim, co sië komu podoba. Na przykîad zaîadowaê jako dane do pamiëci, a z niej przepisaê do tablicy x(). I tu dochodzimy to maîego problemu, poniewaû jeûeli zaîadujemy fragment sampla z dysku prosto do tablicy, np. za pomocâ odczytywania kolejnych bajtów sampla przez INPUT# prosto z pliku, to okaûe sië, ûe dla jëzyków wysokiego poziomu dla ujemnych poîówek sygnaîu otrzymamy bîëdne wartoôci (patrz rys. 5.). Dzieje sië tak dlatego, ûe dúwiëk w wiëkszoôci samplerów zawodowych i komputerów (w tym na Amidze) kodowany jest za pomocâ tzw. kodu U2 czy kodu uzupeînieï do dwóch. Wyjâtkiem jest tu IBM PC, w którym dúwiëki zapisywane sâ inaczej, ale jest to wyjâtek w swojej dziedzinie nie tylko, jeôli chodzi o dúwiëk. Skupmy sië zatem na Amidze. Co jest zatem z tym kodem U2? Otóû charakterystycznâ jego cechâ jest wspólne zero dla liczb ujemnych i dodatnich. Moûna to îatwo zaobserwowaê, wpisujâc jakâô wartoôê z poziomu np. GFA do pamiëci (komendâ POKE) i oglâdajâc jâ nastëpnie za pomocâ asemblera, np. Asm One. Np. dla liczby 3-bitowej peîny zakres oômiu wartoôci wyglâda tak: liczba: U2 3 %011 2 %010 1 %001 0 %000 -1 %111 -2 %110 -3 %101 -4 %100 To byîo dla hipotetycznych liczb 3-bitowych. Trzeba przyznaê, ûe takie kodowanie liczb ze znakiem bardzo uîatwia róûnorakie operacje matematyczne. Teraz przejdúmy do tradycyjnych 8-bitów. Otóû, jeûeli byômy odczytywali z pliku kolejne bajty fali dúwiëkowej, powiedzmy o ksztaîcie piîoksztaîtnym i kolejnych wartoôciach: -4, -3, -2, -1, 0, 1, 2, 3 i przepisywali od razu do tablicy, to okaûe sie, ûe otrzymamy zamiast tego nastëpujâcy ciâg liczb: 252, 253, 254, 255, 0, 1, 2, 3 Chyba kaûdy wie dlaczego. Dla tych, którzy nie wiedzâ: liczba, np. -1, kodowana jest na w ASCII na jednym bajcie za pomocâ: %11111111 co, rzecz jasna, daje zawsze po odczytaniu 255. A wiëc kaûde dane, czy to odczytywane z pamiëci, czy teû z pliku, naleûy zamieniaê na wartoôci ze znakiem przez zanegowanie ostatniego bitu. W wypadku jëzyków wysokiego poziomu sprowadza sië to do odjëcia wartoôci 256 w wypadku wykrycia liczby wiëkszej od 127, co wyglâda mniej wiëcej tak: ' tu czytanie zmiennej 'a&' IF a&>127 SUB a&,256 ENDIF Podany tu przykîad dotyczy jednego z najlepszych BASIC-ów Amigi, czyli GFA-BASIC. Przy czym naleûy tu podkreôliê, ûe GFA, mimo mylâcej nazwy, nie jest wîaôciwie BASIC-iem i ma skîadnië bardziej zbliûonâ do Amiga E czy Pascala. Nie ma numeracji wierszy, ma za to moûliwoôê automatycznego robienia tzw. wciëê w listingu i rozbudowane pëtle: DO..LOOP, WHILE..WEND, UNIT, EXIT, poza tym znakomicie wspópracuje z systemem. Niestety, mimo naprawdë znakomitego narzëdzia, jakim jest GFA, czasami niepoprawnie dziaîa w systemach wyûszych niû 1.3, co sprawdziîem i, maîo tego, znalazîem dwa bîëdy w samym kompilatorze, które od razu poprawiîem, niemniej jednak program, znakomicie dziaîajâcy na 1.3, w wypadku 2.04 czasami dziaîa, a czasami nie. Objawia sië to miëdzy innymi tym, ûe program interpretowany dziaîa bardzo dobrze, a po kompilacji moûe sië zawiesiê, nie wiadomo dlaczego. Na 3.0 czësto teû po prostu pada bez powodu. Dlatego, mimo ûe caîy swój software pisaîem w tym jëzyku z ewentualnymi wstawkami w asemblerze (sam mam 2.04), to mam zamiar stopniowo sië przerzuciê na jakiô inny jëzyk. Podam tu teraz parë oznaczeï, dotyczâcych wîaônie GFA, poniewaû wszystkie programy bëdâ w tej konwencji. Obok podam kilka odpowiedników w innych jëzykach, co uîatwi zrozumienie póúniejszych moich listingów osobom uûywajâcym tych jëzyków: GFA AMOS C assembler Znaczenie a| brak UBYTE .b bajt (liczby 0..255) a& brak WORD .w 16-bitowa liczba ze znakiem a% A LONG .l 32-bitowa liczba ze znakiem a A# float - liczba zmiennoprzecinkowa Jak widaê, GFA stosuje, podobnie jak C, kilka róûnych rodzajów dokîadnoôci przedstawiania liczb. Róûnice sâ w tych jëzykach niewielkie, dlatego nastëpne ewentualne listingi takûe bëdâ juû tylko w konwencji jëzyka GFA lub assemblera. Z kolei jëzyk C nie potrzebuje konwersji, bo moûna sobie od razu zdefiniowaê wskaúnik do pamiëci, na której bëdziemy operowali bezpoôrednio, nie mówiâc juû o asemblerze, który jednakowo widzi liczby ze znakiem, jak i bez znaku, a wiëc odpada problem konwersji. Tak wiëc w wyniku tych dziaîaï otrzymamy z danych typu 0..255 dane sampla 8-bitowego z zakresu -128..127 przepisane do tablicy. Majâc juû w tablicy czy pamiëci (obojëtne) ciâg liczb odpowiadajâcy fragmentowi jakiegoô dúwiëku moûna przystâpiê do analizy. Ale najpierw... Trochë teorii Jest to, niestety, dawka niezbëdna, bez której ciëûko bëdzie cokolwiek dalej zrobiê i zrozumieê. Otóû najpierw przyjrzyjmy sië poprzednim rysunkom. Przedstawione jest tam widmo dúwiëku raczej schematycznie. Oô Y to amplituda poszczególnych skîadowych, a oô X to kolejne czëstotliwoôci. Wszystko gra. Tyle ûe przedstawienie takie dokonane jest w osi czëstotliwoôci (X) w konwencji logarytmicznej, co, przypominam, oznacza, ûe odstëpy (wizualnie) na wykresie sâ miëdzy 16 Hz, a 160 Hz takie same, jak miëdzy 160 Hz a 1600 Hz. Widaê wiëc, ûe wykres ma duûâ rozdzielczoôê dla tonów basowych i ôrednich, natomiast dla tonów wyûszych rozdzialczoôê ta maleje. Takie przedstawienie jest zgodne z dziaîaniem ludzkiego sîuchu, poniewaû przy rosnâcych czëstotliwoôciach ucho nasze sîyszy coraz mniej dokîadnie. Natomiast wszystkie metody analizy dúwiëku, oparte na fourierowskiej analizie, przedstawiajâ pierwotnie widmo sygnaîu w skali liniowej. Taka skala z kolei upoôledza rozdzielczoôê dla basów, natomiast bardzo dobrze pokazuje czëstotliwoôci ôrednie i wyûsze (przy czym odwrotnie do poprzedniego przykîadu -- dokîadnoôê ta roônie z czëstotliwoôciâ). I jest to pierwsza waûna róûnica, którâ naleûy zapamiëtaê. Niedîugo do tego wrócë. Teraz warto krótko wspomnieê o dziaîaniu samej analizy Fouriera, co przedstawia listing 1. Otóû w uproszczeniu za pomocâ mnoûenia sygnaîu przez zmienny w czasie czynnik zespolony powoduje sië, ûe kolejne czëstotliwoôci harmoniczne (f2, f3, f4 itd.) przesuwajâ sië na osi czëstotliwoôci niejako w tyî i stajâ sië podstawowâ czëstotliwoôciâ f1, przez co îatwo juû obliczyê dalej jej zespolonâ amplitudë. Piszë "zespolonâ" i nie naleûy baê sië tego okreôlenia, które oznacza po prostu, ûe amplituda tak wyliczonej harmonicznej nie jest opisywana za pomocâ pojedynczej liczby, lecz za pomocâ liczby zespolonej, oznaczanej zazwyczaj: z=x+jx Liczba ta zaô to po prostu dwie oddzielne liczby. Czytelnik moûe zapytaê, czy amplitudë harmonicznej moûna zapisywaê jednoczeônie dwoma róûnymi liczbami? Otóû moûna. Proszë sobie przypomnieê moj poprzedni artykuî, w którym otrzymywaliômy z szeregów Fouriera dwa równolegîe wspóîczynniki: sinusowy i cosinusowy. Tu jest podobnie. Istnieje nawet daleko posuniëta analogia miëdzy tymi parami wspóîczynników. Nastëpna czëôê programu skaluje otrzymany wynik, poniewaû na skutek kolejnych mnoûeï zespolonych amplituda harmonicznych byîaby za wielka. Dlatego stosuje sië tu dzielenie wszystkich otrzymanych liczb przez czynnik skalujâcy: SQR(n) gdzie 'n' to caîkowita liczba próbek sygnaîu. Co dla 32 próbek daje ok. 5,66, a dla np. 1024 próbek -- wîaônie 32. Innym przykîadem jest maîy podprogramik sortujâcy. Po coû on? Otóû po przeksztaîceniu sygnaîu metodâ FFT otrzymuje sië, owszem, widmo, lecz nieco... pokieîbaszone. Po prostu ze wzglëdu na samâ metodë kolejne czëstotliwoôci nie wystëpujâ po kolei tak, jak f1, f2, f3 itd., lecz sâ pomieszane. Ta procedurka po prostu odpowiednio ukîada je na miejsca. Do prezentowanego dalej programu Czytelnik sam bëdzie musiaî dorobiê odpowiednie do zastosowanego jëzyka programowania instrukcje graficzne, ale liczë na to, ûe nie bëdzie z tym wiëkszego problemu, przy czym proszë sië nie sugerowaê rysunkiem, poniewaû zrobiîem go na podstawie wîasnego dziaîajâcego programu. Przepis na FFT w przykîadach: Otóû po przyjëciu liczby próbek 'n' jako 32 i po policzeniu przykîadowego FFT otrzymujemy taki oto rezultat (patrz rys. 6.). Jak widaê z próbki 32 bajtów sampla, po przepuszczeniu przez FFT otrzymaliômy dwa dziwnie symetryczne wykresy po 32 bajty, zamiast oczekiwanego(?) jednego wykresu, skîadajâcego sië 15 harmonicznych. No i co z tym fantem zrobiê? Otóû otrzymaliômy klasyczny podrëcznikowy rezultat. Na górze mamy wspóîczynniki rzeczywiste x(), a na dole urojone jx(). Przedstawiajâ one po prostu odpowiednie skîadowe cosinusowe i sinusowe. Teraz wystarczy zignorowaê czëôê leûâcâ w "drugiej poîowie" wykresu i juû prawie jesteômy w domu. Dokîadnie rzecz biorâc, naleûy wziâê pod uwagë skîadowe, zawarte pomiëdzy: 1..(n/2)+1 a zignorowaê pozostaîe, czyli w zakresie: (n/2)+2..n oczywiôcie liczâc próbki od 1, a nie od 0 (rys. 7 A). Ale to jeszcze nie wszystko. W celu otrzymania klasycznego widma takiego jak na obrazkach, naleûy obliczyê tzw. >moduî widma zespolonego< wedîug wzoru: a=SQR((x^2)+(jx^2)) dla kolejnych skîadowych. W tym celu moûemy jeszcze zrobiê sobie nastëpnâ tablicë a(), w której umieôcimy wynik tych obliczeï. Teraz po zobrazowaniu tego, co wyszîo i przepisujâc dane, czyli amplitudë harmoniczej jako wysokoôê sîupków, otrzymamy wreszcie normalny wykres widma amplitudowego (rys.7 B), które, jak to wynika ze wzoru, bëdzie miaîo zawsze wartoôci dodatnie. W kolejnym odcinku bëdzie zamieszczony listing programu tym razem w asemblerze. ---------------------------------------------- /podpisy/ Rys.1 Prâd zmienny w sieci zasilajâcej 50 Hz (A) i jego widmo-nieco wyidealizowane (B). Przy kiepskiej aparaturze i marnych stopniach filtrujâcych w zasilaczu do sygnaîu przenika resztka (tëtnienie). Przebieg przykîadowych zakîóceï przedostajâcych sië do sygnaîu (C) i ich widmo (D). Rys.2 Kawaîek muzyki (widmo) z wyraúnymi pikami buczenia (zaznaczone na rysunku). Rysunek przedstawia problem nieco przejaskrawiony poniewaû w rzeczywistoôci tony te majâ znacznie mniejszâ amplitudë i mogâ ginâê wizualnie wôród innych dúwiëków (lecz bëdâc zarazem caîkiem sîyszalnymi). Rys.3 Dwa filtry typu notch poîâczone kaskadowo-schemat blokowy (A) i ich îâczna charakterystyka filtracyjna (B). Rys.4 Muzyka (widmo) po czëôciowym przefiltrowaniu polegajâcym na usuniëciu z pasma dwóch zakîócajâcych czëstotliwoôci. Rys.5 Poprzesuwane poîówki sygnaîu przy odczytywaniu bezpoôrednio z pliku na dysku lub komendâ PEEK() z pamiëci sâ konsekwencjâ sposobu zapisywania liczb ze znakiem. Rys.6 Po przeksztaîceniu sygnaîu przez FFT i zobrazowaniu wnëtrza tablic x() i jx() otrzymamy poczâtkowo doôê zagadkowy rysunek. Rys.7 Najpierw naleûy uciâê niepotrzebnâ czëôê wykresu (A), a nastëpnie obliczyê moduî, co da nareszcie îadny wykresik widma skîadajâcy sië z 17 sîupków (B).