Algorytm motylkowy FFT

[piotr.paw...@gmail.com]

Witam,

Bardzo proszę rozpisanie algorytmu FFT na poszczególne kroki dla 8-punktowej transformaty. Chciałbym w celach edukacyjnych policzyć taka transformatę.
Niestety pomimo zapoznania się z algorytmem motylkowym nie udało mi się samodzielnie takiej transformaty policzyć.
W dokumencie:
https://eti.pg.edu.pl/documents/176593/26763380/Wykl_AlgorOblicz_10.pdf

na stronie ósmej opisano grafem algorytm FFT z podziałem w dziedzinie częstotliwości DIF (ang. decimation in frequency ).
Wszystkie działania są realizowane na liczbach zespolonych dla jednego "motylka" należy wykonać jedno mnożenie, jedno dodawanie i jedno odejmowanie.
Wiem jak policzyć tak zwane współczynniki obrotu. Dla transformaty ośmiopunktowej mamy cztery współczynniki W.

W0 = COS(2pi/1)-jSIN(2pi/1)
W1 = COS(2pi/2)-jSIN(2pi/2)
W2 = COS(2pi/4)-jSIN(2pi/4)
W3 = COS(2pi/8)-jSIN(2pi/8)

Jakie konkretnie działania matematyczne opisuje ten graf? Co oznaczają węzły (miejsca w których spotykają się trzy proste) tego grafu? Co oznacza (-1) przy niektórych węzłach?
Po prostu co do szczegółu nie wiem jak przeprowadzić te obliczenia.

Wiem że są gotowe biblioteki liczące FFT. Ja chciałbym w celu edukacyjnym policzyć taką transformatę na pieszo.

Bardzo proszę o pomoc w postaci konkretnego równania/ń jakie należy policzyć dla grafu opisującego wyżej wymieniony algorytm motylkowy ( proszę opisać przynajmniej początek czyli pierwszy etap).

Za pomoc góry dziękuję.

Piotr Olczyk

[J.F.]

Użytkownik piotr.pawel.olczyk napisał w wiadomości grup
dyskusyjnych:422d4014-e006-45b6...@googlegroups.com...

>Bardzo proszę rozpisanie algorytmu FFT na poszczególne kroki dla
>8-punktowej transformaty. Chciałbym w celach edukacyjnych policzyć
>taka transformatę.
>Niestety pomimo zapoznania się z algorytmem motylkowym nie udało mi
>się samodzielnie takiej transformaty policzyć.
>W dokumencie:
>https://eti.pg.edu.pl/documents/176593/26763380/Wykl_AlgorOblicz_10.pdf

>na stronie ósmej opisano grafem algorytm FFT z podziałem w dziedzinie
>częstotliwości DIF (ang. decimation in frequency ).
>Wszystkie działania są realizowane na liczbach zespolonych dla
>jednego "motylka" należy wykonać jedno mnożenie, jedno dodawanie i
>jedno odejmowanie.
>Wiem jak policzyć tak zwane współczynniki obrotu. Dla transformaty
>ośmiopunktowej mamy cztery współczynniki W.

>W0 = COS(2pi/1)-jSIN(2pi/1)
>W1 = COS(2pi/2)-jSIN(2pi/2)
>W2 = COS(2pi/4)-jSIN(2pi/4)
>W3 = COS(2pi/8)-jSIN(2pi/8)

wygodniej sie poslugiwac e^(-j*...)

P.S. nadmien ze elektronicza/elektryczna - inne branze uzywaja "i"
jako jednostki urojonej.

>Jakie konkretnie działania matematyczne opisuje ten graf?
>Co oznaczają węzły (miejsca w których spotykają się trzy proste) tego
>grafu?

Jesli dobrze widze - to mnozymy przez wspolczynnik przy strzalce i
dodajemy

>Co oznacza (-1) przy niektórych węzłach?

Mnozenie przez -1.

spojrz na wzor 9.1.i i 9.1.j - w jednym masz sume wyrazow, w drugim
roznice.
-1 w motylku to odzwierciedla ..

>Po prostu co do szczegółu nie wiem jak przeprowadzić te obliczenia.

rys 9.1.4 wyglada mi prawidlowo.

Czyli majac macierz na elementy x[i] (zespolone) liczysz

