Amiga MA Magazine 1997 #3/amigamamagazinepolishissue03-1/ma_1995/07/ami26

Text File | 1997-04-15 | 17KB | 391 lines

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

<lead>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.

<a>William Mobius

<txt>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:

<l>n^2

<txt>mnoûeï i dodawaï podczas gdy zastosowana tu wersja algorytmu FFT
wymaga tylko:

<l>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ê:

<l>256^2=65536

<txt>dziaîaï, a dla metody FFT tylko:

<l>256*Log2(256)=2048

<txt>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:

<l>(1024^2)/(1024*Log2(1024))= ponad 102 krotny!

<txt>Û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.

<sr>Kod U2

<txt>W trzecim odcinku DSP przyrównaîem tablicë liczb, uûywanâ w
jëzykach wysokiego poziomu i zapisywanâ w sposób nastëpujâcy:

<l>x()

<txt>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:

<r>
liczba:      U2

 3      %011

 2      %010

 1      %001

 0      %000

-1      %111

-2      %110

-3      %101

-4      %100

<txt>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:

<l>-4, -3, -2, -1, 0, 1, 2, 3

<txt>i przepisywali od razu do tablicy, to okaûe sie, ûe
otrzymamy zamiast tego nastëpujâcy ciâg liczb:

<l>252, 253, 254, 255, 0, 1, 2, 3

<txt>Chyba kaûdy wie dlaczego. Dla tych, którzy nie wiedzâ:
liczba, np.  -1, kodowana jest na w ASCII na jednym bajcie za
pomocâ:

<l>%11111111

<txt>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:

<l>
' tu czytanie zmiennej 'a&'

IF a&>127

  SUB a&,256

ENDIF

<txt>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:

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

<txt>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...

<sr>Trochë teorii

<txt>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:

<l>z=x+jx

<txt>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:

<l>SQR(n)

<txt>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:

<l>1..(n/2)+1

<txt>a zignorowaê pozostaîe, czyli w zakresie:

<l>(n/2)+2..n

<txt>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:

<l>a=SQR((x^2)+(jx^2))

<txt>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).