Algorytmy: FFT

Algorytmy: FFT

  Cooley, Lewis i Welch (1969) przedstawili bardzo prostą i dość efektywną realizacje algorytmu FFT do wyliczenia DFT ( The Fast Fourier Transform and Its Applications, IEEE Transactions on Education, vol.12, no.1,pp.27-34 )

Procedurę tą przytoczono też w dziele „Theory and Application of Digital Signal Processing ”, L.Rabiner – Bell Laboratories, B.Gold – MIT Lincol Laboratory, Prentice Hall , 1975.

Adaptowanie tekstu Fortranu na inny język programowania jest proste. Zaleta Fortranu jest dostępność typu zespolonego COMPLEX. W innych językach operacje na liczbach zespolonych musimy wykonywać na piechote. Zamiast zespolonego wektora A() dajemy dwa wektory Ar() i Ai(), jak r-real oraz i-image.

Bardzo ważne jest przetestowanie algorytmu aby nabrac pewności ze daje zawsze dobre wyniki i jest przewidywalny ( czas dzialania w funkcji rozmiaru danych wejsciowych ) oraz stabilny. Do tego celu potrzebne są dane testowe. Możemy je sporzadzić samodzielnie lub lepiej skorzystać z danych testowych udostępnianych przez choćby IBM. Z systemu IBM 360/370 nie ma prostego sposobu przeniesienia plików do mikrokomputera.

Przed przystąpieniem do testu musimy znać wynik. W psychike ludzką wbudowany jest efekt racjonalizowania. I tak wyniki błednie dzialającej procedury zracjonalizujemy jako dobre.

Fortranowska procedura ponizej jest nierekurencyjna i działa „w miejscu” ( in place ) bez rezerwacji dodatkowego miejsca. Ma więc zalety. Kod DO 7 do etykiety 7 służy do inversyjno - dwójkowego przestawienia elementów wejściowego wektora A(N). Zauważmy ze brak instrukcji contunue obniża czytelność kodu.

Rysunki w powołanym dziele „Theory and Application...” pokazują efekt przestawień wektora. Aby zobaczyć czy przestawianie działa jak należy wpisujemy do A() kolejne liczby naturalne i drukujemy A() po etykiecie 7 nie wykonujac reszty procedury. Dla wartośći N=8 i N=32 przestawianie działa jak należy. W praktyce wektory są dużo większe.

Dla M=10, N wynosi 1024 a ilośc przestawień wynosi 496. Za instrukcją A(I) = T umieścmy wstepnie wyzerowaną inkrementowa zmienną licznika aby policzyć ilosc przestawień. Wykonywana ilość przestawień jest poprawna.

Dla M=10, N=1024 ilość „motylków” wynosi ( M/2 ) x N = 5120. Licznik umieszczamy za instrukcja „10 A(I) = A(I)+T”. Wynik testu jest znów poprawny.

   SUBROUTINE FFT(A,M,N)

   COMPLEX A(N),U,W,T

   N = 2**M

   NV2 = N/2

   NM1 = N-1

   J=1

   DO 7 I=1,NM1

   IF(I.GE.J) GO TO 5

   T = A(J)

   A(J) = A(I)

   A(I) = T

5 K=NV2

6 IF(K.GE.J) GO TO 7

   J = J-K

   K=K/2

   GO TO 6

7 J = J+K

   PI = 3.141592653589793

   DO 20 L=1,M

   LE = 2**L

   LE1 = LE/2

   U = (1.0,0.)

   W = CMPLX(COS(PI/LE1),SIN(PI/LE1))

   DO 20 J=1,LE1

   DO 10 I=J,N,LE

   IP = I+LE1

   T = A(IP)*U

   A(IP) = A(I)-T

10 A(I) = A(I)+T

20 U=U*W

   RETURN

   END

Potrzebna jest mala dygresja. Koncern IBM bardzo szybko dostrzegł że nakład pracy na wykonanie programu może być bardzo duży. Stąd pomysł na prymitywny translator Fortran ( = Formula Translation ), który pierwszy powstał już w połowie lat pięćdziesiątych. Konstruktorzy mieli świadomość ze kod generowany przez translator musi być efektywny to znaczy i w miare mały i szybki bowiem inaczej uzytkownicy go nie zaakceptuja ! 