temp=x[0]+x[4]
x[4] = x[0] * W^0 -x[4]
x[0] = temp

i tak 4 razy zwiekszajac potege przy W.

>Bardzo proszę o pomoc w postaci konkretnego równania/ń jakie należy
>policzyć dla grafu opisującego wyżej wymieniony algorytm motylkowy
>( proszę opisać przynajmniej początek czyli pierwszy etap).

Hi hi - jak sobie to rozpiszesz, to otrzymasz zwykle DFT a nie FFT.
To wlasnie trzeba robic motylkami.


P.S.
A potem spojrz na ostatnia strone, dwa pierwsze wzory

X(k) = C(k) + W_k_8* D(k)
X(k+4) = C(k) - W_k_8*D(k)

Tu czynnik W_k_8*D(k) obliczamy raz, a wykorzystujemy dwa razy.
Zakladajac ze jego obliczenie jest kosztowne, a dodawanie i
odejmowanie "tanie" - tu tez sie kryje zysk.

J.

[piotr.paw...@gmail.com]

> spojrz na wzor 9.1.i i 9.1.j - w jednym masz sume wyrazow, w drugim
> roznice.
> -1 w motylku to odzwierciedla ..
>
> >Po prostu co do szczegółu nie wiem jak przeprowadzić te obliczenia.
>
> rys 9.1.4 wyglada mi prawidlowo.
>
> Czyli majac macierz na elementy x[i] (zespolone) liczysz
>
> temp=x[0]+x[4]
> x[4] = x[0] * W^0 -x[4]
> x[0] = temp
>
> i tak 4 razy zwiekszajac potege przy W.

czyli po Etapie 1 na wyjściu

x[0] = x[0] + x[4]

x[4] = x[0]*W^0 - x[4]

Czy na pewno?
A jak to się ma do podanego wzoru 9.1.o tam współczynnik obrotu W jest przy wszystkich składowych.

> P.S.
> A potem spojrz na ostatnia strone, dwa pierwsze wzory
>
> X(k) = C(k) + W_k_8* D(k)
> X(k+4) = C(k) - W_k_8*D(k)
>
> Tu czynnik W_k_8*D(k) obliczamy raz, a wykorzystujemy dwa razy.
> Zakladajac ze jego obliczenie jest kosztowne, a dodawanie i
> odejmowanie "tanie" - tu tez sie kryje zysk.
>
> J.

To dotyczy algorytmu FFT z podziałem w dziedzinie czasu DIT czyli innego podejścia, bo jak wiadomo jest:

algorytm FFT z podziałem w dziedzinie częstotliwości DIF

algorytm FFT z podziałem w dziedzinie czasu DIT

[J.F.]

Dnia Mon, 7 Oct 2019 12:35:36 -0700 (PDT),
piotr.paw...@gmail.com napisał(a):

>> spojrz na wzor 9.1.i i 9.1.j - w jednym masz sume wyrazow, w drugim
>> roznice.
>> -1 w motylku to odzwierciedla ..
>>
>>>Po prostu co do szczegółu nie wiem jak przeprowadzić te obliczenia.
>>
>> rys 9.1.4 wyglada mi prawidlowo.
>>
>> Czyli majac macierz na elementy x[i] (zespolone) liczysz
>>
>> temp=x[0]+x[4]
>> x[4] = x[0] * W^0 -x[4]
>> x[0] = temp
>>
>> i tak 4 razy zwiekszajac potege przy W.
>
> czyli po Etapie 1 na wyjściu
>
> x[0] = x[0] + x[4]
> x[4] = x[0]*W^0 - x[4]
>
> Czy na pewno?
> A jak to się ma do podanego wzoru 9.1.o tam współczynnik obrotu W jest przy wszystkich składowych.

Masz racje. Cos mi to smierdzialo, bo powinny byc dwa czynniki
mnozone, ale nie moj przyklad, spojrzalem sobie na rysunek, zaczelo mi
sie zgadzac ... widac zle spojrzalem

Czyli
x[4] = (x[0] - x[4] )*W^0

>> P.S.
>> A potem spojrz na ostatnia strone, dwa pierwsze wzory
>>
>> X(k) = C(k) + W_k_8* D(k)
>> X(k+4) = C(k) - W_k_8*D(k)
>>
>> Tu czynnik W_k_8*D(k) obliczamy raz, a wykorzystujemy dwa razy.
>> Zakladajac ze jego obliczenie jest kosztowne, a dodawanie i
>> odejmowanie "tanie" - tu tez sie kryje zysk.
>>

