Transformacja fouriera - prośba o wytłumaczenie

Question

---

Cześć, chciałem napisać program, który przyjmowałby pewną funkcję i zwracałby jej widmo, jednakże jak szukałem, to wszelkie biblioteki(np. ta http://www.aforgenet.com/framework/docs/html/e1932533-0745-940d-c998-140e8ee23927.htm) bazują albo na tablicy liczb zespolonych, albo na wektorach. Chciałbym prosić o wytłumaczenie, najlepiej jak dla idioty, jak to działa, w sensie, czym są te dane(w tym wypadku liczby zespolone i wektory. Znaczy, wiem czym są liczby zespolone i wektory, ale czym są te konkretne, używane w tych funkcjach) i jak je "wyciągnąć" z "zwykłej" funkcji, np. sinusoidy

Comments

Możesz tego dokonać niemal w każdym języku, ale chyba najłatwiej będzie w Pythonie (z użyciem np. scipy i numpy), FORTRAN-ie lub jeśli masz dostęp, w MATLAB-ie (mogą być jego odpowiedniki FOSS, szczególnie Octave).

Rzuć okiem na stronę teoretyczną:

https://www.youtube.com/watch?v=spUNpyF58BY
But what is the Fourier Transform? A visual introduction. - YouTube

https://www.youtube.com/watch?v=A9A_EyMF7JE
Transformata Fouriera - YouTube

https://www.youtube.com/watch?v=8_S6XKK53c4
Czym jest transformata Fouriera? - YouTube

Answer 1

Matematyka zna metody rozbijania dowolnych  funkcji na składowe o ustalonym przebiegu. Używa się tego w obliczeniach komputerowych, gdzie różne funkcje dostępne w komputerze liczy się poprzez odpowiadające im ciągi wielomianowe itp.

Transformata Fouriera to metoda rozbicia dowolnej funkcji na składowe sinusoidalne.

Istnieje analogowa transformata Fouriera, która działa na "funkcjach/przebiegach" i jej wynikiem sa funkcje oraz istnieje cyfrowa transformata Fouriera (DFT), która działa na ciągach liczb zespolonych i takie też wyniki daje.

Jeśli argumentem analogowej transformaty jest funkcja okresowa, wynikiem jest funkcja która jest różna od zera jedynie w pewnych ustalonych punktach.

Jeśli chcesz zamienić jakąś "analogową" funkcję matetyczną na ciąg liczb, musisz przeprowadzić coś co się nazywa "próbkowaniem", czyli obliczyć kolejne wartości funkcji w równoodległych punktach zmiennej x. Tutaj istotne jest twierdzenie o próbkowaniu, które określa minimalną "gęstość próbkowania", przy którym ciąg liczb jest równoważny z funkcją "analogową". Możesz np. wziąć sinus i probkować go np. 20 (duży zapas) razy w ramach każdego okresu (0-2*PI) - takiej jednej "falki", wrzucić te danę na DFT i dostaniesz 1 pik - ale z przyległościami wynikającymi ze skończonej długości ciągu.

FFT to jest pewien algorytm liczenia DTF.

Comments

Ok, coś mi tam zaczyna świtać, jednakże co będzie częścią urojoną, a co rzeczywistą? W sensie, to próbkowanie mam przeprowadzić w ten sposób, że biorę liczbę zespoloną i np. częścią rzeczywistą jest punkt X, a próbką (czyli wynikiem na osi Y z funkcji) jest część urojona?

---

Ostatnio musiałem odświeżyć moje wiadomości z DSP , a nawet nauczyć się nowych rzeczy dotyczących sygnałów zespolonych. Ale nic, liczby rzeczywiste są podzbiorem liczb zespolowych, dane wpisujesz jako Re, składową Im dajesz zero. Jednak wynik będzie zespolony, będzie miał niezerowe obie składowe.

---

Dobra, wykombinowałem coś takiego:

double func(double t)
            {
                return 10 * Math.Sin(200 * Math.PI * t) + 10 * Math.Sin(800 * Math.PI * t) + 10 * Math.Sin(2000 * Math.PI * t);      
            }

            Complex[] data = new Complex[64];

            for (int i = 0; i < 64; ++i)
            {
                data[i] = new Complex(func(0.02 * (i / 64.0)), 0);
            }

            FourierTransform.DFT(data, FourierTransform.Direction.Forward);

            for(int i = 0; i < 64; ++i)
            {
                Console.Write(data[i] + "\n");
            }

I otrzymuję jakiś tam wynik, teraz ostatnie pytanie, jak go rozczytać? W sensie, jak mógłbym go przełożyć na wykres? Po prostu uznać jedną wartość za oś X, a drugą za oś Y? Jak chcę zrobić coś takiego, to wychodzi bardzo szalony wykres, więc zgaduję, że nie

---

Na osi poziomej masz "częstotliwości". Zrób dwie osie pionowe lub dwa wykresy. Na jednym nanieś moduły poszczególnych liczba zespolonych (re^2 + im^2), na drugim fazę (atan2(re, im)).

Pierwszy wykres poda jak "dużo" składowej sinusoidalnej o danej częstotliwości należy wziąć do sumy by odtworzyć pierwotny sygnał/dane. Drugi wykres to informacja o ile należy przesunąć poszczególne sinusy.

Jeśli, jak w przyładowym kodzie, dane wejściowe są rzeczywiste, wtedy pierwszy wykres powinien wyjść parzysty, a drugi nieparzysty (takie symetrie).

W sumie realne dane zawiera połowa wyniku czyli te w zakresie 0-31.

Dałeś trochę małe tablice, by dostać ładny wykres daj co najmniej 512, będzie pasować do rozdzielczości ekranu.

---

No, wygląda to już lepiej, tylko nie wiem dlaczego wynik jest zdublowany i odbity

Martwi mnie też skala Y. Nie wiem czy się gdzieś nie kopsnąłem i nie ma przecieków

Kod wygląda tak w chwili obecnej:

 

const double period = 2 * Math.PI;
            const double intervals = 3;
            const int howMany = 512;

            double func(double t)
            {
                //return 10 * Math.Sin(200 * Math.PI * t) + 10 * Math.Sin(800 * Math.PI * t) + 10 * Math.Sin(2000 * Math.PI * t);
                return Math.Sin(t * Math.PI * 2);
            }

            Complex[] data = new Complex[howMany];

            for (int i = 0; i < howMany; ++i)
            {
                data[i] = new Complex(func(period * intervals * (i / Convert.ToDouble(howMany))), 0);
            }

            FourierTransform.DFT(data, FourierTransform.Direction.Forward);

            for (int i = 0; i < howMany; ++i)
            {             
                yValues.Add(data[i].Re * data[i].Re + data[i].Im * data[i].Im);               
                xValues.Add(period * intervals * (i / Convert.ToDouble(howMany)));            
            }

---

"Odbicie" to właśnie ta parzysta symetria. Dla sygnałów rzeczywistych tylko połowa wyniku jest znacząca, reszta to "odbicie". Skala Y wymaga normalizacji, ale nie pamiętam jaki jest współczynnik - może ta funkcja już to robi? Tak czy inaczej powinien być tam jeszcze pierwiastek - jak w tw. Pitagorasa sqrt(x^2+y^2), który wstępnie pominałem. Cześto się jeszcze stosuje logarytm by uzyskać "ludzką" skale.

Rozumiem, że drugi wykres pokazuje już wszystkie trzy składowe sinusoidalne (200, 800 i 2000)?

---

Tak, drugi wykres pokazuje 3 składowe. A ta skala właśnie mnie martwi, ponieważ na przykładzie, który otrzymałem od profesora, w tym wzorze z 3 składowymi, to transformacja pokazuje właśnie to, tylko, że do 10 a nie do 25. Za to jak patrzyłem na przykłady z sinusem, to co przykład to inna skala, więc może tak to powinno być? (xD)

---

Wskazana amplituda zależy często od doboru parametrów. Poszczególne "prążki" wyniku odpowiadają pewnym częstotliwościom na wejściu. Jeśli jednak sygnał wejściowy nie trafia dokładnie w żadną z częstotliwości wyjściowych to zostanie podzielony na dwa kolejne prążki, jednak ich amplituda będzie odpowiednio mniejsza. Częstotliwości prążków są to kelejne wielokrotności częstotliwości próbkowania podzielonej przez rozmiar transformaty (wielkość tablicy). Może się trafić, że przy jednych parametrach uda się trafić i wyjdzie jeden prążek o dużej amplitudzie, a w innych dostanie się dwa mniejsze.

---

Dobra, udało mi się to wszystko rozwiązać jakiś czas temu i zapomniałem zedytować komentarza. Tak więc, dzięki za pomoc, bo bez ciebie bym tego nie zapisał A dla potomnych jeszcze zostawię ten wzór na oś X:

(double)i * sr / howMany

gdzie:
i, jest iteracją pętli(tej samej co iteruje oś y)
sr jest częstotliwością próbkowania w hercach
howMany to ilość próbek I iterację należy zacząć od 1 a nie od 0, inaczej nastąpi przesunięcie

I cały kod ostatecznie wygląda tak:

 

const double period = Math.PI * 2;
            const double intervals = 3;
            const int howMany = 512;

            double func(double t)
            {
                //return 10 * Math.Sin(200 * Math.PI * t) + 10 * Math.Sin(800 * Math.PI * t) + 10 * Math.Sin(2000 * Math.PI * t);
                /*
                if(Math.Sin(t * Math.PI * 2) > 0)
                {
                    return 1;
                }
                else
                {
                    return -1;
                }
                */
                return Math.Sin(t * Math.PI * 2);
            }

            Complex[] data = new Complex[howMany];

            for (int i = 0; i < howMany; ++i)
            {
                data[i] = new Complex(func(period * intervals * (i / Convert.ToDouble(howMany))), 0);
            }

            FourierTransform.DFT(data, FourierTransform.Direction.Forward);

            double[] bins = new double[howMany];

            double maxBin = 0;

            for (int i = 0; i < howMany; ++i)
            {
                bins[i] = Math.Sqrt(data[i].Re * data[i].Re + data[i].Im * data[i].Im);

                if(bins[i] > maxBin)
                {
                    maxBin = bins[i];
              }
            }

            double fs = 0.5 / (maxBin * 2);
            double sr = howMany / (period * intervals);

            for (int i = 1; i < howMany / 2; ++i)
            {
                yValues.Add(bins[i] * 2);
                xValues.Add((double)i * sr / howMany);
            }