Prawie każdy translator Fotranu optymalizuje kod. Instrukcja K=K/2 zostanie pewnie przetłumaczona jako ASR czyli Artmetic Shift Right czyli artmetyczne przesunięcie w prawo, równoważna dzieleniu przez 2 bowiem zmienna K jest calkowitoliczbowa. Dzielenie liczb stałoprzecinkowych bywa zdumiewajaco wolne.

Fortran nadaje się wyłącznie do obliczeń naukowo-technicznych. Z reguły program wynikowy jest istotnie szybszy niż z innych kompilatorów.

W Fortranie macierze pamiętane sa kolumnami a nie wierszami jak w pozostałych językach programowania. Bywa to zaleta i były powody aby zastosować taką organizacje.

Irytująca instrukcja „GO TO” sluzy do poruszania się po programie. Wydaje się że taka instrukcja we wpolczesnych jezykach powinna być zabroniona. Logika programu Fortranowskiego jest jednak bardzo trudna do zrozumienia i to jest poważna wada. Pożniejsze wersje Fortranu maja już instrukcje IF..THEN … ELSE...END IF.

Zewnętrzne instrukcje z sekwencji

6 IF(K.GE.J) GO TO 7

J = J-K

K=K/2

GO TO 6

to przecież pętla While. Wiedząc jaka sekwencja instrukcji fortranowskich odpowiada instrukcjom strukturalnym można sobie w myslach czy na papierze przetłumaczyc program i probować go analizować.

Kontynujmy testowanie procedury. Niech pierwszy element A będzie liczbą rzeczywistą a pozostale zerowe. Transformata Fouriera impulsu Diraca jest rzeczywista i wszystkie elementy stransformowanego wektora A powinny mieć ta sama wartośc co wstawiliśmy do A. I tak jest . Jesli tą wartosc wstawimy na kolejne miejsce to moduły dalej bedą mialy ta wartosc ale pojawią się przesunięcia fazy.

Przeksztalcenie FFT jest liniowe i przeksztalcenie z wazonej sumy dwóch dowolnych wektorów jest ważoną sumą przekształceń wektorow.

W procedurze U nie jest wektorem i wyliczane jest na bieżąco. Autorzy podają ze spowalnia to program o mniej niż 15%. Wyliczenia U są bardzo mało dokładne i to jest spora wada tej procedury. Koztem niewielkiej rozbudowy można wade usunąc.

W praktyce wektor A(N) jest rzeczywisty a nie zepolony. Aby zaoszczędzic czasu i miejsca w pamięci stosuje się odpowiednie prosciutkie podstawienia - transformacje.

Są one podane na stronie 225 w „Przetwarzanie sygnałow metodami analogowymi i cyfrowymi, K.G.Beauchamp, WNT 1978. Podano tam też ich ciekawe wyprowadzenie.

W najprostszym razie transformacje FFT stosuje się do wyliczania widma sygnałów. Z uwagi na zjawisko przecieku spektralnego ( leakage ) wpierw zebrane próbki sygnału mnożymy przez odpowiednie okno: Kaiser, Hamming, Hann, Blackmann-Harris, Gauss, Trójkątne - Bartlett, Trapezoidalne, Welch, itd. Wyliczenie okien jest trywialne. Różne cechy okna mają kompromisowe własności. Rabiner i Gold dają wystarczająco dużo informacji pozwalajacych wybrać własciwe do zastosowania okno. Transformata FFT już ma bardzo szerokie zastosowania profesjonalne ( radary, echosondy, medycyna, geofizyka, badania kosmosu, badania maszyn ) i w miare popularyzacji komputerow obszar będzie wzrastał.

Komputerek Spectrum nadaje się tylko do transformacji malych wektorów ( zalezy od tego ile mamy ciepliwosci w oczekiwaniu na wyniki ) ale wyniki są jak najbardziej poprawne.