> To dotyczy algorytmu FFT z podziałem w dziedzinie czasu DIT czyli innego podejścia, bo jak wiadomo jest:
> algorytm FFT z podziałem w dziedzinie częstotliwości DIF
> algorytm FFT z podziałem w dziedzinie czasu DIT

Szumne nazwy, a ciagle to samo :-)
Tylko podejscie do podzialu odwrotne ...

J.

[piotr.paw...@gmail.com]

Chcę dopytać o współczynniki obrotu. Jak policzyć wykładnik potęgi tego współczynnika na danym etapie?
Dlaczego na przykład w etapie pierwszym wykładnik przy współczynniku W to 0,1,2,3 na etapie drugim wykładnik współczynnika W to 0 i 2? Dlaczego są w etapie drugim dwa współczynniki a na etapie 1 cztery?
Jak wyglądał by schemat dla FFT 16-punktowego ile byłoby współczynników i jakie by były przy nich wykładniki?

[J.F.]

Użytkownik piotr.pawel.olczyk napisał w wiadomości grup

dyskusyjnych:4618e16e-c2b5-4bc7...@googlegroups.com...

W dniu wtorek, 8 października 2019 00:41:15 UTC+2 użytkownik J.F.
napisał:

[...]

>Chcę dopytać o współczynniki obrotu. Jak policzyć wykładnik potęgi
>tego współczynnika na danym etapie?
>Dlaczego na przykład w etapie pierwszym wykładnik przy współczynniku
>W to 0,1,2,3 na etapie drugim wykładnik współczynnika W to 0 i 2?
>Dlaczego są w etapie drugim dwa współczynniki a na etapie 1 cztery?
>Jak wyglądał by schemat dla FFT 16-punktowego ile byłoby
>współczynników i jakie by były przy nich wykładniki?

Na pierwszym etapie stosujesz W, moze lepiej nazwijmy go
W0=exp(-j*2*pi/N), gdzie N to ilosc punktow transformaty, czyli np 16.
I potegi od 0 do polowy, tzn N/2-1, czyli do 7.

(Bo drugą połowę potęg załatwia to -1 - czyli W^(N/2) = exp(-j*pi)=-1

Na kolejnym etapie ... to juz juz kwestia nazewnictwa - tam sa dwie
transformaty o połowie punktow,
wiec nalezaloby uzyc wspolczynnika W1=exp(-j*2*pi/(N/2)) co jest
rowne W0^2.

I mozna mowic, ze to jest W1 do kolejnych poteg 0, 1, 2, 3, albo W0
do poteg parzystych : 0, 2, 4, 6

Na kolejnym bedzie W2=W1^2=W0^4 i kolejne potegi W2 0, 1,
albo jak kto woli W0 do poteg wielokrotnosci 4, czyli W0^0 i W0^4


To sie ladnie sklada informatycznie - w etapie uzywasz wspolczynnika
W, ktory miedzy etapami podnosisz do kwadratu.
A w danym etapie uzywasz jego kolejnych poteg - tylko ze z kazdym
etapem spada nam ilosc tych poteg o polowe,
za to rosnie ilosc motylkow z jedna potegą.

J.

[piotr.paw...@gmail.com]

Czyli dla 16-punktowej FFT w podejściu (z podziałem w dziedzinie częstotliwości DIF), współczynniki obrotu W^x są następujące:

Etap1
W^0,W^1,W^2,W^3,W^4,W^5,W^6,W^7

Etap2
W^0,W^2,W^4,W^8

Etap3
W^0,W^8

Etap4
-

Czy to się zgadza?

[piotr.paw...@gmail.com]

Tutaj znalazłem dość dokładny opis algorytmu FFT z podzałem w dziedzinie częstotliwości DIF.

https://web.archive.org/web/20171114115729/http://en.dsplib.org/content/fft_dec_in_freq.html

[J.F.]

Użytkownik piotr.pawel.olczyk napisał w wiadomości grup

dyskusyjnych:f888c1df-dbdb-41d1...@googlegroups.com...

W dniu wtorek, 8 października 2019 17:16:34 UTC+2 użytkownik J.F.
napisał:
> Użytkownik piotr.pawel.olczyk napisał w wiadomości grup

> >Chcę dopytać o współczynniki obrotu. Jak policzyć wykładnik potęgi
> >tego współczynnika na danym etapie?

[...]

>Czyli dla 16-punktowej FFT w podejściu (z podziałem w dziedzinie
>częstotliwości DIF), współczynniki obrotu W^x są następujące:
>Etap1
>W^0,W^1,W^2,W^3,W^4,W^5,W^6,W^7

>Etap2
>W^0,W^2,W^4,W^8

W^0,W^2,W^4, W^6

>Etap3
>W^0,W^8

W^0, W^4

>Etap4
>-

J.

[piotr.paw...@gmail.com]

>W dniu poniedziałek, 7 października 2019 18:13:38 UTC+2 użytkownik J.F. napisał
>
> wygodniej sie poslugiwac e^(-j*...)

Zapytam o sposób liczenia współczynnika obrotu W z postaci wykładniczej piszesz że wygodniej jest się posłużyć postacią wykładniczą.
Obliczam współczynniki przykładowo dla 8-punktowej FFT, etap 1

Licząc z postaci trygonometrycznej wychodzi tak:

W0 = cos(2pi/1) - sin(2pi/1)= 1
W1 = cos(2pi/2) - sin(2pi/2)= -1
W2 = cos(2pi/4) - sin(2pi/4)= -i
W4 = cos(2pi/8) - sin(2pi/8)= 0.707 - i0.707

Sugerujesz aby posługiwać się postacią wykładniczą, jak policzyć wykładniki z postaci wykładniczej e(-j*…)?

[piotr.paw...@gmail.com]

Przepraszam w ostatnim współczynniku jest pomyłka w indeksie. Powinno być
W3 = cos(2pi/8) - sin(2pi/8)= 0.707 - i0.707

[J.F.]

Użytkownik piotr.pawel.olczyk napisał w wiadomości grup

dyskusyjnych:8988e2e2-a635-4822...@googlegroups.com...

na papierze wygodniej e^-j* ... :-)

W programie raz wyliczasz W = cos(2pi/N) - j*sin(2pi/N)
a potem tylko wyliczasz kolejne potegi za pomoca mnozenia
(zespolonego), a miedzy etapami jak widziales - wystarczy W^2 zrobic.

No chyba, ze uzywany jezyk ma zgrabna funkcje exp(..) zespoloną.

Ewentualnie tablicujesz sobie te wspolczynniki - co ma sens, jak
transformaty krotkie, pamiec nie boli, a transformat bedzie duzo.

J.

[piotr.paw...@gmail.com]

>W dniu środa, 9 października 2019 16:20:33 UTC+2 użytkownik J.F. napisał:
>
> >Obliczam współczynniki przykładowo dla 8-punktowej FFT, etap 1
> >Licząc z postaci trygonometrycznej wychodzi tak:
>
> >W0 = cos(2pi/1) - sin(2pi/1)= 1
> >W1 = cos(2pi/2) - sin(2pi/2)= -1
> >W2 = cos(2pi/4) - sin(2pi/4)= -i
> >W4 = cos(2pi/8) - sin(2pi/8)= 0.707 - i0.707
>

>W programie raz wyliczasz W = cos(2pi/N) - j*sin(2pi/N)

Czy możesz potwierdzić że dla transformaty 8-punktowej prawidłowo policzyłem współczynniki obrotu W(N) dla pierwszego etapu? Zwróć proszę szczególnie uwagę na N.
Podajesz wzór W(N) = cos(2pi/N) - j*sin(2pi/N) czy on jest poprawny?
Czy nie powinno być W(N) = cos(2pi/2^N) - j*sin(2pi/2^N)

[J.F.]

Dnia Fri, 11 Oct 2019 12:04:46 -0700 (PDT),
piotr.paw...@gmail.com napisał(a):

U mnie N jest iloscia punktow transformaty. czyli np 8, 16, 32,
a nie 3, 4, 5.

czyli np dla 16 pkt zaczynasz od
W4 = cos(2pi/16) - j*sin(2pi/16)

w kolejnym etapie potrzebny bedzie W3=W4^2,
w kolejnym W2=W3^2
itd.

J.