Cyfrowe przetwarzanie sygnałów

To nie jest miejsce tylko dla początkujących, wszyscy jesteśmy w czymś początkujący i wymieniamy się doświadczeniami.
Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » wtorek 11 sty 2022, 00:06

Jean Baptiste Joseph Fourier (ur. 21 marca 1768
w Auxerre, zm. 16 maja 1830 r. w Paryżu) – francuski
matematyk, który jest sprawcą
całego zamieszania.


Wstęp

sig01.png


W dobie cyfrowego dźwięku, cyfrowych obrazów (w sensie obrazów nieruchomych → zdjęć oraz ruchomych → filmy) cyfrowej automatyki warto wypracować sobie elementarne pojęcie dotyczące „filozofii” cyfrowego przetwarzania sygnałów. W ogólnym przypadku cała ta gama zastosowań i rozwiązań określana jest jako DSP (od ang. Digital Signal Processing – cyfrowe przetwarzania sygnałów). To bardzo nowoczesna idea, która niejako przebojem wbiła się we współczesność. Na razie nie będziemy wybierać się w kosmos. Nie wyobrażam sobie, by świeżo wypuszczony w kosmiczną nieskończoność Teleskop James'a Webb'a nie posługiwał się ta technologią. Ten sprzęt, jak sądzę, służy do przetwarzania i badania sygnałów „bardzo anemicznych”. Dodając do tego szum generowany przez tło, potrafi skutecznie utrudnić życie. Jak uzyskać dobry rezultat, gdzie niechciane przeszkody są często większe od informacji użytecznej. Jednak pomysłowość ludzka nie zna granic. W historii wymyślono różne filtry pozwalające na znaczącą poprawę procesu badawczo-poznawczego. No ale „nie od razu Kraków zbudowano”. Metodologią drobnych kroczków idziemy do przodu. Na razie kosmos jest poza zasięgiem, ale nie znaczy, że tak zawsze będzie. Per aspera ad astra → nikt nie twierdzi, że będzie łatwo, ale również nikt nie twierdzi, że się nie da.
W poniższych rozważaniach postaram się przedstawić w prosty sposób (mam taką nadzieję) własne doświadczenia w temacie. Nie jest wykluczone, że dla niektórych Czytelników poruszane tematy nie wnoszą nic nowego, dlatego proszę o wyrozumiałość.
Jedną z istotnych trudności ze zrozumieniem „filozofii” przetwarzania sygnałów jest zaawansowany aparat matematyczny związany z tematem. W ogólnym przypadku są to złożone rozważania matematyczne realizowane z przestrzeni liczb zespolonych. Samo pojęcie liczb zespolonych dla wielu może być tajemnicze a u innych wywołuje „opór” przed dalszym zagłębianiem się w tematykę. Jednak możliwe jest „zbliżenie się” do tematu przetwarzania sygnałów bez zagłębiania się z tematykę liczb zespolonych (przynajmniej na początek). Niestety pewne rozważania matematyczne „z wyższej półki” nadal pozostają, jednak postaram się „nie kaleczyć” przesadną matematyką.
Rozpatrzmy przykładowo zwykły wzmacniacz audio. Na jego wejście doprowadzony jest jakiś sygnał dźwiękowy. Sam wzmacniacz z punktu widzenia elektroniki jest jakimś układem, który przetwarza ten sygnał. W sumie żadna wielka filozofia. Z takiego radyjka sączy się jakaś fajna muza. Wciskamy jakiś magiczny przycisk w radyjku i z całego sygnału audio wycięło wokal, została tylko orkiestra. Czy jest to możliwe? Może tak, może nie. Ja nie wiem, ale to nie powód by tego nie spróbować.
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » wtorek 11 sty 2022, 00:45

Sygnał analogowy

Sygnał analogowy – sygnał, który może przyjmować dowolną wartość z ciągłego przedziału (nieskończonego lub ograniczonego zakresem zmienności). Jego wartości mogą zostać określone w dowolnej chwili czasu. Zwyczajowo ten sygnał określa się terminem: funkcja czasu.
sig02.png

Jeżeli mamy dowolny przebieg sygnału f(t) (jako funkcja czasu: ilustracja 1), w sensie funkcji matematycznej z określonym okresem (zwyczajowo oznaczanym literką T), to tą funkcję można aproksymować szeregiem trygonometrycznym (aproksymować, czyli zastąpić daną funkcję zbiorem innych, które wystarczająco dokładnie odwzorują oryginalną funkcję, a które dają się później łatwo obrabiać matematycznie, szczególnie użyteczne są tu funkcje trygonometryczne), ilustracja 2.
sig03.png

Oczywiście rodzi się tu pewna wątpliwość (z tą powtarzalnością): jest to w dużej mierze świat wyidealizowany natomiast rzeczywistość praktycznie zawsze odbiega od tego wyidealizowanego obrazu, niemniej można próbować zbliżyć się tej idealności. Po cichu załóżmy, że tak właśnie jest → sygnał powtarza się z określonym okresem.
Szereg ten (trygonometryczny) można przedstawić jako sumę pewnej liczby funkcji cos i sin o odpowiednich amplitudach określonych poprzez wartości współczynników an i bn wraz z określonymi pulsacjami oraz wartości a0 będącej wartością średnią funkcji w przedziale. W sensie matematycznym można to zapisać następująco:
part01_form001.png

Powyższa zależność pozwala aproksymować dowolną funkcję okresową sygnału zestawem odpowiednich funkcji trygonometrycznych. Jednak w tej aproksymacji cała trudność polega na wyznaczeniu odpowiednich współczynników. Oczywiście jest to trudność „pozorna” gdyż „wielcy tego świata” (Fourier) pokazali jak to można zrobić. W powyższym związku najłatwiej jest wyjaśnić znaczenie współczynnika a0, który odzwierciedla wartość średnią funkcji w przedziale okresowym. Wyjaśnia to poniższa ilustracja 3.
sig04.png

Wyobraźmy sobie, że na funkcję naciska w góry prasa, która rozgniecie funkcję do postaci prostokąta. Wszystkie górki wypełnią dolinki i powstanie plaski teren. Ten poziom jaki się ustali to właśnie będzie wartość średnia. Wracając do świata matematyki, należy wyznaczyć taką wartość współczynnika a0, by sygnał miał wartości średnią. Z punktu widzenia matematyki jest to wartość średnia funkcji w zadanym przedziale. Ten przedział to (-π,π) (dlaczego taki, będzie w dalszej części). Matematyka potrafi takie rzeczy policzyć używając do tego jednego z najbardziej genialnych swoich wynalazków: całki. Mamy do obliczenia następująca całkę:
part01_form002.png

W ogólnym ujęciu kolejne poszczególne harmoniczne uzyskuje się w wyniku całkowania funkcji opisującej sygnał w iloczynie z funkcjami trygonometycznymi dla poszczególnych harmonicznych. Dochodzi do „wyfiltrowania” określonych harmonicznych z sygnału. By wyjaśnić mechanizm działania, w oparciu o który obliczane są poszczególne współczynniki (te współczynniki an i bn tworzą razem amplitudę) każdej z harmonicznych należy przyjrzeć się pewnym operacjom matematycznym. Wspomniana amplituda to (oczywiście każda harmoniczna ma swoją własną):
part01_form003.png

Niech rzeczywisty sygnał będzie aproksymowany następującym szeregiem trygonometrycznym:
part01_form004.png

Wprowadzenie pulsacji kołowej ω „synchronizuje” rzeczywisty przebieg funkcji okresowej z funkcją trygonometryczną (ilustracja 4). Po takim zabiegu, z punktu widzenia funkcji trygonometrycznej, funkcja okresowa ma również okres . Podobnie funkcje trygonometryczne sin(ωt) jak i cos(ωt) mają również swój okres wynoszący T. Łatwo to sprawdzić podstawiając do tych zależności parametry przedziału okresowości (funkcja sin rozpina się na krańcach przedziału całym swoim okresem):
part01_form005.png

Identycznie jest z funkcją cos.
Wszystkie rozkminy rachunkowe sprowadzają się do przedziału [-T/2 , T/2].
sig05.png

Dokonane zostanie obliczenie wartości całki oznaczonej w granicach całkowania obejmujących pojedynczy okres iloczynu powyższej funkcji i funkcji trygonometrycznych sin(kωt) i cos(kωt) pełniącymi rolę „filtru” dla sygnału, gdzie odpowiada określonej harmonicznej (jest całkowitoliczbową krotnością harmonicznej podstawowej ω). Mamy do obliczenia:
part01_form006.png

Warto tu zauważyć, że powyższą sumę można obliczyć „po kawałku”. Na początek zostaje obliczona:
part01_form007.png

co oznacza, że powyższa operacja „filtrowania” wartości średniej (reprezentowanej przez współczynnik a0) daje wartość zerową bez względu na numer „wycinanej” harmonicznej reprezentowanej przez parametr .
Kolejnym fragmentem jest:
part01_form008.png

Z punktu widzenia rachunkowego koniecznością staje się rozpatrzenie dwóch wariantów, gdzie k≠n oraz k=n. Pierwszy przypadek w postaci:
part01_form009.png

gdzie jest k≠n, co oznacza „wyfiltrowanie” innej niż k-tej harmonicznej z całego zestawu harmonicznych:
part01_form010.png

gdyż dowolna wielokrotność połowy pełnego okresu funkcji sin jest zerem sin(kπ)=0.
W przypadku gdy k=n, czyli gdy zachodzi „filtrowanie” własnej harmonicznej, obliczenia są następujące:
part01_form011.png

Kolejny fragment całki
part01_form012.png

również wymaga rozpatrzenia dwóch przypadków, gdzie następuje „filtrowanie” innej harmonicznej (k≠n) oraz własnej harmonicznej (k=n). Pierwszy przypadek daje następujący wynik:
part01_form013.png

drugi przypadek:
part01_form014.png

W przypadku wydzielenia własnej harmonicznej wynik jest różny od zera.
Wracając do całej wyjściowej całki, jej wartość jest następująca:
part01_form015.png

lub
part01_form016.png

Jest to kwestia punktu widzenia: czy brany jest pod uwagę okres T funkcji okresowej, czy okres (wynoszący ) funkcji trygonometrycznej → ogólnie bn pomnożone przez połowę długości okresu.
Analogiczne obliczenia wartości całki oznaczonej z iloczynu funkcji czasu i funkcji cos(kωt), gdzie k jest liczbą całkowitą pełniącą funkcję „filtru” dla sygnału dają wartość amplitudy odpowiedniej harmonicznej (ilustracja 5).
sig06.png

part01_form017.png

Gdzie poszczególne elementy to:
part01_form018.png

part01_form019.png

part01_form020.png

part01_form021.png

part01_form022.png

Ostatecznie, skoro
part01_form023.png

part01_form024.png

part01_form025.png

Powyższe rozważania umożliwiają obliczenie wartości odpowiednich współczynników w szeregu trygonometrycznym pozwalającym aproksymować dowolny przebieg okresowej funkcji opisującej sygnał. Ten szereg można zapisać następująco:
part01_form026.png

A poszczególne jego współczynniki wiadomo jak obliczyć.
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
Zegar
User
User
Posty: 316
Rejestracja: wtorek 02 lip 2019, 14:42

Re: Cyfrowe przetwarzanie sygnałów

Postautor: Zegar » wtorek 11 sty 2022, 10:26

gaweł pisze:jednak postaram się „nie kaleczyć” przesadną matematyką.

gaweł pisze:A poszczególne jego współczynniki wiadomo jak obliczyć.

No jasne. Teraz już będzie z górki. :lol:
Przypomniały mi się szkolne lata. Ostatni raz liczyłem całki w ubiegłym wieku. :D
"If A = success, then the formula is A = X + Y + Z.
X is work. Y is play. Z is keep your mouth shut."
A. Einstein

Awatar użytkownika
GrumpyRez
User
User
Posty: 224
Rejestracja: poniedziałek 04 cze 2018, 09:19

Re: Cyfrowe przetwarzanie sygnałów

Postautor: GrumpyRez » wtorek 11 sty 2022, 15:08

W sumie ja w okolicach 2001-2004... przydało by się najpierw jakieś materiały wprowadzające przedstawić :mrgreen:

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » wtorek 11 sty 2022, 20:11

Zegar pisze:Ostatni raz liczyłem całki w ubiegłym wieku. :D

No ja to robią ciągle.
To nie jest jeszcze przesadna matma, to jeszcze w miarę proste rzeczy :D

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » wtorek 11 sty 2022, 23:14

Granice całkowania, coś o przedziałach

W dotychczasowych rozważaniach były stosowanie granice całkowania -π..π. Wcześniej uważałem, że zagadnienie należy rozpatrywać w granicach 0..2π, gdyż funkcja czasu dla argumentów ujemnych ma słabą interpretację fizyczną. To tak, jakbyśmy patrzyli od teraz w przyszłość. Obecnie raczej skłaniam się do poglądu, że dobrym rozwiązaniem jest przedział -π..π, gdyż bardziej uwidacznia się charakter parzystości i nieparzystości funkcji. To tak jakby spojrzeć w przeszłość i w przyszłość. W sumie, z punktu widzenia rachunkowego nie ma to żadnego znaczenia. By nie wdawać się z zawiłe rozkminy matematyczne, na tematykę można spojrzeć z następującej perspektywy. W rachunkach występuje całkowanie w określonych granicach. Niech będzie to przedział ( -T/2 .. T/2 ). Odpowiada to ilustracji 1.
sig07.png

Można „wyciąć” odpowiedni fragment i przemieścić go w inne miejsce, to uzyskuje się ponownie ten sam przebieg (ilustracja 2).
sig08.png

Tym razem rozpatrywany jest przedział (0,T).
Biorąc pod uwagę interpretację całki, że jest to suma iloczynów wartości funkcji okresowej i odpowiedniej funkcji trygonometrycznej oraz szerokości przedziału czasowego (dążącego do zera), jak na ilustracji 3, to
sig09.png

jeżeli cały przedział ( -T/2 .. T/2 ) zostanie poszatkowany na drobne odcinki czasowe i utworzone zostaną odpowiednie sumy, to powstaje pytanie: czym różnią się przebiegi między ilustracją 2 a ilustracją 1? Po chwili zastanowienia mamy w głowie odpowiedź: różnią się kolejnością składników sumy całki, nie różnią się wynikiem (sumowanie jest przemienne). Powstaje głębokie przekonanie, że nie jest istotne gdzie znajduje się przedział, z którego wynikają granice całkowania.
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » czwartek 13 sty 2022, 01:31

W stronę abstrakcji

Poprzez „zabieg rozłożenia” dowolnego sygnału okresowego na poszczególne harmoniczne (poszczególne częstotliwości) wchodzące w skład sygnału uzyskuje się przedstawienie tego sygnału jako sumy poszczególnych harmonicznych (częstotliwości składowych). Ponieważ „rachunki” dotyczące reakcji wzmacniacza na każdą częstotliwość (w sensie czystego sygnału sinusoidalnego o określonej częstotliwości) są znane, to można określić jak zachowa się dany układ (wzmacniacz) na wymuszenie dowolnym sygnałem. Wymaga to na początku rozłożenia sygnału na poszczególne „czyste” harmoniczne, przeprocesowanie każdej częstotliwości (odpowiada to reakcji układu [wzmacniacza] na każdą częstotliwość harmoniczną) i złożeniu poszczególnych harmonicznych do finalnego sygnału na wyjściu. Przetwarzanie sygnałów (w sensie DSP) podlega dokładnie tym samym regułom. Z tą różnicą, że w klasycznym (analogowym) wzmacniaczu procesowanie poszczególnych harmonicznych było realizowane w sposób analogowy przez elementy wchodzące w skład urządzenia (wzmacniacza), na które składają się przykładowo tranzystory, wzmacniacze itp. W ujęciu DSP, sygnał nie jest przetwarzany przez elementy wzmacniacza lecz jest przetwarzany przez procesor (czyli jakiś mały komputerek).
Wracając jeszcze do rozkładania sygnału analogowego na szereg trygonometryczny, na zagadnienie można spojrzeć z nieco innej perspektywy. Formułę:
part03_form001.png

można zapisać w postaci (gdzie b0=0, choć w sumie to nieistotne):
part03_form002.png

Nie zmienia to wyniku, gdyż cos(0 ω t)=1 oraz sin(0 ω t)=0, natomiast wpływa na pewną filozofię tematu: wprowadzona zostaje „zerowa” harmoniczna sygnału, gdyż całość można zapisać jako:
part03_form003.png

Trudno sobie wyobrazić harmoniczną o częstotliwości zero, to takie trochę mało fizyczne. Całe zjawisko zaczyna migrować do świata abstrakcji. Można tu puścić wodze fantazji: skoro jest częstotliwość zerowa, to może istnieją częstotliwości ujemne? Okazuje się, że jakkolwiek dziwnie to brzmi, to wszystko jest możliwe.
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » czwartek 13 sty 2022, 19:55

Kilka dobrze znanych przypadków

Weźmy dobrze znaną z świecie elektroniki funkcję piły. Niech jej szczegóły przedstawia ilustracja 1:
sig10.png

Dla prostoty rachunków, niech jej okres wynosi sekund, czyli pulsacja podstawowa ω=1 rad/s f=0.3183 hz. Składowa stała a0 to:
part04_from001.png

pozostałe współczynniki an szeregu trygonometrycznego to:
part04_from002.png

Skoro współczynniki an dla wszystkich n, z przedziału <0,+∞), to prowadzi do sugestii, że funkcja jest nieparzysta i jak popatrzeć na rysunek funkcji, to tak jest. Pozostałe współczynniki szeregu to:
part04_from003.png

Ponieważ cos(nπ)=1 dla n parzystych oraz cos(nπ)=-1 dla n nieparzystych, to ciąg współczynników bn przyjmuje wartości (kilka początkowych):
part04_from004.png

Funkcję piły można aproksymować następującym szeregiem trygonometrycznym:
part04_from005.png

Ograniczając się do początkowych kilku harmonicznych, powyższa funkcja ma następujący przebieg (ilustracja 2). Im więcej składników jest wziętych pod uwagę, tym szereg trygonometryczny dokładniej odtwarza funkcję oryginalną.
sig11.png

Innym „popularnym” sygnałem, który ma całkiem proste rachunki, jest sygnał fali prostokątnej. Niech sygnał przyjmuje wartości od -5V do 5V. Szczegóły pokazuje ilustracja 3.
sig12.png

Z parametrów funkcji wynika, że T=628 ms, f=1.59Hz, ω=10 rad/s.
Trochę rachunków. Składowa stała:
part04_from006.png

Współczynniki an są obliczone z formuły (A=5 jest amplitudą):
part04_from007.png

Współczynniki bn, to:
part04_from008.png

Wartość funkcji cos(nπ) przyjmuje wartości 1 dla n parzystych i -1 dla n nieparzystych.
Uwzględniając powyższą zależność współczynniki bn przyjmują wartość
part04_from009.png

i finalnie szereg ten można zapisać w postaci (przy założonych parametrach):
part04_from010.png

Funkcja fali prostokątnej aproksymowana kilkoma początkowymi harmonicznym ma przebieg pokazany na ilustracji 4.
sig13.png

Oczywiście, im tych harmonicznych jest więcej, tym szereg bardziej zbliża się do ideału.
Kilka innych „dobrze znanych” funkcji:
  • sygnał fali prostokątnej parzystej (ilustracja 5)
sig14.png

part04_from011.png

  • sygnał fali trójkątnej (ilustracja 6)
sig15.png

part04_from012.png

  • sygnał piły (ilustracja 7)
sig16.png

part04_from013.png
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » piątek 14 sty 2022, 00:48

Numeryczne wyznaczenie współczynników szeregu Fouriera

Pomimo, że analitycznie zostały określone współczynniki wszystkich wyrazów szeregu Fouriera, to wcale nie znaczy, że jest już z górki. Jeżeli jest znana wyjściowa funkcja czasu (w sensie zapisu analitycznego), to współczynniki szeregu dają się policzyć na papierze. Gorzej jest jak nie jest znana postać analityczna funkcji czasu a jedynie znane są jej wartości w ściśle określonych chwilach czasowych. To może być wynik próbkowania ze stałym interwałem czasowym przez przetwornik ADC, w wyniku którego uzyskuje się ciąg próbek. Mając tak „zadaną” funkcję (znane są jej wartości w obrębie jednego okresu w określonych chwilach czasowych) możliwe jest osiągnięcie celu, jakim jest postać szeregu trygonometrycznego. Można zaproponować przynajmniej dwie metody uzyskania wartości współczynników szeregu.
Pierwsze metoda to utworzenie układu równań i rozwiązanie go. To wyjątkowo „niesympatyczna” metoda, gdyż wymaga ogromnych nakładów i generuje olbrzymie koszty. Oprócz samych wad ma jedną istotną zaletę: wyjaśnia pewien istotny szczegół w obliczaniu współczynników, którego nie potrafię inaczej wyjaśnić a element jest istotny.
sig20.png

Rozpatrzmy sygnał analogowy o parametrach jak na ilustracji 1. Jeżeli sygnał został spróbkowany jak na rysunku, to siłą rzeczy musi zawierać zerową i pierwszą harmoniczną (twierdzenie Kotielnikowa-Shannona). Sprowadza się to do następującego zapisu:
part05_form01.png

W chwili t0,t1 i t2 zachodzą związki:
part05_form02.png

powstaje układ równań o trzech niewiadomych (a0, a1 i b1). Parametr Δp jest znany gdyż wynika z okresu próbkowania → znane są wartości funkcji trygonometrycznych. By układ o 3 niewiadomych miał jednoznaczne rozwiązanie muszą być 3 próbki.
W konkretnym przypadku (sygnału z ilustracji 1), jest to:
part05_form03.png

Jego rozwiązanie nie należy do skomplikowanych działań (użyłem arkusza kalkulacyjnego):
  • a_0 = 0.3
  • a_1 = 0.2
  • b_1 = 0.333333334
Jednak nie rozwiązanie jest tu najistotniejsze. Warto zauważyć, że by rozwiązać powyższy układ równań opisujący dwie harmoniczne (zerowa i pierwsza), konieczne są trzy próbki. Jeżeli tych harmonicznych będzie 2 (zerowa, pierwsza i druga), to koniecznych jest 5 próbek. Ekstrapolując zagadnienie, by określić n harmonicznych, to wymagane są 2n+1 próbek. Zawsze jest to nieparzysta liczba i to jest najważniejsza informacja wynikająca z tego sposobu rachunków. Oczywiście praktyczne wykorzystanie tej metody jest niewskazane, gdyż wymaga ogromnych nakładów (dla 50 harmonicznych konieczna jest kwadratowa macierz układu równań zawierająca 10201 elementów: 101X101) oraz ogromny koszt (obliczenie ponad 100 wyznaczników ogromnej macierzy).
Wróćmy do formuły pozwalającej na obliczenie współczynników rozwinięcia funkcji w szereg trygonometryczny:
part05_form04.png

Pozostaje jeszcze metoda obliczeń numerycznych. Wychodząc z interpretacji samej całki oznaczonej: pole figury zawartej pomiędzy osią czasu, krzywą funkcji w określonych granicach całkowania (jak na ilustracji 2), to
sig21.png

w ogólnym przypadku (całkowania jakiejś funkcji), należy utworzyć pola prostokątów i je zsumować.
part05_form05.png

W tym szczególnym przypadku (dotyczącym szeregu trygonometrycznego) sprawa się troszkę i komplikuje ale i upraszcza. Uproszczeniem jest przyrost argumentu Δt=1 jako, że próbki są ponumerowane niejako z krokiem 1. Komplikuje, bo należy brać pod uwagę iloczyn wartości próbki z funkcją sin oraz cos ale o częstotliwościach odpowiadającym określonym harmonicznych. W przypadku pierwszej harmonicznej dla współczynników stojących przy sin jest to iloczyn wartości próbki z jednookresową funkcją (ilustracja 3).
sig22.png

W przypadku drugiej harmonicznej jest to iloczyn wartości próbki z dwuokresową funkcją (ilustracja 4).
sig23.png

I tak dalej. W przypadku współczynników stojących przy funkcji cos, w iloczynie występuje odpowiedniookresowa funkcja cos. (przykładowo 3-okresowa funkcja cos dla wyznaczenia 3-iej harmonicznej, ilustracja 5).
sig25.png

Łatwo porachować, że tych całek jest 2n+1 (tyle ile wynosi liczba próbek).
Występująca tu metoda całkowania opiera się o prostokąty. Można wymyślić bardziej finezyjną metodę opierającą się o trapezy. Jej idea polega na tym, że zamiast pola prostokątów liczone są pola trapezów. Koncepcję pokazuje ilustracja 6.
sig24.png

Nawet nieuzbrojonym okiem widać, że jej dokładność powinna być większa, ale jest to okupione większym nakładem obliczeń i w sumie nie warta skórka za wyprawkę (sprawdziłem, nie warto).
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » piątek 14 sty 2022, 16:43

Coś, co tygryski lubią najbardziej - programy

Czas na słowo pisane (w języku C), eksperymenty i badania. Do powyższego algorytmu napisałem drobny programik w C (użyłem Dev-C++ ver. 5.11).
Program „wczytuje” (czytaj: ma podstawione gotowe dane) tablicę próbek i dokonuje odpowiednich rachunków. Jako przykład użyta została „dobrze znana funkcja”. Założyłem sobie, że pakiet danych stanowi 257 próbek indeksowanych w tablicy od indeksu=0 do indeksu=256. Próbki o numerach 0, 128 i 256 mają wartość 0, pozostałe wartość +5 lub -5 (jak na ilustracji 1). Praktycznie jest to funkcja prostokątna a zarazem ciągła (nie ma „skoków wartości”) i okresowa (koniec jednego okresu „styka” z początkiem kolejnego).
sig26.png

Program dokonuje analizy (rozkłada funkcję czasu na harmoniczne) i syntezy (składa z harmonicznych do funkcji czasu).
Nieparzystość liczby próbek jest osiągnięta następująco (deklaracja stałych i typów):
[code#define SampleNumDiv2 128
#define SampleNum (2*SampleNumDiv2)
#define Pi 3.141592654
#define DoublePi 6.283185307

typedef double real ;
//typedef float real ;

typedef real SampleArrayType [ SampleNum + 1 ] ;
typedef real CoofArrType [ SampleNumDiv2 + 1 ] ;[/code]
„Wczytanie” próbek realizuje poniższa funkcja (uzyskuje rezultat prezentowany przez powyższy rysunek → ilustracja 1):

Kod: Zaznacz cały

void InputSampleTable ( void )
{
  unsigned short Loop ;   
  /*------------------------------------------------------------------------*/   
  for ( Loop = 1 ; Loop < SampleNumDiv2 ; Loop ++ )
    InpSampleTable [ Loop ] = 5.0 ;
  InpSampleTable [ SampleNumDiv2 ] = 0.0 ;
  for ( Loop = SampleNumDiv2 + 1 ; Loop < SampleNum ; Loop ++ )
    InpSampleTable [ Loop ] = -5.0 ;
  InpSampleTable [ 0 ] = 0.0 ;
  InpSampleTable [ SampleNum ] = 0.0 ;
  T = 0.1 ; /* 100ms */
} /* InputSampleTable */

Funkcja analizy:

Kod: Zaznacz cały

void FourierSeries ( real ATable [ ] ,
                     real BTable [ ] ,
                     real SampleTable [ ] )
{
  unsigned short HarmInx ;
  unsigned short SampleInx ;
  real DAngle ;
  real Angle ;
  /*------------------------------------------------------------------------*/   
  for ( HarmInx = 0 ; HarmInx <= SampleNumDiv2 ; HarmInx ++ )
  {
    ATable [ HarmInx ] = 0.0 ;
    BTable [ HarmInx ] = 0.0 ;
  } /* for */ ;
  for ( SampleInx = 0 ; SampleInx <= SampleNum ; SampleInx ++ )
    ATable [ 0 ] = ATable [ 0 ] + SampleTable [ SampleInx ] ;
  for ( HarmInx = 1 ; HarmInx <= SampleNumDiv2 ; HarmInx ++ )
  {
    DAngle = DoublePi * ( real ) HarmInx / ( real ) SampleNum ;
    for ( SampleInx = 0 ; SampleInx <= SampleNum ; SampleInx ++ )
    {
      Angle = DAngle  * ( real ) SampleInx ;
      ATable [ HarmInx ] = ATable [ HarmInx ] + SampleTable [ SampleInx ] * cos ( Angle ) ;
      BTable [ HarmInx ] = BTable [ HarmInx ] + SampleTable [ SampleInx ] * sin ( Angle ) ;
    } /* for */ ;
  } /* for */
  ATable [ 0 ] = ATable [ 0 ] / ( real ) SampleNum ;
  BTable [ 0 ] = 0.0 ;
  for ( HarmInx = 1 ; HarmInx <= SampleNumDiv2 ; HarmInx ++ )
  {
    ATable [ HarmInx ] = ATable [ HarmInx ] / ( real ) SampleNumDiv2 ;
    BTable [ HarmInx ] = BTable [ HarmInx ] / ( real ) SampleNumDiv2 ;
  } /* for */ ; 
} /* FourierSeries */

W powyższym przykładzie oddzielnie jest porachowana składowa stała:

Kod: Zaznacz cały

void FourierSeries ( real ATable [ ] ,
                     real BTable [ ] ,
                     real SampleTable [ ] )
{
( . . )
  for ( SampleInx = 0 ; SampleInx <= SampleNum ; SampleInx ++ )
    ATable [ 0 ] = ATable [ 0 ] + SampleTable [ SampleInx ] ;
  for ( HarmInx = 1 ; HarmInx <= SampleNumDiv2 ; HarmInx ++ )
  {
( . . )
  ATable [ 0 ] = ATable [ 0 ] / ( real ) SampleNum ;
  BTable [ 0 ] = 0.0 ;
( . . )
} /* FourierSeries */

chociaż z „innego punktu widzenia” (potraktować wartość średnią jako zerową harmoniczną) można ją „włączyć” w główną pętlę algorytmu (usunąć pętlę sumującą wszystkie próbki do elementu ATable[0] oraz zacząć pętlę algorytmu obliczającego współczynniki dla harmonicznej od indeksu 0):

Kod: Zaznacz cały

void FourierSeries ( real ATable [ ] ,
                     real BTable [ ] ,
                     real SampleTable [ ] )
{
( . . )
  for ( HarmInx = 0 ; HarmInx <= SampleNumDiv2 ; HarmInx ++ )
  {
( . . )
  ATable [ 0 ] = ATable [ 0 ] / ( real ) SampleNum ;
  BTable [ 0 ] = 0.0 ;
( . . )
} /* FourierSeries */

Wynik wychodzi taki sam, a jedyna różnica dotyczy „kosztów” wykonania programu [nie zachodzi potrzeba wielokrotnego liczenia sin(0) oraz cos (0)].
Funkcja syntezy:

Kod: Zaznacz cały

void TestFourierSeries ( void )
{
  unsigned short SampleInx ;
  unsigned short HarmInx ;
  real SigVal ;
  real DAngle ;
  real Angle ;
  real DeltaVal ;
  real Percent ;
  /*------------------------------------------------------------------------*/   
  printf ( "\n\nTesty\n" ) ;
  for ( SampleInx = 0 ; SampleInx <= SampleNum ; SampleInx ++ )
  {
    DAngle = DoublePi * ( real ) SampleInx / ( real ) SampleNum ;
    SigVal = OutATable [ 0 ] ;
    for ( HarmInx = 1 ; HarmInx <= SampleNumDiv2 ; HarmInx ++ )
    {
      Angle  = DAngle * ( real ) HarmInx ;
      SigVal = SigVal + OutATable [ HarmInx ] * cos ( Angle ) + OutBTable [ HarmInx ] * sin ( Angle ) ;
    } /* for */ ;
    DeltaVal = ( real ) InpSampleTable [ SampleInx ] - SigVal ;
    printf ( "Probka [ %d ] = %f, obl. wart.=%f [rozn.=%f" , SampleInx , InpSampleTable [ SampleInx ] , SigVal , DeltaVal ) ;
    if ( InpSampleTable [ SampleInx ] )
    {
      Percent = 100.0 * DeltaVal / ( real ) InpSampleTable [ SampleInx ] ;
      printf ( " :%f%%" , Percent ) ;       
   } /* if */ ;
    printf ( "]\n" ) ;
} /* for */ ;
} /* TestFourierSeries */

Funkcja main to:

Kod: Zaznacz cały

int main ( int argc , char ** argv )
{
  unsigned short Loop ;
  /*------------------------------------------------------------------------*/   
  InputSampleTable ( ) ;
  for ( Loop = 0 ; Loop <= SampleNum ; Loop ++ )
    printf ( "Probka [ %d ] = %f\n" , Loop , InpSampleTable [ Loop ] ) ;
  printf ( "Liczba probek: %d\n" , SampleNum + 1 ) ;
  printf ( "Okres sygnalu: %f ms\n" , 1000.0 * T ) ;
  printf ( "**********************************\r\n" ) ;
  printf ( "Okres probkowania: %f ms\n" , 1000.0 * T / ( real ) SampleNum ) ;
  printf ( "Pulsacja podstawowa: %f rad/s\n" ,  DoublePi / T ) ;
  printf ( "Czestotliwosc podstawowa: %f Hz\n\n" ,  1.0 / T ) ;
  FourierSeries ( OutATable , OutBTable , InpSampleTable ) ;
  printf ( "Skladowa stala: A.0= %f\n" ,  OutATable [ 0 ] ) ;
  for ( Loop = 1 ; Loop <= SampleNumDiv2 ; Loop ++ )
  {
    printf ( "Harm. %d:    A.%d=%f    B.%d=%f\n" , Loop , Loop , OutATable [ Loop ] , Loop , OutBTable [ Loop ] ) ;
  } /* for */ ;
  TestFourierSeries ( ) ;
} /* main */

Wyniki badań
Jeżeli przyjąć, że pobranie wszystkich próbek zajęło T=0.1s, to

Kod: Zaznacz cały

Okres sygnalu: 100.000000 ms
Okres probkowania: 0.390625 ms
Pulsacja podstawowa: 62.831853 rad/s
Czestotliwosc podstawowa: 10.000000 Hz/s

Z założonej liczby próbek daje się uzyskać 128 harmonicznych (od harmonicznej zerowej do harmonicznej o numerze 127). Są one następujące:

Kod: Zaznacz cały

Skladowa stala: A.0= 0.000000
Harm. 1:      A.1= 0.000000     B.1= 6.365878
Harm. 2:      A.2=-0.000000     B.2=-0.000000
Harm. 3:      A.3= 0.000000     B.3= 2.121107
Harm. 4:      A.4=-0.000000     B.4=-0.000000
Harm. 5:      A.5= 0.000000     B.5= 1.271641
Harm. 6:      A.6= 0.000000     B.6=-0.000000
Harm. 7:      A.7= 0.000000     B.7= 0.907219
Harm. 8:      A.8=-0.000000     B.8=-0.000000
Harm. 9:      A.9= 0.000000     B.9= 0.704477
Harm. 10:    A.10=-0.000000    B.10=-0.000000
Harm. 11:    A.11= 0.000000    B.11= 0.575226
Harm. 12:    A.12=-0.000000    B.12=-0.000000
Harm. 13:    A.13= 0.000000    B.13= 0.485546
Harm. 14:    A.14=-0.000000    B.14=-0.000000
Harm. 15:    A.15= 0.000000    B.15= 0.419609
Harm. 16:    A.16= 0.000000    B.16=-0.000000
Harm. 17:    A.17= 0.000000    B.17= 0.369034
Harm. 18:    A.18= 0.000000    B.18=-0.000000
Harm. 19:    A.19= 0.000000    B.19= 0.328969
Harm. 20:    A.20=-0.000000    B.20=-0.000000
Harm. 21:    A.21= 0.000000    B.21= 0.296411
Harm. 22:    A.22= 0.000000    B.22=-0.000000
Harm. 23:    A.23= 0.000000    B.23= 0.269402
Harm. 24:    A.24=-0.000000    B.24=-0.000000
Harm. 25:    A.25= 0.000000    B.25= 0.246608
Harm. 26:    A.26= 0.000000    B.26=-0.000000
Harm. 27:    A.27= 0.000000    B.27= 0.227093
Harm. 28:    A.28=-0.000000    B.28=-0.000000
Harm. 29:    A.29= 0.000000    B.29= 0.210177
Harm. 30:    A.30=-0.000000    B.30=-0.000000
Harm. 31:    A.31= 0.000000    B.31= 0.195357
Harm. 32:    A.32= 0.000000    B.32=-0.000000
Harm. 33:    A.33= 0.000000    B.33= 0.182252
Harm. 34:    A.34= 0.000000    B.34=-0.000000
Harm. 35:    A.35= 0.000000    B.35= 0.170566
Harm. 36:    A.36= 0.000000    B.36=-0.000000
Harm. 37:    A.37= 0.000000    B.37= 0.160069
Harm. 38:    A.38= 0.000000    B.38=-0.000000
Harm. 39:    A.39= 0.000000    B.39= 0.150578
Harm. 40:    A.40=-0.000000    B.40=-0.000000
Harm. 41:    A.41= 0.000000    B.41= 0.141944
Harm. 42:    A.42=-0.000000    B.42=-0.000000
Harm. 43:    A.43= 0.000000    B.43= 0.134047
Harm. 44:    A.44=-0.000000    B.44=-0.000000
Harm. 45:    A.45= 0.000000    B.45= 0.126789
Harm. 46:    A.46=-0.000000    B.46=-0.000000
Harm. 47:    A.47= 0.000000    B.47= 0.120087
Harm. 48:    A.48= 0.000000    B.48=-0.000000
Harm. 49:    A.49= 0.000000    B.49= 0.113872
Harm. 50:    A.50= 0.000000    B.50=-0.000000
Harm. 51:    A.51= 0.000000    B.51= 0.108087
Harm. 52:    A.52=-0.000000    B.52=-0.000000
Harm. 53:    A.53= 0.000000    B.53= 0.102681
Harm. 54:    A.54=-0.000000    B.54=-0.000000
Harm. 55:    A.55= 0.000000    B.55= 0.097614
Harm. 56:    A.56= 0.000000    B.56=-0.000000
Harm. 57:    A.57= 0.000000    B.57= 0.092848
Harm. 58:    A.58= 0.000000    B.58=-0.000000
Harm. 59:    A.59= 0.000000    B.59= 0.088353
Harm. 60:    A.60= 0.000000    B.60=-0.000000
Harm. 61:    A.61= 0.000000    B.61= 0.084100
Harm. 62:    A.62=-0.000000    B.62=-0.000000
Harm. 63:    A.63= 0.000000    B.63= 0.080066
Harm. 64:    A.64= 0.000000    B.64=-0.000000
Harm. 65:    A.65= 0.000000    B.65= 0.076231
Harm. 66:    A.66= 0.000000    B.66=-0.000000
Harm. 67:    A.67= 0.000000    B.67= 0.072574
Harm. 68:    A.68= 0.000000    B.68=-0.000000
Harm. 69:    A.69= 0.000000    B.69= 0.069081
Harm. 70:    A.70=-0.000000    B.70=-0.000000
Harm. 71:    A.71= 0.000000    B.71= 0.065736
Harm. 72:    A.72=-0.000000    B.72=-0.000000
Harm. 73:    A.73= 0.000000    B.73= 0.062527
Harm. 74:    A.74= 0.000000    B.74=-0.000000
Harm. 75:    A.75= 0.000000    B.75= 0.059441
Harm. 76:    A.76=-0.000000    B.76=-0.000000
Harm. 77:    A.77= 0.000000    B.77= 0.056469
Harm. 78:    A.78= 0.000000    B.78=-0.000000
Harm. 79:    A.79= 0.000000    B.79= 0.053600
Harm. 80:    A.80= 0.000000    B.80=-0.000000
Harm. 81:    A.81= 0.000000    B.81= 0.050826
Harm. 82:    A.82=-0.000000    B.82=-0.000000
Harm. 83:    A.83= 0.000000    B.83= 0.048139
Harm. 84:    A.84=-0.000000    B.84=-0.000000
Harm. 85:    A.85= 0.000000    B.85= 0.045533
Harm. 86:    A.86= 0.000000    B.86=-0.000000
Harm. 87:    A.87= 0.000000    B.87= 0.043000
Harm. 88:    A.88=-0.000000    B.88=-0.000000
Harm. 89:    A.89= 0.000000    B.89= 0.040534
Harm. 90:    A.90=-0.000000    B.90=-0.000000
Harm. 91:    A.91= 0.000000    B.91= 0.038130
Harm. 92:    A.92= 0.000000    B.92=-0.000000
Harm. 93:    A.93= 0.000000    B.93= 0.035784
Harm. 94:    A.94=-0.000000    B.94=-0.000000
Harm. 95:    A.95= 0.000000    B.95= 0.033489
Harm. 96:    A.96= 0.000000    B.96=-0.000000
Harm. 97:    A.97= 0.000000    B.97= 0.031243
Harm. 98:    A.98=-0.000000    B.98=-0.000000
Harm. 99:    A.99= 0.000000    B.99= 0.029040
Harm. 100:   A.100=-0.000000   B.100=-0.000000
Harm. 101:   A.101= 0.000000   B.101= 0.026877
Harm. 102:   A.102= 0.000000   B.102=-0.000000
Harm. 103:   A.103= 0.000000   B.103= 0.024750
Harm. 104:   A.104= 0.000000   B.104=-0.000000
Harm. 105:   A.105= 0.000000   B.105= 0.022656
Harm. 106:   A.106=-0.000000   B.106=-0.000000
Harm. 107:   A.107= 0.000000   B.107= 0.020591
Harm. 108:   A.108=-0.000000   B.108=-0.000000
Harm. 109:   A.109= 0.000000   B.109= 0.018553
Harm. 110:   A.110=-0.000000   B.110=-0.000000
Harm. 111:   A.111= 0.000000   B.111= 0.016539
Harm. 112:   A.112= 0.000000   B.112=-0.000000
Harm. 113:   A.113= 0.000000   B.113= 0.014546
Harm. 114:   A.114= 0.000000   B.114=-0.000000
Harm. 115:   A.115= 0.000000   B.115= 0.012570
Harm. 116:   A.116=-0.000000   B.116=-0.000000
Harm. 117:   A.117= 0.000000   B.117= 0.010611
Harm. 118:   A.118= 0.000000   B.118=-0.000000
Harm. 119:   A.119= 0.000000   B.119= 0.008664
Harm. 120:   A.120= 0.000000   B.120=-0.000000
Harm. 121:   A.121= 0.000000   B.121= 0.006728
Harm. 122:   A.122= 0.000000   B.122=-0.000000
Harm. 123:   A.123= 0.000000   B.123= 0.004800
Harm. 124:   A.124= 0.000000   B.124=-0.000000
Harm. 125:   A.125= 0.000000   B.125= 0.002878
Harm. 126:   A.126= 0.000000   B.126=-0.000000
Harm. 127:   A.127= 0.000000   B.127= 0.000959
Harm. 128:   A.128=-0.000000   B.128=-0.000000

Nawet nie trzeba wkładać okularów, by zauważyć, że wszystkie współczynniki an mają wartość zerową oraz współczynniki bn zerują się w przypadku parzystych numerów harmonicznych. Ich wartości są zgodne z przewidywaniami teoretycznymi, kilka początkowych (dla sygnału o założonych parametrach: amplitudzie =5 → bn=20/(n π) ):
  • jest B.1= 6.365878 powinno być: 20/π=6.366197724 (w zaokrągleniu:6.366198)
  • jest B.3= 2.121107 powinno być: 20/3π=2.122065908 (w zaokrągleniu:2.122066)
  • jest B.5= 1.271641 powinno być: 20/5π=1.273239545 (w zaokrągleniu:1.273240)
  • jest B.7= 0.907219 powinno być: 20/7π=0.909456817 (w zaokrągleniu:0.909457)
  • jest B.9= 0.704477 powinno być: 20/9π=0.707355302 (w zaokrągleniu:0.707355)
  • jest B.11= 0.575226 powinno być: 20/11π=0.578745247 (w zaokrągleniu:0.578745)
Sądzę, że niewielkie różnice są spowodowane tym, że funkcja nie ma idealnie stromych zboczy, ale jest „ciągła” i okresowa.
Synteza (złożenie z harmonicznych sygnału wyjściowego), program policzył wartość próbki, obliczył różnicę od wartości założonej oraz określił błąd w procentach:

Kod: Zaznacz cały

Probka [ 0 ] = 0.000000, obl. wart.=0.000000 [rozn.=-0.000000]
Probka [ 1 ] = 5.000000, obl. wart.=5.000000 [rozn.=0.000000 :0.000000%]
Probka [ 2 ] = 5.000000, obl. wart.=5.000000 [rozn.=-0.000000 :-0.000000%]
Probka [ 3 ] = 5.000000, obl. wart.=5.000000 [rozn.=0.000000 :0.000000%]
Probka [ 4 ] = 5.000000, obl. wart.=5.000000 [rozn.=-0.000000 :-0.000000%]
Probka [ 5 ] = 5.000000, obl. wart.=5.000000 [rozn.=0.000000 :0.000000%]
Probka [ 6 ] = 5.000000, obl. wart.=5.000000 [rozn.=-0.000000 :-0.000000%]
Probka [ 7 ] = 5.000000, obl. wart.=5.000000 [rozn.=0.000000 :0.000000%]
Probka [ 8 ] = 5.000000, obl. wart.=5.000000 [rozn.=-0.000000 :-0.000000%]
Probka [ 9 ] = 5.000000, obl. wart.=5.000000 [rozn.=0.000000 :0.000000%]
Probka [ 10 ] = 5.000000, obl. wart.=5.000000 [rozn.=-0.000000 :-0.000000%]
Probka [ 11 ] = 5.000000, obl. wart.=5.000000 [rozn.=0.000000 :0.000000%]
Probka [ 12 ] = 5.000000, obl. wart.=5.000000 [rozn.=-0.000000 :-0.000000%]
( . . . )
Probka [ 248 ] = -5.000000, obl. wart.=-5.000000 [rozn.=0.000000 :-0.000000%]
Probka [ 249 ] = -5.000000, obl. wart.=-5.000000 [rozn.=-0.000000 :0.000000%]
Probka [ 250 ] = -5.000000, obl. wart.=-5.000000 [rozn.=0.000000 :-0.000000%]
Probka [ 251 ] = -5.000000, obl. wart.=-5.000000 [rozn.=-0.000000 :0.000000%]
Probka [ 252 ] = -5.000000, obl. wart.=-5.000000 [rozn.=0.000000 :-0.000000%]
Probka [ 253 ] = -5.000000, obl. wart.=-5.000000 [rozn.=-0.000000 :0.000000%]
Probka [ 254 ] = -5.000000, obl. wart.=-5.000000 [rozn.=0.000000 :-0.000000%]
Probka [ 255 ] = -5.000000, obl. wart.=-5.000000 [rozn.=-0.000000 :0.000000%]
Probka [ 256 ] = 0.000000, obl. wart.=-0.000000 [rozn.=0.000000]

Generalnie nie należy popadać w przesadną euforię, gdyż obliczenia były przeprowadzone dla liczb o podwójnej precyzji. W programie można „przestawić” rachunki na liczby o pojedynczej precyzji (wystarczy zakometarzować i odkomentarzować odpowiednie deklaracje). Nasuwa się prosty wniosek: zjawisko jest bardzo subtelne.

Kod: Zaznacz cały

//typedef double real ;
typedef float real ;

Informacje wygenerowane przez program dla wariantu pojedynczej precyzji.

Kod: Zaznacz cały

Probka [ 0 ] = 0.000000, obl. wart.=-0.000015 [rozn.=0.000015]
Probka [ 1 ] = 5.000000, obl. wart.=4.999980 [rozn.=0.000020 :0.000391%]
Probka [ 2 ] = 5.000000, obl. wart.=5.000004 [rozn.=-0.000004 :-0.000076%]
Probka [ 3 ] = 5.000000, obl. wart.=5.000036 [rozn.=-0.000036 :-0.000715%]
Probka [ 4 ] = 5.000000, obl. wart.=4.999989 [rozn.=0.000011 :0.000219%]
Probka [ 5 ] = 5.000000, obl. wart.=5.000013 [rozn.=-0.000013 :-0.000257%]
Probka [ 6 ] = 5.000000, obl. wart.=4.999978 [rozn.=0.000022 :0.000439%]
Probka [ 7 ] = 5.000000, obl. wart.=5.000003 [rozn.=-0.000003 :-0.000067%]
Probka [ 8 ] = 5.000000, obl. wart.=4.999963 [rozn.=0.000037 :0.000744%]
Probka [ 9 ] = 5.000000, obl. wart.=5.000055 [rozn.=-0.000055 :-0.001097%]
Probka [ 10 ] = 5.000000, obl. wart.=4.999961 [rozn.=0.000039 :0.000782%]
Probka [ 11 ] = 5.000000, obl. wart.=5.000040 [rozn.=-0.000040 :-0.000801%]
Probka [ 12 ] = 5.000000, obl. wart.=4.999970 [rozn.=0.000030 :0.000591%]
Probka [ 13 ] = 5.000000, obl. wart.=5.000015 [rozn.=-0.000015 :-0.000305%]
Probka [ 14 ] = 5.000000, obl. wart.=4.999963 [rozn.=0.000037 :0.000734%]
Probka [ 15 ] = 5.000000, obl. wart.=5.000038 [rozn.=-0.000038 :-0.000763%]
Probka [ 16 ] = 5.000000, obl. wart.=4.999952 [rozn.=0.000048 :0.000954%]
Probka [ 17 ] = 5.000000, obl. wart.=5.000053 [rozn.=-0.000053 :-0.001059%]
Probka [ 18 ] = 5.000000, obl. wart.=4.999915 [rozn.=0.000085 :0.001707%]
Probka [ 19 ] = 5.000000, obl. wart.=5.000084 [rozn.=-0.000084 :-0.001678%]
Probka [ 20 ] = 5.000000, obl. wart.=4.999926 [rozn.=0.000074 :0.001478%]
Probka [ 21 ] = 5.000000, obl. wart.=5.000067 [rozn.=-0.000067 :-0.001335%]
Probka [ 22 ] = 5.000000, obl. wart.=4.999931 [rozn.=0.000069 :0.001373%]
Probka [ 23 ] = 5.000000, obl. wart.=5.000038 [rozn.=-0.000038 :-0.000753%]
Probka [ 24 ] = 5.000000, obl. wart.=4.999985 [rozn.=0.000015 :0.000296%]
Probka [ 25 ] = 5.000000, obl. wart.=5.000007 [rozn.=-0.000007 :-0.000143%]
Probka [ 26 ] = 5.000000, obl. wart.=4.999986 [rozn.=0.000014 :0.000286%]
Probka [ 27 ] = 5.000000, obl. wart.=5.000038 [rozn.=-0.000038 :-0.000763%]
Probka [ 28 ] = 5.000000, obl. wart.=4.999957 [rozn.=0.000043 :0.000868%]
Probka [ 29 ] = 5.000000, obl. wart.=5.000041 [rozn.=-0.000041 :-0.000811%]
Probka [ 30 ] = 5.000000, obl. wart.=4.999983 [rozn.=0.000017 :0.000343%]
Probka [ 31 ] = 5.000000, obl. wart.=4.999973 [rozn.=0.000027 :0.000544%]
Probka [ 32 ] = 5.000000, obl. wart.=5.000018 [rozn.=-0.000018 :-0.000362%]
Probka [ 33 ] = 5.000000, obl. wart.=4.999975 [rozn.=0.000025 :0.000505%]
Probka [ 34 ] = 5.000000, obl. wart.=5.000023 [rozn.=-0.000023 :-0.000467%]
Probka [ 35 ] = 5.000000, obl. wart.=4.999992 [rozn.=0.000008 :0.000162%]
Probka [ 36 ] = 5.000000, obl. wart.=5.000018 [rozn.=-0.000018 :-0.000362%]
Probka [ 37 ] = 5.000000, obl. wart.=4.999981 [rozn.=0.000019 :0.000381%]
Probka [ 38 ] = 5.000000, obl. wart.=5.000016 [rozn.=-0.000016 :-0.000324%]
Probka [ 39 ] = 5.000000, obl. wart.=4.999985 [rozn.=0.000015 :0.000305%]
Probka [ 40 ] = 5.000000, obl. wart.=5.000018 [rozn.=-0.000018 :-0.000353%]
Probka [ 41 ] = 5.000000, obl. wart.=4.999984 [rozn.=0.000016 :0.000315%]
Probka [ 42 ] = 5.000000, obl. wart.=5.000009 [rozn.=-0.000009 :-0.000172%]
Probka [ 43 ] = 5.000000, obl. wart.=4.999979 [rozn.=0.000021 :0.000429%]
Probka [ 44 ] = 5.000000, obl. wart.=4.999986 [rozn.=0.000014 :0.000277%]
Probka [ 45 ] = 5.000000, obl. wart.=5.000032 [rozn.=-0.000032 :-0.000648%]
Probka [ 46 ] = 5.000000, obl. wart.=4.999997 [rozn.=0.000003 :0.000067%]
Probka [ 47 ] = 5.000000, obl. wart.=5.000021 [rozn.=-0.000021 :-0.000429%]
Probka [ 48 ] = 5.000000, obl. wart.=5.000003 [rozn.=-0.000003 :-0.000067%]
Probka [ 49 ] = 5.000000, obl. wart.=4.999976 [rozn.=0.000024 :0.000486%]
Probka [ 50 ] = 5.000000, obl. wart.=5.000015 [rozn.=-0.000015 :-0.000296%]
Probka [ 51 ] = 5.000000, obl. wart.=4.999980 [rozn.=0.000020 :0.000391%]
Probka [ 52 ] = 5.000000, obl. wart.=5.000044 [rozn.=-0.000044 :-0.000887%]
Probka [ 53 ] = 5.000000, obl. wart.=4.999967 [rozn.=0.000033 :0.000658%]
Probka [ 54 ] = 5.000000, obl. wart.=5.000015 [rozn.=-0.000015 :-0.000305%]
Probka [ 55 ] = 5.000000, obl. wart.=4.999960 [rozn.=0.000040 :0.000792%]
Probka [ 56 ] = 5.000000, obl. wart.=5.000023 [rozn.=-0.000023 :-0.000458%]
Probka [ 57 ] = 5.000000, obl. wart.=4.999984 [rozn.=0.000016 :0.000324%]
Probka [ 58 ] = 5.000000, obl. wart.=5.000022 [rozn.=-0.000022 :-0.000439%]
Probka [ 59 ] = 5.000000, obl. wart.=4.999972 [rozn.=0.000028 :0.000553%]
Probka [ 60 ] = 5.000000, obl. wart.=5.000062 [rozn.=-0.000062 :-0.001230%]
Probka [ 61 ] = 5.000000, obl. wart.=4.999938 [rozn.=0.000062 :0.001230%]
Probka [ 62 ] = 5.000000, obl. wart.=5.000040 [rozn.=-0.000040 :-0.000801%]
Probka [ 63 ] = 5.000000, obl. wart.=4.999969 [rozn.=0.000031 :0.000620%]
Probka [ 64 ] = 5.000000, obl. wart.=5.000017 [rozn.=-0.000017 :-0.000334%]
Probka [ 65 ] = 5.000000, obl. wart.=5.000000 [rozn.=0.000000 :0.000000%]
Probka [ 66 ] = 5.000000, obl. wart.=4.999994 [rozn.=0.000006 :0.000114%]
Probka [ 67 ] = 5.000000, obl. wart.=4.999985 [rozn.=0.000015 :0.000296%]
Probka [ 68 ] = 5.000000, obl. wart.=5.000026 [rozn.=-0.000026 :-0.000525%]
Probka [ 69 ] = 5.000000, obl. wart.=4.999955 [rozn.=0.000045 :0.000896%]
Probka [ 70 ] = 5.000000, obl. wart.=5.000064 [rozn.=-0.000064 :-0.001278%]
Probka [ 71 ] = 5.000000, obl. wart.=4.999962 [rozn.=0.000038 :0.000763%]
Probka [ 72 ] = 5.000000, obl. wart.=5.000012 [rozn.=-0.000012 :-0.000238%]
Probka [ 73 ] = 5.000000, obl. wart.=5.000021 [rozn.=-0.000021 :-0.000410%]
Probka [ 74 ] = 5.000000, obl. wart.=4.999968 [rozn.=0.000032 :0.000648%]
Probka [ 75 ] = 5.000000, obl. wart.=5.000039 [rozn.=-0.000039 :-0.000772%]
Probka [ 76 ] = 5.000000, obl. wart.=4.999985 [rozn.=0.000015 :0.000296%]
Probka [ 77 ] = 5.000000, obl. wart.=4.999987 [rozn.=0.000013 :0.000267%]
Probka [ 78 ] = 5.000000, obl. wart.=5.000042 [rozn.=-0.000042 :-0.000839%]
Probka [ 79 ] = 5.000000, obl. wart.=4.999978 [rozn.=0.000022 :0.000448%]
Probka [ 80 ] = 5.000000, obl. wart.=4.999986 [rozn.=0.000014 :0.000277%]
Probka [ 81 ] = 5.000000, obl. wart.=4.999986 [rozn.=0.000014 :0.000277%]
Probka [ 82 ] = 5.000000, obl. wart.=5.000002 [rozn.=-0.000002 :-0.000038%]
Probka [ 83 ] = 5.000000, obl. wart.=5.000032 [rozn.=-0.000032 :-0.000639%]
Probka [ 84 ] = 5.000000, obl. wart.=4.999983 [rozn.=0.000017 :0.000343%]
Probka [ 85 ] = 5.000000, obl. wart.=4.999993 [rozn.=0.000007 :0.000134%]
Probka [ 86 ] = 5.000000, obl. wart.=5.000000 [rozn.=-0.000000 :-0.000010%]
Probka [ 87 ] = 5.000000, obl. wart.=5.000035 [rozn.=-0.000035 :-0.000696%]
Probka [ 88 ] = 5.000000, obl. wart.=4.999969 [rozn.=0.000031 :0.000629%]
Probka [ 89 ] = 5.000000, obl. wart.=5.000029 [rozn.=-0.000029 :-0.000572%]
Probka [ 90 ] = 5.000000, obl. wart.=4.999993 [rozn.=0.000007 :0.000143%]
Probka [ 91 ] = 5.000000, obl. wart.=5.000000 [rozn.=-0.000000 :-0.000010%]
Probka [ 92 ] = 5.000000, obl. wart.=4.999980 [rozn.=0.000020 :0.000401%]
Probka [ 93 ] = 5.000000, obl. wart.=5.000029 [rozn.=-0.000029 :-0.000582%]
Probka [ 94 ] = 5.000000, obl. wart.=4.999952 [rozn.=0.000048 :0.000963%]
Probka [ 95 ] = 5.000000, obl. wart.=5.000021 [rozn.=-0.000021 :-0.000420%]
Probka [ 96 ] = 5.000000, obl. wart.=5.000000 [rozn.=0.000000 :0.000000%]
Probka [ 97 ] = 5.000000, obl. wart.=5.000009 [rozn.=-0.000009 :-0.000181%]
Probka [ 98 ] = 5.000000, obl. wart.=4.999989 [rozn.=0.000011 :0.000219%]
Probka [ 99 ] = 5.000000, obl. wart.=5.000036 [rozn.=-0.000036 :-0.000725%]
Probka [ 100 ] = 5.000000, obl. wart.=4.999961 [rozn.=0.000039 :0.000772%]
(...)
Probka [ 241 ] = -5.000000, obl. wart.=-5.000075 [rozn.=0.000075 :-0.001507%]
Probka [ 242 ] = -5.000000, obl. wart.=-4.999925 [rozn.=-0.000075 :0.001507%]
Probka [ 243 ] = -5.000000, obl. wart.=-5.000043 [rozn.=0.000043 :-0.000858%]
Probka [ 244 ] = -5.000000, obl. wart.=-4.999953 [rozn.=-0.000047 :0.000935%]
Probka [ 245 ] = -5.000000, obl. wart.=-5.000027 [rozn.=0.000027 :-0.000534%]
Probka [ 246 ] = -5.000000, obl. wart.=-5.000006 [rozn.=0.000006 :-0.000124%]
Probka [ 247 ] = -5.000000, obl. wart.=-4.999997 [rozn.=-0.000003 :0.000067%]
Probka [ 248 ] = -5.000000, obl. wart.=-4.999988 [rozn.=-0.000012 :0.000238%]
Probka [ 249 ] = -5.000000, obl. wart.=-5.000024 [rozn.=0.000024 :-0.000486%]
Probka [ 250 ] = -5.000000, obl. wart.=-4.999988 [rozn.=-0.000012 :0.000238%]
Probka [ 251 ] = -5.000000, obl. wart.=-5.000017 [rozn.=0.000017 :-0.000343%]
Probka [ 252 ] = -5.000000, obl. wart.=-4.999979 [rozn.=-0.000021 :0.000429%]
Probka [ 253 ] = -5.000000, obl. wart.=-4.999996 [rozn.=-0.000004 :0.000086%]
Probka [ 254 ] = -5.000000, obl. wart.=-5.000052 [rozn.=0.000052 :-0.001040%]
Probka [ 255 ] = -5.000000, obl. wart.=-4.999939 [rozn.=-0.000061 :0.001211%]
Probka [ 256 ] = 0.000000, obl. wart.=0.000034 [rozn.=-0.000034]

Poza tym, można dostrzec w raporcie dotyczącym wartości współczynników kolejnych harmonicznych, że pojawiają się wartości niezerowe (co prawda, pierwsza niezerowa cyfra jest praktycznie na 6 miejscu po przecinku). Biorąc pod uwagę, że w systemach prockowych próbkowanie jest całkowitoliczbowe, to jestem przekonany, że nie ma to finalnego znaczenia. Warto mieć świadomość, że pewne zjawiska zachodzą.
Dla porównania warto sprawdzić zachowanie się programu w sytuacji przetwarzania parzystej liczby próbek. Różnic w algorytmie nie ma, jedyna różnica wynika z deklaracji typu

Kod: Zaznacz cały

typedef real SampleArrayType [ SampleNum ] ;   // wariant parzysty
typedef real SampleArrayType [ SampleNum+1 ] ; // wariant nieparzysty

oraz organizacji pętli obliczeniowych

Kod: Zaznacz cały

for ( SampleInx = 0 ; SampleInx < SampleNum ; SampleInx ++ ) // wariant parzysty
for ( SampleInx = 0 ; SampleInx <= SampleNum ; SampleInx ++ ) // wariant nieparzysty

Raport wygenerowany przez program:

Kod: Zaznacz cały

Skladowa stala: A.0= 0.000000
Harm. 1:    A.1=0.078125    B.1=6.365878
Harm. 2:    A.2=-0.000000    B.2=-0.000000
Harm. 3:    A.3=0.078125    B.3=2.121107
Harm. 4:    A.4=-0.000000    B.4=-0.000000
Harm. 5:    A.5=0.078125    B.5=1.271641
Harm. 6:    A.6=0.000000    B.6=-0.000000
Harm. 7:    A.7=0.078125    B.7=0.907219
Harm. 8:    A.8=-0.000000    B.8=-0.000000
Harm. 9:    A.9=0.078125    B.9=0.704477
Harm. 10:    A.10=-0.000000    B.10=-0.000000
Harm. 11:    A.11=0.078125    B.11=0.575226
Harm. 12:    A.12=-0.000000    B.12=0.000000
Harm. 13:    A.13=0.078125    B.13=0.485546
Harm. 14:    A.14=-0.000000    B.14=-0.000000
Harm. 15:    A.15=0.078125    B.15=0.419609
Harm. 16:    A.16=0.000000    B.16=0.000000
Harm. 17:    A.17=0.078125    B.17=0.369034
Harm. 18:    A.18=0.000000    B.18=-0.000000
Harm. 19:    A.19=0.078125    B.19=0.328969
Harm. 20:    A.20=-0.000000    B.20=-0.000000
Harm. 21:    A.21=0.078125    B.21=0.296411
Harm. 22:    A.22=0.000000    B.22=0.000000
Harm. 23:    A.23=0.078125    B.23=0.269402
Harm. 24:    A.24=-0.000000    B.24=-0.000000
Harm. 25:    A.25=0.078125    B.25=0.246608
Harm. 26:    A.26=0.000000    B.26=-0.000000
Harm. 27:    A.27=0.078125    B.27=0.227093
Harm. 28:    A.28=-0.000000    B.28=0.000000
Harm. 29:    A.29=0.078125    B.29=0.210177
Harm. 30:    A.30=-0.000000    B.30=-0.000000
Harm. 31:    A.31=0.078125    B.31=0.195357
Harm. 32:    A.32=0.000000    B.32=0.000000
Harm. 33:    A.33=0.078125    B.33=0.182252
Harm. 34:    A.34=0.000000    B.34=-0.000000
Harm. 35:    A.35=0.078125    B.35=0.170566
Harm. 36:    A.36=0.000000    B.36=-0.000000
Harm. 37:    A.37=0.078125    B.37=0.160069
Harm. 38:    A.38=0.000000    B.38=0.000000
Harm. 39:    A.39=0.078125    B.39=0.150578
Harm. 40:    A.40=-0.000000    B.40=0.000000
Harm. 41:    A.41=0.078125    B.41=0.141944
Harm. 42:    A.42=-0.000000    B.42=0.000000
Harm. 43:    A.43=0.078125    B.43=0.134047
Harm. 44:    A.44=-0.000000    B.44=0.000000
Harm. 45:    A.45=0.078125    B.45=0.126789
Harm. 46:    A.46=-0.000000    B.46=0.000000
Harm. 47:    A.47=0.078125    B.47=0.120087
Harm. 48:    A.48=0.000000    B.48=0.000000
Harm. 49:    A.49=0.078125    B.49=0.113872
Harm. 50:    A.50=0.000000    B.50=-0.000000
Harm. 51:    A.51=0.078125    B.51=0.108087
Harm. 52:    A.52=-0.000000    B.52=0.000000
Harm. 53:    A.53=0.078125    B.53=0.102681
Harm. 54:    A.54=-0.000000    B.54=0.000000
Harm. 55:    A.55=0.078125    B.55=0.097614
Harm. 56:    A.56=0.000000    B.56=-0.000000
Harm. 57:    A.57=0.078125    B.57=0.092848
Harm. 58:    A.58=0.000000    B.58=-0.000000
Harm. 59:    A.59=0.078125    B.59=0.088353
Harm. 60:    A.60=0.000000    B.60=-0.000000
Harm. 61:    A.61=0.078125    B.61=0.084100
Harm. 62:    A.62=-0.000000    B.62=-0.000000
Harm. 63:    A.63=0.078125    B.63=0.080066
Harm. 64:    A.64=0.000000    B.64=-0.000000
Harm. 65:    A.65=0.078125    B.65=0.076231
Harm. 66:    A.66=0.000000    B.66=-0.000000
Harm. 67:    A.67=0.078125    B.67=0.072574
Harm. 68:    A.68=0.000000    B.68=0.000000
Harm. 69:    A.69=0.078125    B.69=0.069081
Harm. 70:    A.70=-0.000000    B.70=0.000000
Harm. 71:    A.71=0.078125    B.71=0.065736
Harm. 72:    A.72=-0.000000    B.72=0.000000
Harm. 73:    A.73=0.078125    B.73=0.062527
Harm. 74:    A.74=0.000000    B.74=0.000000
Harm. 75:    A.75=0.078125    B.75=0.059441
Harm. 76:    A.76=-0.000000    B.76=-0.000000
Harm. 77:    A.77=0.078125    B.77=0.056469
Harm. 78:    A.78=0.000000    B.78=-0.000000
Harm. 79:    A.79=0.078125    B.79=0.053600
Harm. 80:    A.80=0.000000    B.80=-0.000000
Harm. 81:    A.81=0.078125    B.81=0.050826
Harm. 82:    A.82=-0.000000    B.82=-0.000000
Harm. 83:    A.83=0.078125    B.83=0.048139
Harm. 84:    A.84=-0.000000    B.84=0.000000
Harm. 85:    A.85=0.078125    B.85=0.045533
Harm. 86:    A.86=0.000000    B.86=0.000000
Harm. 87:    A.87=0.078125    B.87=0.043000
Harm. 88:    A.88=-0.000000    B.88=0.000000
Harm. 89:    A.89=0.078125    B.89=0.040534
Harm. 90:    A.90=-0.000000    B.90=-0.000000
Harm. 91:    A.91=0.078125    B.91=0.038130
Harm. 92:    A.92=0.000000    B.92=-0.000000
Harm. 93:    A.93=0.078125    B.93=0.035784
Harm. 94:    A.94=-0.000000    B.94=0.000000
Harm. 95:    A.95=0.078125    B.95=0.033489
Harm. 96:    A.96=0.000000    B.96=0.000000
Harm. 97:    A.97=0.078125    B.97=0.031243
Harm. 98:    A.98=-0.000000    B.98=0.000000
Harm. 99:    A.99=0.078125    B.99=0.029040
Harm. 100:    A.100=-0.000000    B.100=-0.000000
Harm. 101:    A.101=0.078125    B.101=0.026877
Harm. 102:    A.102=0.000000    B.102=-0.000000
Harm. 103:    A.103=0.078125    B.103=0.024750
Harm. 104:    A.104=0.000000    B.104=0.000000
Harm. 105:    A.105=0.078125    B.105=0.022656
Harm. 106:    A.106=-0.000000    B.106=0.000000
Harm. 107:    A.107=0.078125    B.107=0.020591
Harm. 108:    A.108=-0.000000    B.108=-0.000000
Harm. 109:    A.109=0.078125    B.109=0.018553
Harm. 110:    A.110=-0.000000    B.110=-0.000000
Harm. 111:    A.111=0.078125    B.111=0.016539
Harm. 112:    A.112=0.000000    B.112=0.000000
Harm. 113:    A.113=0.078125    B.113=0.014546
Harm. 114:    A.114=0.000000    B.114=-0.000000
Harm. 115:    A.115=0.078125    B.115=0.012570
Harm. 116:    A.116=-0.000000    B.116=0.000000
Harm. 117:    A.117=0.078125    B.117=0.010611
Harm. 118:    A.118=0.000000    B.118=0.000000
Harm. 119:    A.119=0.078125    B.119=0.008664
Harm. 120:    A.120=0.000000    B.120=0.000000
Harm. 121:    A.121=0.078125    B.121=0.006728
Harm. 122:    A.122=0.000000    B.122=-0.000000
Harm. 123:    A.123=0.078125    B.123=0.004800
Harm. 124:    A.124=0.000000    B.124=-0.000000
Harm. 125:    A.125=0.078125    B.125=0.002878
Harm. 126:    A.126=0.000000    B.126=-0.000000
Harm. 127:    A.127=0.078125    B.127=0.000959

Pojawiają się niezerowe współczynniki szeregu trygonometrycznego, które z przewidywań teoretycznych powinny być zerowe. W wielu wypadkach wartość współczynnika an jest porównywalna z wartością współczynnika bn co znacząco wpływa na finalną amplitudę harmonicznej.
Soft (oczywiście nieparzysty).
Fourier_Series.7z
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.
Ostatnio zmieniony wtorek 25 sty 2022, 15:17 przez gaweł, łącznie zmieniany 1 raz.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » niedziela 16 sty 2022, 18:13

Przykład analizy przebiegów okresowych

Rozpatrzmy następującą funkcję okresową czasu: f(t)=-200 + 1000 sin ( t ) + 100 cos ( t ) + 500 sin ( 3 t ) + 300 cos ( 3 t ) (ot taka sobie fantazja). Jej przebieg pokazuje ilustracja 1.
sig27.png

Ten przebieg zostaje „spróbkowany” przez program w następujący sposób:

Kod: Zaznacz cały

void InputSampleTable ( void )
{
  unsigned short Loop ;   
  real Angle ;
  /*------------------------------------------------------------------------*/   
  for ( Loop = 0 ; Loop <= SampleNum ; Loop ++ )
  {
    Angle = DoublePi * ( real ) Loop / ( real ) SampleNum ;
    InpSampleTable[Loop]=100.0*(-2.0+10.0*sin(Angle)+cos(Angle)+
                         5.0*sin(3.0*Angle)+3.0*cos(3.0*Angle));
  } /* for */ ;
  T = 0.1 ; /* 100ms */
} /* InputSampleTable */

W wyniku przepuszczenia tych danych przez program został wyprodukowany następujący raport (w raporcie znajdują się dane M.<liczba> → jest to amplituda danej harmonicznej i alfa.<liczba> → jest to faza danej harmonicznej):

Kod: Zaznacz cały

Harm. 0:  A.0=-199.218750  B.0=   0.000000  M.0= 199.218750  alfa.0=-3.141593
Harm. 1:  A.1= 101.562500  B.1=1000.000000  M.1=1005.144239  alfa.1= 1.469581
Harm. 2:  A.2=   1.562500  B.2=  -0.000000  M.2=   1.562500  alfa.2=-0.000000
Harm. 3:  A.3= 301.562500  B.3= 500.000000  M.3= 583.900626  alfa.3= 1.028082
Harm. 4:  A.4=   1.562500  B.4=  -0.000000  M.4=   1.562500  alfa.4=-0.000000
Harm. 5:  A.5=   1.562500  B.5=  -0.000000  M.5=   1.562500  alfa.5=-0.000000
Harm. 6:  A.6=   1.562500  B.6=  -0.000000  M.6=   1.562500  alfa.6=-0.000000

W oparciu o dane można utworzyć następujące rysunki („trzymają” skalę). Zależność amplitudy od numeru harmonicznej pokazuje ilustracja 2. Przedstawia on widmo sygnału (ściślej widmo częstotliwościowe sygnału) – przedstawienie sygnału w dziedzinie częstotliwości lub pulsacji, otrzymane poprzez rozłożenie sygnału na poszczególne harmoniczne. Z wykresu tego można przykładowo odczytać, jakie składowe harmoniczne wchodzą w skład danego sygnału, czy sygnał ma ograniczone pasmo, jaka jest szerokość jego pasma itp.
Ponieważ każda harmoniczna jest charakteryzowana przez dwa współczynniki (an i bn), a wykres płaski pozwala na przedstawienie jednej informacji, to można pokazać amplitudy, jako
part07_form01.png

lub fazy (jako kąta utworzonego przez składowe harmonicznych, ilustracja 3).
sig28.png

sig30.png

W przypadku badanej funkcji, widmo fazowe pokazuje ilustracja 4.
sig31.png

Finalnie wyjściowa funkcja składa się z następujących elementów (ilustracja 5).
sig29.png
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
j23
Expert
Expert
Posty: 506
Rejestracja: czwartek 08 paź 2015, 18:40

Re: Cyfrowe przetwarzanie sygnałów

Postautor: j23 » wtorek 18 sty 2022, 01:04

Ja te wzory, zagadnienia i wykresy, prawdopodobnie jak większość rozumiem (co innego zastosować to w praktyce, pewnie byłoby dużo gorzej). W moim przypadku potrafię wyobrazić sobie to przestrzennie (jako tako, dopóki nie wchodzą całki potrójne, kołowe, etc. - bo wtedy trudniej), ALE... po tym co tutaj zobaczyłem, jeszcze te wszystkie przykłady i to w kodzie C, co ja mały, biedny jednak w wiedzę (w tym kontekście) mogę napisać.. Chylę czoła przed takim wspaniałym opracowaniem zagadnienia Kolego Gaweł. Dziękuję za ciekawe przybliżenie tematyki cyfrowego przetwarzania (w tym analizy) sygnałów. :like:
Pozdrawiam! J23
Internet łączy ludzi, którzy dzielą się swoimi zainteresowaniami, pomysłami i potrzebami, bez względu na geograficzne (przeciwności).
BOB TAYLOR, PARC

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » środa 19 sty 2022, 12:52

j23 pisze:dopóki nie wchodzą całki potrójne, kołowe, etc. - bo wtedy trudniej

Jak do tej pory to tylko był wstęp do wstępu. Mowa była o szeregach Fouriera a to nie jest fransformata Fouriera. By wdepnąć w DFT lub FFT to niestety musimy wdepnąć w liczby zespolone. Obiecuję, że całek wielokrotnych nie będzie :lol:

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » środa 19 sty 2022, 20:42

Rozkminy na temat okresu próbkowania

Opisany poprzednio przypadek rozkładu funkcji na harmoniczne jest w dużej mierze niejako akademicki → wszystko w nim się zgodziło, znaczy okres próbkowania. Weźmy pod uwagę jakiś „dobrze” znany przebieg, przykładowo następujący: 100 [ sin ( 3 t ) + cos ( 3 t ) ]. Odpowiednio „spróbkowany” w programie (dla odmiany w języku PASCAL):

Kod: Zaznacz cały

const
  SampleNumDiv2         = 128 ;
  SampleNum             = 2  * SampleNumDiv2 ;
  Pi                    = 3.141592654 ;
  DoublePi              = 6.283185307 ;

type

  SampleArrayType      = array [ 0 .. SampleNum ] of double ;
  CoofArrType          = array [ 0 .. SampleNumDiv2 ] of double ;

( . . . )

procedure TFourSeriesForm.InputSampleTable ;
var
  Loop                  : word ;
  Angle                 : double ;
  DAngle                : double ;
begin (* TFourSeriesForm.InputSampleTable *)
  DAngle := DoublePi / SampleNum ;
  for Loop := 0 to SampleNum do
  begin (* 1 *)
    Angle := DAngle * Loop ;
    InpSampleTable [ Loop ] := 100.0*(sin(3.0*Angle)+cos(3.0*Angle)) ;
  end (* 1 *) ;
  T := 0.1 ; (* 100ms *)
end (* TFourSeriesForm.InputSampleTable *) ;

Elementem istotnym w tym przypadku jest to, że okres próbkowania jest „za długi”. Pokazuje to ilustracja 1 (w obrębie okresu próbkowania w rzeczywistości wystąpiły trzy okresy funkcji sin oraz cos).
sig33.png

Takie dane z próbkowania zostały przetworzone przez program. Algorytm jest identyczny jak w przykładzie napisanym w języku C, sposób użycia jest również identyczny.

Kod: Zaznacz cały

procedure FourierSeries ( var ATable      : CoofArrType ;
                          var BTable      : CoofArrType ;
                          var SampleTable : SampleArrayType ) ;
var
  HarmInx               : word ;
  SampleInx             : word ;
  DAngle                : double ;
  Angle                 : double ;
begin (* FourierSeries *)
  for HarmInx := 0 to SampleNumDiv2 do
  begin (* 1 *)
    ATable [ HarmInx ] := 0.0 ;
    BTable [ HarmInx ] := 0.0 ;
  end (* 1 *) ;
  for SampleInx := 0 to SampleNum do
    ATable [ 0 ] := ATable [ 0 ] + SampleTable [ SampleInx ] ;
  for HarmInx := 1 to SampleNumDiv2 do
  begin (* 1 *)
    DAngle := DoublePi * HarmInx / SampleNum ;
    for SampleInx := 0 to SampleNum do
    begin (* 2 *)
      Angle := DAngle  * SampleInx ;
      ATable [ HarmInx ] := ATable [ HarmInx ] + SampleTable [ SampleInx ] * cos ( Angle ) ;
      BTable [ HarmInx ] := BTable [ HarmInx ] + SampleTable [ SampleInx ] * sin ( Angle ) ;
    end (* 2 *) ;
  end (* 1 *) ;
  ATable [ 0 ] := ATable [ 0 ] / SampleNum ;
  BTable [ 0 ] := 0.0 ;
  for HarmInx := 1 to SampleNumDiv2 do
  begin (* 1 *)
    ATable [ HarmInx ] := ATable [ HarmInx ] / SampleNumDiv2 ;
    BTable [ HarmInx ] := BTable [ HarmInx ] / SampleNumDiv2 ;
  end (* 1 *) ;
end (* FourierSeries *) ;     

Program wyrachował następujące elementy:

Kod: Zaznacz cały

Okres probkowania: 0,781 ms
Czestotliwosc probkowania: 1280,000 Hz
Czestotliwosc podstawowa (pierwsza harmoniczna): 5,000 Hz
Pulsacja podstawowa (pierwsza harmoniczna): 31,415927 rad/s
Harm. 0:  A.0=  0,3906249970  B.0=  0,0000000000  Czest. =0,000 Hz
Harm. 1:  A.1=  0,7812499941  B.1= -0,0000000001  Czest. =5,000 Hz
Harm. 2:  A.2=  0,7812499941  B.2= -0,0000000001  Czest. =10,000 Hz
Harm. 3:  A.3=100,7812499969  B.3=100,0000000026  Czest. =15,000 Hz
Harm. 4:  A.4=  0,7812499941  B.4= -0,0000000003  Czest. =20,000 Hz
Harm. 5:  A.5=  0,7812499941  B.5= -0,0000000004  Czest. =25,000 Hz

Analizując te wyniki powstaje sugestia, że program nie wykrył w tych próbkach harmonicznych o mniejszych numerach niż 3, no to by się zgadzało z ilustracją 1.
Wykres amplitud harmonicznych pokazuje ilustracja 2 oraz
sig34.png

wykres fazy owych harmonicznych pokazuje ilustracja 3.
sig35.png

Pamiętając o tożsamości trygonometrycznej
part08_form01.png

zrozumiałe staje się dlaczego faza sygnału wynosi 45º. Ogólny wniosek płynący z tego eksperymentu jest taki, że okres próbkowania nie odpowiada okresowi rzeczywistego sygnału.

Dla tropicieli, implementacja w Lazarus: (uwaliłem, by głupota się nie rozniosła)
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.
Ostatnio zmieniony piątek 21 sty 2022, 16:56 przez gaweł, łącznie zmieniany 2 razy.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
Zegar
User
User
Posty: 316
Rejestracja: wtorek 02 lip 2019, 14:42

Re: Cyfrowe przetwarzanie sygnałów

Postautor: Zegar » środa 19 sty 2022, 22:35

gaweł pisze:Dla tropicieli, implementacja w Lazarus:

Trudno coś wytropić...
Reputacja.png
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.
"If A = success, then the formula is A = X + Y + Z.
X is work. Y is play. Z is keep your mouth shut."
A. Einstein

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » środa 19 sty 2022, 23:04

Że niby fseries.exe ma złą reputację i musi być poddany kwarantannie? :lol:
Poniekąd temat na czasie :lol:
Jak uwaliło ci binarniaka, to mając źródła możesz sobie kompilnąć.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » piątek 21 sty 2022, 17:25

Najważniejsze jest próbkowanie

Analizując sygnały spróbkowane z określoną częstotliwością, przy niewłaściwym podejściu, można odnieść mylne przekonanie o otaczającej rzeczywistości. To tak, jakby popatrzeć na świat przez zakratowane okno (jakby ktoś przebywał w więzieniu) lub mieć widok przez płot ze sztachetami. Część widoku staje się niedostępna, przez co odnosi się mylny pogląd na otaczającą rzeczywistość (a ta często potrafi zaskoczyć). Rozpatrzmy układ próbkujący przebieg sinusoidalny o częstotliwości 25Hz z okresem próbkowania 10ms (czyli częstotliwość próbkowania wynosi 100Hz). Takie próbkowanie jest poprawne, gdyż zgodnie z twierdzeniem Kotielnikowa-Shannona, częstotliwość próbkowania jest minimum dwukrotnie większa od największej częstotliwości występującej w sygnale. Spróbkowany sygnał o częstotliwości 25Hz może być poprawnie odtworzony. Teraz zostaje zmieniona częstotliwość analizowanego sygnału z 25Hz na 125Hz. Okazuje się, że układ próbkujący nie zauważy żadnych różnic. Identycznie, jeżeli częstotliwość próbkowanego sygnału zostanie zmieniona na 225Hz, również układ próbkujący nie jest w stanie tego odróżnić. Oczywiście w obu tych przypadkach (125Hz i 225Hz) nie są spełnione warunki wymagane przez wspomniane przed chwilą twierdzenie, gdyż częstotliwość sygnału 125Hz musi być próbkowana minimum z częstotliwością 250Hz oraz 225Hz wymaga minimum 450Hz by jednoznacznie odtworzyć sygnał.
Przyczynę takiego stanu wyjaśnia ilustracja 1 (rysunek trzyma skalę).
sig38.png

Dla każdej z próbkowanych sinusoid, tak się jakoś dziwnie zbiegło, że ich wartości w chwili próbkowania są identyczne. Układ pomiarowy, którego zadaniem jest „obserwacja” zjawiska nie jest w stanie rozróżnić tych przypadków. Doskonale to widać na sumarycznym złożeniu wszystkich próbkowanych sygnałów. Gdyby można było podejrzeć, co się dzieje w międzyczasie, to można wyrobić sobie jakiś pogląd na zdarzenia ale to oznacza, że musi zostać podniesiona częstotliwość próbkowania.
Popatrzmy, co na to program. Zostało spróbkowanych pięć okresów funkcji sin o częstotliwości 25Hz, po 4 próbki na każdy okres, co dało finalnie 21 próbek (w sumie czas pomiaru jest dwa razy dłuższy od tego na ilustracji 1). Proces generacji próbek pokazuje poniższy kawałek kodu [sin(5.0*...) oznacza wzięcie pod uwagę 5 harmoniczną):

Kod: Zaznacz cały

procedure TFourSeriesForm.InputSampleTable ;
var
  Loop                  : word ;
  Angle                 : double ;
  DAngle                : double ;
begin (* TFourSeriesForm.InputSampleTable *)
  DAngle := DoublePi / SampleNum ;
  for Loop := 0 to SampleNum do
  begin (* 1 *)
    Angle := DAngle * Loop ;
    InpSampleTable [ Loop ] := sin ( 5.0 * Angle ) ;
  end (* 1 *) ;
  T := 0.2 ; (* Okres calego sprobkowanego sygnalu *)
end (* TFourSeriesForm.InputSampleTable *) ;

Wyniki programu są następujące:

Kod: Zaznacz cały

Probka [ 0 ] = 0,00000000
Probka [ 1 ] = 1,00000000
Probka [ 2 ] = 0,00000000
Probka [ 3 ] = -1,00000000
Probka [ 4 ] = 0,00000000
Probka [ 5 ] = 1,00000000
Probka [ 6 ] = 0,00000000
Probka [ 7 ] = -1,00000000
Probka [ 8 ] = 0,00000000
Probka [ 9 ] = 1,00000000
Probka [ 10 ] = 0,00000000
Probka [ 11 ] = -1,00000000
Probka [ 12 ] = 0,00000000
Probka [ 13 ] = 1,00000000
Probka [ 14 ] = 0,00000000
Probka [ 15 ] = -1,00000000
Probka [ 16 ] = 0,00000000
Probka [ 17 ] = 1,00000000
Probka [ 18 ] = 0,00000000
Probka [ 19 ] = -1,00000000
Probka [ 20 ] = 0,00000000
Liczba probek: 21
Okres calego sygnalu: 200,000 ms
*****************************************************************************
Okres probkowania: 10,000 ms
Czestotliwosc probkowania: 100,000 Hz
Czestotliwosc podstawowa (pierwsza harmoniczna): 5,000 Hz
Pulsacja podstawowa (pierwsza harmoniczna): 31,415927 rad/s
Harm. 0:    A.0=0,0000000000    B.0=0,0000000000   Czest. =0,000 Hz
Harm. 1:    A.1=0,0000000000    B.1=0,0000000000   Czest. =5,000 Hz
Harm. 2:    A.2=0,0000000000    B.2=0,0000000000   Czest. =10,000 Hz
Harm. 3:    A.3=0,0000000000    B.3=0,0000000000   Czest. =15,000 Hz
Harm. 4:    A.4=0,0000000000    B.4=0,0000000000   Czest. =20,000 Hz
Harm. 5:    A.5=0,0000000000    B.5=1,0000000000   Czest. =25,000 Hz
Harm. 6:    A.6=0,0000000000    B.6=0,0000000000   Czest. =30,000 Hz
Harm. 7:    A.7=0,0000000000    B.7=0,0000000000   Czest. =35,000 Hz
Harm. 8:    A.8=0,0000000000    B.8=-0,0000000001  Czest. =40,000 Hz
Harm. 9:    A.9=0,0000000000    B.9=-0,0000000001  Czest. =45,000 Hz
Harm. 10:  A.10=0,0000000000   B.10=-0,0000000001  Czest. =50,000 Hz

Patrząc na dane pomiarowe, to pokrywają się one z rysunkiem (przyjmowane są wartości 0, 1 i -1). Program oczywiście zakłada, że cała paczka pomiarowa dotyczy jednego okresu, który trwa 200ms. Odpowiada to częstotliwości podstawowej (czyli pierwszej harmonicznej) wynoszącej 5Hz. Oczywiście takiej tam nie ma, ale zostaje stwierdzona harmoniczna numer pięć o amplitudzie 1 (odpowiadająca częstotliwości 25Hz).
Jeżeli teraz zostanie spróbkowany sygnał o częstotliwości 125Hz lub 225Hz, to co się zmieni? Każdy może sam wyciągnąć właściwy wniosek → nic się nie zmieni. Program nadal wykryje jedynie częstotliwość 25Hz. Innych nie jest w stanie dostrzec.
Może warto to sprawdzić zamiast brać coś na wiarę (bywa, że można czasami się pomylić). Generacja danych pomiarowych (25 w funkcji sin oznacza 25-tą harmoniczną o wartości 125Hz):

Kod: Zaznacz cały

procedure TFourSeriesForm.InputSampleTable ;
var
  Loop                  : word ;
  Angle                 : double ;
  DAngle                : double ;
begin (* TFourSeriesForm.InputSampleTable *)
  DAngle := DoublePi / SampleNum ;
  for Loop := 0 to SampleNum do
  begin (* 1 *)
    Angle := DAngle * Loop ;
    InpSampleTable [ Loop ] := sin ( 25.0 * Angle ) ;
  end (* 1 *) ;
  T := 0.2 ; (* Okres calego sprobkowanego sygnalu *)
end (* TFourSeriesForm.InputSampleTable *) ;

Po uruchomieniu programu wyszło dokładnie to samo.
Można w tym miejscu zadać proste pytanie (pytanie jest w rzeczy samej proste, gorzej jest ze znalezieniem prostej odpowiedzi): dlaczego tak się dzieje jak się dzieje?
By uzyskać prostą odpowiedź, koniecznej jest zagłębienie się w liczby zespolone i ten temat jeszcze powróci. Teraz można przyjąć następującą odpowiedź:
jeżeli jest próbkowany sygnał o częstotliwości fs i częstotliwość próbkowania wynosi fp, to nierozróżnialne są częstotliwości sygnału spełniające zależność fs + k * fp, gdzie k jest dowolną liczną naturalną.

PS
Stało się :( , wykryłem wadę w programie (trochę prowadzi błędne obliczenia). Robić pomyłki jest rzeczą ludzką, jak również je naprawiać. Cóż, przyznaję się i biorę na klatę...(walnąłem się w pierś). Więcej błędów nie zauważyłem :D .

Dla tropicieli aktualny soft:
fourseries.7z
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » niedziela 23 sty 2022, 01:30

Meandry próbkowania

Zawsze można się zastanawiać, na ile sposobów można coś zepsuć lub prawie zepsuć. Niewłaściwe próbkowanie jest właśnie takim psuciem. Poprzedni przykład pokazywał, że coś można ukryć, próbkowany przebieg sygnału stworzył prawdopodobną, jednak fałszywą iluzję samego siebie. Wszystko jest kwestią punktu widzenia i doboru odpowiednich parametrów w całym tym procesie. Przedstawię inny ciekawy wariant „fałszowania rzeczywistości”, gdzie występuje sygnał o częstotliwości 48Hz. Sygnał ten jest próbkowany z okresem 10ms (jak łatwo wyrachować, odpowiada to częstotliwości próbkowania 100Hz). Jeżeli komuś przyszło do głowy by sprawdzić czy spełnione są kryteria poprawnego próbkowania, to jest na dobrej drodze. Parametry próbkowania spełniają wymogi: sygnał o częstotliwości 48Hz wymaga minimum 96Hz częstotliwości próbkowania (a jest 100Hz). Natomiast ciekawa jest interakcja sygnału mierzonego z układem pomiarowym, pokazuje to ilustracja 1.
sig36.png

To co ujrzy światło dzienne w tym układzie pomiarowym pokazuje wyekstrahowanych przebieg w kolorze czerwonym. Widok ten powinien być dobrze znajomy wszystkim użytkownikom oscyloskopów i oznacza, że częstotliwość próbkowania jest trochę za mała w stosunku do częstotliwości sygnału mierzonego.
W tym przypadku wystarczy zwiększyć częstotliwość próbkowania (w oscylku rozciągnąć podstawę czasu) i wszystko wskakuje na lepsze tory. Zwiększenie częstotliwości próbkowania do 1kHz znacząco pomaga, sygnał ma około 20 próbek na swój jeden okres (ilustracja 2).
sig37.png

Podsumowując już te wszystkie rozkminy, to naturalnym jest, że nasuwa się wniosek, który sprowadza się do prostego stwierdzenia: parametry próbkowania sygnału muszą spełniać pewne wymogi i właściwie powinny być indywidualnie dobierane do sytuacji. Tu łatwo „skręcić w krzaki” i utknąć na walce z eliminacją pewnych efektów ubocznych.
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » poniedziałek 24 sty 2022, 16:40

Pasmo przenoszenia

Kiedyś w historii zaintrygował mnie pewien efekt, jaki daje się zaobserwować korzystając z oscyloskopu. Zdarzało się śledząc sygnał prostokątny zobaczyć pewne zafalowania występujące na „kantach” sygnału prostokątnego (ilustracja 1).
sig39.png

Przyszło mi wtedy do głowy, że może sprzęt czasami sobie nie radzi z gwałtownymi zmianami. Jednak to by jedynie tłumaczyło wystąpienie efektu „przeregulowania” jedynie dla zbocza narastającego. Wystąpiła gwałtowna zmiana, więc stała się przyczyną jakichś tam efektów. No ale ta koncepcja nie pasuje do zbocza opadającego. Skąd oscylek wie, że należy zacząć „falować” bo wkrótce wystąpi gwałtowna zmiana (zbocze opadające sygnału). To tak jakby przewidywał przyszłość?
Ta swoista prekognicja ma solidne podstawy w matematyce i fizyce i nosi nazwę efektu Gibbsa. Jak sobie przypomnimy z wcześniejszych rozważań, falę prostokątną można rozwinąć w nieskończony szereg trygonometryczny:
part11_form01.png

By ten sygnał został przeniesiony przez sprzęt w niezmienionej formie, to sprzęt musi mieć nieskończenie wielkie pasmo przenoszenia, przenosić (nie tłumić) wszystkie częstotliwości (a te w swych wartościach dążą do nieskończoności). To trochę mało realizowalne w rzeczywistości, od pewnego miejsca należy sobie uzmysłowić, że maszyneria więcej nie może i to w jakiś sposób odciska swoje piętno na tym, co daje się zaobserwować. Jeżeli przykładowo nasz oscylek ma pasmo do 50MHz, to nie znaczy, że sygnał o częstotliwości 51MHz już nie przejdzie. Przejdzie, tylko zostanie lekko zniekształcony i każda następna harmoniczna również.
Przykładowo: jeżeli mierzony sygnał ma 8MHz a oscylek pasmo przenoszenia 50MHz, to obserwowany efekt może być następujący (ilustracja 2). Przy tych założeniach, sprzęt przenosi bez zniekształceń harmoniczne o numerach 0..5 (w tym konkretnym wypadku jako nieparzyste 1, 3 i 5). Każde następne ulegają jakiejś postępującej degradacji i tracą na wartości amplitudy aż do sytuacji, gdzie przestają być przenoszone. Tu dla uproszczenie nie biorę pod uwagę efektów wynikających z przesunięć fazowych a jedynie wpływ na amplitudę. W rzeczywistości trzeba liczyć się z każdą możliwością i wszystko może się zdarzyć.
sig40.png

Całe zajście można rozpatrywać w kategorii filtru dolnoprzepustowego. W realnym sprzęcie, od jakiegoś momentu harmoniczne są tłumione w określony sposób (jednak odczuwalne jest ich echo w finalnym sygnale, w teraźniejszości). W świecie cyfrowego przetwarzania sygnału wszystko może okazać się możliwe, toteż nikt nie zabroni, by od harmonicznej o określonej częstotliwości pozbyć się ich wpływu, po prostu się odciąć.
sig41.png

Na efekty nie trzeba długo czekać (ilustracja 3), bo w świecie DSP wszystko jest możliwe.
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » środa 26 sty 2022, 01:00

Ulepszona wersja programu

Przyznam szczerze, że dotychczasowe programy nie były najlepsze (w kwestii interfejsu użytkownika). Procedura generująca próbki pomiarowe miała tą specyfikę, że generowała dane jedynie w obrębie jednego okresu funkcji. Do tej pory to wystarczało, gdyż rozważania skupiały się jedynie na pojedynczym okresie próbkowanej funkcji. Jak to zwykle bywa, po jakimś czasie przestaje to wystarczać i zachodzi potrzeba „więcej”. Właśnie do czegoś takiego doszło, spojrzenia ze znacząco szerszej perspektywy, trzeba radzić sobie z bardziej „rozmytym” środowiskiem, jak przykładowo spróbkowanie funkcji o nieokreślonym z góry okresie. Tak jak w prawdziwej rzeczywistości, robimy pomiary i przymiarki i nikt nie da gwarancji, że wszystko będzie dopięte na ostatni guzik. Dopiero z jakiejś perspektywy można dostrzec tematykę w szerszym zakresie. Podam przykład. Niech spróbkowany sygnał dotyczy EKG (określa pracę serca). Ile jest niezbędnych próbek, by coś można wywnioskować? A liczba uderzeń na minutę? Czy trzeba próbkować przez minutę? A może wystarczy jakaś aproksymacją po kilku uderzeniach? Takich pytań można postawić znacząco więcej.
Dotychczas program dzielił jeden okres próbkowanej funkcji na ileś interwałów czasu, co prowadzi do różnych okresów próbkowania (jak i również różnych częstotliwości próbkowania) i daje w sumie jeden okres funkcji. Obecnie wprowadziłem zmianę, że program ma stały okres próbkowania i „jedzie” w czasie aż skończy się miejsce w tablicy próbek. To oznacza, że może zostać spróbkowana funkcja, która zawiera w sobie wiele (nie koniecznie całkowitych) okresów funkcji.
Z tego powodu funkcja przewidziana do generowania próbek dostaje w parametrze okres próbkowania (wyrażony w sekundach). Chcąc teraz wygenerować określony przebieg należy to odpowiednio ująć w „pobraniu próbki”. Przykładowo:

Kod: Zaznacz cały

InpSampleTable [ Loop ] := 10.0 + 5.0 * sin ( DoublePi * 15.0 * Time ) +
                           7.0 * cos ( DoublePi * 25.0 * Time ) ;

generuje 15Hz i 25Hz.
Jednocześnie program generuje kilka plików danych dla programu OCTAVE do robienia ładnych rysunków. W zakresie rozwiązań numerycznych, program nie zawiera zmian.
Dla powyższego spróbkowanego przebiegu powstają następujące rysunki przebiegów wybranych sygnałów.
Próbkowany sygnał w funkcji numeru próbki (ilustracja 1).
sig42.png

Próbkowany sygnał w funkcji czasu (ilustracja 2).
sig43.png

Rozkład harmonicznych w funkcji częstotliwości harmonicznych (ilustracja 3).
sig44.png


Dla tropicieli:
fourseries.7z
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » środa 26 sty 2022, 01:09

Oczywiście można sobie zadać pytanie: no fajnie, ale po jaką chusteczkę to wszystko?
Ano mając znane częstotliwości, to można wyznaczyć transformatę Laplace'a dla sygnału i finalnie wyznaczyć transmitancję układu. No i tu się zaczyna już fajna zabawa, tylko muszę się udać na swoją uczelnię i wydobyć własną pracę dyplomową. Tam było kilka całkiem fajnych wynalazków...
Co prawda, trzeba się jeszcze trochę napracować nad tym, bo ... wszystko było napisane w FORTRAN'ie. To taki prawie starożytny język programowania, ale w moim przekonaniu nadal niedościgły w obliczeniach numerycznych. Może być ciekawie, zapraszam chętnych do wspólnej zabawy.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » sobota 29 sty 2022, 01:45

No więc udałem się na swoją uczelnię i spotkała mnie przykra niespodzianka: uczelnia zgubiła część mojej pracy :( . Mogłem sobie pooglądać "książkę", która w gruncie rzeczy jest najmniej istotnym elementem.
sig_t.png

Coś co było najbardziej cenne, czyli załącznik w postaci wydruku programu umyka. To tam były numeryczne rozwiązania przykładowo transformaty Laplace'a. No generuje to jakieś dodatkowe problemy, bo wymaga dodatkowej pracy by te coś odzyskać. Trudno, ale da się to zrobić.
Od tego ma się kumpli na uczelni, którzy podrzucą parę rozwiązań numerycznych.
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » sobota 29 sty 2022, 02:01

Transformata Fouriera

Nadszedł czas, by przejść do transformaty Fouriera. To, co do tej pory było przedstawione, nie jest w żadnym stopniu transformatą Fouriera (nawet nigdzie nie zostało użyte takie sformułowanie, mowa była o szeregach trygonometrycznych, które są tożsame z szeregami Fouriera). Ważnym jest, by uzmysłowić sobie istotne elementy z tym związane. Szereg trygonometryczny w punktu widzenia matematycznego jest szeregiem funkcji (finalnie nadal funkcją) zmiennej rzeczywistej (czasu) o wartościach rzeczywistych. Transformata Fouriera jest przekształceniem (całkowym), które przenosi funkcję czasu (zmiennej rzeczywistej) zwanej oryginałem do przestrzeni zespolonej i jej zadaniem jest „uwypuklenie” pewnych cech funkcji oryginału. Te „wynalazki” są autorstwa monsieur Fourier. Jest on ojcem tego, ale te dzieci to niezależne byty. Mają wiele wspólnego, ale to są różne rzeczy. Wspomniana przestrzeń zespolona jest zbiorem liczb zespolonych. To takie dziwne liczby, które zamiast być umiejscowione na osi liczbowej, są umiejscowione na płaszczyźnie. Te liczby w ujęciu „elektrycznym” mają troszkę odmienną symbolikę niż dla „reszty świata”, chodzi o oznaczenie jednostki urojonej. Elektrycy na jej oznaczenie używają literki j, natomiast wszyscy inni (szczególnie matematycy) używają literki i. Ponieważ uważam się za elektryka, w swoich rozważaniach stosuję oznaczenie j.
Zanim przejdziemy do matematyki, trochę filozoficznych rozważań. Wyobraźmy sobie źródło białego światła (białego, bo zawiera w sobie każdą barwę). To światło przechodzi przez jakąś kolorową szybkę (ilustracja 1). Powstaje pytanie: co z tego wyszło?
sig50.png

To dosyć proste pytanie generuje parę trudnych odpowiedzi. Po pierwsze co to znaczy „co wyszło”, jak to zmierzyć i tak naprawdę co jest mierzone. W przypadku światła to można w miarę prosto odpowiedzieć. Światło jest mieszanką wszystkiego (wszystkich kolorów), które są identyfikowane przez swoją barwę czyli częstotliwość (lub pulsację, bo zagadnienie jest proporcjonalne). Po przejściu przez jakąś barierę doznaje swoistej przemiany, niektóre barwy zanikają, zawartość barwowa jest inna. Można w miarę prosto zmierzyć poziom sygnału (jako przykładowo procent w stosunku do tego co zostało wyemitowane). Jednak my chcemy więcej, interesuje nas ile jest każdej barwy. Naturalnym rozwiązaniem nasuwającym się każdemu jest: przepuścić światło przez pryzmat (ilustracja 2).
sig51.png

Podobnie działa transformacja Fouriera: rozkłada wejściową funkcję na jakieś składowe, pozwala określić z jakich elementów się składa. W przypadku światła z jakich kolorów się składa (no i oczywiście ile tego jest w danym kolorze).
Posłużę się jeszcze jedną analogią. Mamy w szufladzie dużo kolorowych klocków. Pytanie jest następujące: ile jest tych detali w każdym kolorze? Bierzemy monochromatyczną lupę, przez którą widać jedynie jedną barwę. Patrząc z tej perspektywy można policzyć widoczne detale (tylko te, które leżą w obszarze zainteresowania, pozostałych nie widać). Po zmianie lupy na inny kolor można powtórzyć operację. Teraz wyobraźmy sobie, że lupa ma w sobie takie pokrętełko, która zmienia barwę w sposób analogowy. Można w sposób ciągły przeskanować całe spektrum barwowe.
Wracając do koncepcji związanej ze światłem padającym na pryzmat, to niech przesłona będzie zmienna (jak obracające się kółko zmieniające skład kolorów) i niech lampka będąca źródłem światła będzie zmieniała jego natężenie (ilustracja 3).
sig52.png

Dochodzimy do koncepcyjnej istoty działania transformaty Fouriera. By określić ile jest czegoś w całości, to należy podzielić całość przez te „coś”. Rolę całości spełnia funkcja f(t), natomiast rolę tego „czegoś”, takiej lupy, pełni e. Ile jakiejś ω zawiera się w f(t) określa wyrażenie:
part13_from01.png

Oczywiście to tylko pokazuje tą „magiczną” zawartość tylko w jednym konkretnym miejscu na osi czasu dla funkcji oryginału. By uzyskać „komplet” danych należy zsumować te wszystkie przyczynki, do których przykładany jest filtr zbudowany z lupy (ilustracja 4), czyli przeskanować całą oś liczbową argumentu funkcji oryginału.
sig53.png

Sumowanie tak małych elementów (taki monochromatyczny prążek jest nieskończenie cienki) wymaga zastosowania całki. Tu warto pamiętać, że granice całkowania dotyczą argumentu funkcji oryginału, czyli dotyczą czasu.
part13_from02.png

Powyższa formuła stanowi definicję całkowego przekształcenia Fouriera funkcji czasu (oryginału) f(t) do funkcji zespolonej F(jω) popularnie zwanej transformatą Fouriera. Jest to już funkcja zespolona. Tu warto zdać sobie sprawę z kolejnej istotnej różnicy między szeregiem Fouriera a transformatę Fouriera. Funkcja czasu nie musi być okresowa i może rozciągać się od minus nieskończoności do plus nieskończoności. Można na to spojrzeć z takiej perspektywy: funkcja oryginału jest okresowa, tylko jej okres jest nieskończenie duży.

Robienie wała z funkcji, czyli nawijanie na wałek

Element formuły e-jω daje się zapisać zgodnie z trygonometryczną reprezentacją liczb zespolonych jako cos(ω) - j sin(ω) jest operatorem obrotu na płaszczyźnie zespolonej (znak przed jω określa kierunek obrotu). Przykładowo jeżeli:
part13_from03.png

to mnożenie czegoś przez j realizuje na płaszczyźnie zespolonej obrót o kąt 90°, pokazuje to ilustracja 5.
sig54.png

Zastosowanie tej technologii na funkcję oryginału, czyli nawijanie funkcji na zespolone koło pokazuje ilustracja 6. Oczywiście zawsze wystąpi problem dopasowania długości (długości dziedziny funkcji do długości okręgu zespolonego). Rolę dopasowywacza pełni właśnie pulsacja ω.
sig55.png
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » sobota 29 sty 2022, 20:41

Dyskretna transformata Fouriera

Dyskretna transformata Fouriera to pewna odmiana klasycznej transformaty Fouriera. Ta odmienność polega przede wszystkim na dyskretyzacji całego zjawiska i nie sprowadza się to jedynie do dyskretyzacji sygnału funkcji oryginału, czyli sprowadzeniu funkcji do ciągu próbek. Typowo tę transformatę określa się jako DFT (od ang. Discrete Fourier Transform).
Klasyczna transformata Fouriera jest zjawiskiem ciągłym, zarówno czas w funkcji oryginału jak i wyjściowe pulsacje są sygnałem ciągłym, analogowym. Przy próbkowaniu sygnału powstaje ciąg próbek, który jest przetwarzany przez transformatę dając również ciąg. W przeciwieństwie do szeregu trygonometrycznego wynikiem transformaty nie są amplitudy kolejnych harmonicznych. Można uzyskać z nich amplitudy, ale to wymaga dodatkowych rachunków. Można zadać sobie pytanie: co jest wynikiem transformaty DFT?
DFT jest chwytem, procedurą matematyczną używaną do wyznaczenia zawartości harmonicznych. Tak jak sygnał wejściowy określa próbki w dziedzinie czasu, tak wynik DFT określa próbki w dziedzinie częstotliwości. W klasycznym (ciągłym) ujęciu transformaty rezultat jest funkcją ciągłą. W ujęciu dyskretnym, rezultat jest ciągiem próbek. Próbkując sygnał wejściowy może się zdarzyć, że źle trafiliśmy z próbkowaniem. To skutkuje również złym trafieniem w rozkład widma i zamiast fajnych wskazań na częstotliwości harmoniczne można uzyskać rozmyte dane. To zjawisko nawet ma swoją nazwę: wycieki widma.
Generalnie temat jest bardzo szeroki i złożony. Trudno jest posiąść wiedzę „za jednym zasiadem” i raczej sprowadza się to do ciągłego doświadczania zjawiska. Najlepszym rozwiązaniem jest metoda „drobnych kroczków”. Wystarczy zrobić coś niewielkiego, sprawdzić, zbadać i wyciągnąć wnioski. Dobrym przybliżeniem zjawiska jest „rolowanie funkcji”. Takie rolowanie to dosyć efektywny (i efektowny) sposób na przekonanie się, że ta rzeczywiście działa. Można sprawę odwrócić i realizować działania w dziedzinie czasu, ale wtedy trudno jest ocenić, że wszystko się powtarza. Nakładanie obrazu na samego siebie pozwala bardzo precyzyjnie dostrzec okresowość.
Weźmy przykładowo sygnał, którego kilka okresów pokazuje ilustracja 1.
sig56.png

Będziemy ten sygnał nawijać na zespolone koło z prędkością ω=0.1 (ilustracja 2), czyli rozciągać lub zagęszczać czas z osi czasu w funkcji oryginału. Z rysunku można wywnioskować, że sygnał ma coś wspólnego z pulsacją ω=0.1 rad/s. Przeskalowanie osi czasu z tym współczynnikiem daje nakładający się na siebie obraz.
sig57.png

Odmiennie sprawa wygląda jeżeli wartości współczynników lekko będą się różnić od wyżej wymienionego.
sig58.png

Z racji, że sygnał na ilustracji 1 ma swoją „długość”, co daje 1,5 lub prawie 2 obroty na zespolonym kole (ilustracja 3). Gdyby sygnał był dłuższy, to wygeneruje bardziej lub mniej „rozmazaną oponkę”.
Jeżeli sygnał w dziedzinie czasu zostanie przeskalowany przez ω=2.0 (ilustracja 4), to znów ujawni swoje powinowactwo.
sig59.png

Niewielkie „zejście z kursu” daje wyraźny sygnał, że coś jest nie tak. Pokazuje to ilustracja 5.
sig60.png

Podobnie dla ω=3.0 (ilustracja 6 i 7).
sig61.png

sig62.png

Stworzyłem sobie program do badania tego zjawiska. Można samodzielnie poeksperymentować i pobadać to zjawisko. Przykładowo: co ma wspólnego pulsacja ω=3 oraz ω=1/3 (w kolejnych przybliżeniach dziesiętnych)?
sshot14_01.png

sshot14_02.png

Kolejne rysunki pokazują, jak subtelne jest to zjawisko.
sshot14_03.png


Program do rolowania:
rollover.7z
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse

Awatar użytkownika
gaweł
Geek
Geek
Posty: 1259
Rejestracja: wtorek 24 sty 2017, 22:05
Lokalizacja: Białystok

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » wtorek 01 lut 2022, 00:09

Czas na rachunki

Najwyższy czas na rozpoczęcie prawdziwej pracy związanej z transformatą Fouriera. Wspomniałem, że funkcja oryginału może rozciągać się od -∞ do +∞. Ta fajna koncepcja ma jednak jedną poważną wadę: jest nierealizowalna fizycznie, tzn. nie da się stworzyć takiego programu, który będzie rachował w wyżej wymienionych obszarach. By całość zagadnienia urealnić, ludzkość wymyśliła coś pozwalające na pokonanie owych problemów. Jak mówi pewne porzekadło: na każdą koncepcję zawsze da się zaproponować antykoncepcję. W przypadku przetwarzania sygnałów taką rolę pełni okno do próbkowania. Jako okno rozumiana jest funkcja, której zadaniem jest wycięcie z ciągu próbek jakiegoś jej fragmentu. Użyte sformułowanie „ciągu próbek” sugeruje właśnie zdyskretyzowaną funkcję. Biorąc pod uwagę w miarę proste rozwiązanie funkcji okienkującej, takiej, że w przedziale <0,T) funkcja daje wartość jeden oraz wartość zero poza tym przedziałem. Wymnożenie funkcji oryginału przez funkcję okienkującą da oczekiwany ciąg próbek. Ideę pokazuje ilustracja 1.
sig65.png

Wyokienkowaną funkcje można już bez problemów próbkować, nie jest to zagadnienie trwające nieskończenie długo. Poza „okienkiem” funkcja przybiera wartości zerowe, a transformata Fouriera dla nich również ma wartość zerową. Z punktu widzenia matematyki transformata Fouriera zmieniła się do postaci:
part15_from01.png

Dla przypadku dyskretnego operacja całkowania zamienia się na sumę. Przed rozpoczęciem rachunków konieczne staje się uświadomienie sobie pewnej własności. W przypadku szeregów trygonometrycznych zachodziła potrzeba uzyskania nieparzystej liczby próbek. W nowej koncepcji to wymaganie zmienia się na: parzysta liczba próbek. Powodem tego jest to, że harmoniczna zerowa (wartość stała funkcji) zmieniła swój charakter z rzeczywistej na zespoloną, a ta wymaga dwóch wartości by jednoznacznie ją określić. Implementacja transformaty DFT wynika z jej definicji: będziemy kręcić ciągiem próbek na wyimaginowanym kole. Ponieważ w języku PASCAL nie ma standardowego typu do rachunków zespolonych niezbędne staje się dorobienie paru detali:

Kod: Zaznacz cały

const
  Pi               = 3.141592654 ;
  DoublePi         = 6.283185307 ;

type
  ComplexT           = record
                         Re           : double ;
                         Im           : double ;
                       end (* record *) ;

Bazując na tych podstawach przydatne staje się powołanie typów do reprezentacji widma zespolonego oraz samego spróbkowanego sygnału (warto zauważyć parzystość liczby elementów obu tablic):

Kod: Zaznacz cały

const
  SampleNumDiv2         = 256 ;
  SampleNum             = 2  * SampleNumDiv2 ;

type
  SampleArrayType       = array [ 0 .. SampleNum - 1 ] of double ;
  HarmArrayTable        = array [ 0 .. SampleNum - 1 ] of ComplexT ;

Sama procedura do rachowania DFT jest następująca (to mój wariant, bo wariantów jest tyle i jest autorów):

Kod: Zaznacz cały

procedure DFT ( var DFTTable      : array of ComplexT ;
                var SampleTable   : array of double ;
                    SampleTabSize : word  ) ;
var
  HarmInx               : word ;
  SampleInx             : word ;
  DAngle                : double ;
  Angle                 : double ;
begin (* DFT *)
  SampleTabSize := SampleTabSize - 1 ;
  for HarmInx := 0 to SampleTabSize do
  begin (* 1 *)
    DFTTable [ HarmInx ] := CAssign ( 0.0 , 0.0 ) ;
  end (* 1 *) ;
  for HarmInx := 0 to SampleTabSize do
  begin (* 1 *)
    DAngle := DoublePi * HarmInx / SampleNum ;
    for SampleInx := 0 to SampleTabSize do
    begin (* 2 *)
      Angle := DAngle  * SampleInx ;
      DFTTable[HarmInx].Re := DFTTable[HarmInx].Re + SampleTable[SampleInx]*cos(Angle) ;
      DFTTable[HarmInx].Im := DFTTable[HarmInx].Im - SampleTable[SampleInx]*sin(Angle) ;
    end (* 2 *) ;
  end (* 1 *) ;
end (* DFT *) ; 

Monsieur Fourier oprócz transformacji funkcji czasu do przestrzeni zespolonej, zaproponował również operację odwrotną pozwalającą na powrót do normalności. Ta formuła dla klasycznej (ciągłej) operacji to:
part15_from02.png

W wersji komputerowej następuje kolejna zamiana całki na sumę. Jak się przyjrzeć, to odwrotna transformata Fouriera daje się implementować identycznym algorytmem. Różnica jest taka, że w jedną stroną jest łatwiej (nie występuje czynnik skalujący przed całką). Przy powrocie już jest trudniej, występuje jakiś czynnik, który należy uwzględnić (ciągną się zaszłości z transformaty „tam”). Tu po raz kolejny można się przekonać, że każdy robi po swojemu. Często można spotkać wariant, że „tam” i „z powrotem” ma przed sobą czynnik skalujący typu
part15_from03.png

co pozwala na wykorzystanie dokładnie tego samego kodu procedury do rachunków w obie strony (temat jeszcze powróci). Na razie są dwa warianty procedur, gdyż jest znacząco odmienna reprezentacja danych wejściowych i wyjściowych.

Kod: Zaznacz cały

procedure IDFT ( var SampleTable   : array of double ;
                 var DFTTable      : array of ComplexT ;
                     SampleTabSize : word  ) ;
var
  HarmInx               : word ;
  SampleInx             : word ;
  DAngle                : double ;
  Angle                 : double ;
begin (* IDFT *)
  SampleTabSize := SampleTabSize - 1 ;
  for SampleInx := 0 to SampleTabSize do
  begin (* 1 *)
    SampleTable [ SampleInx ] := 0.0 ;
  end (* 1 *) ;
  for SampleInx := 0 to SampleTabSize do
  begin (* 1 *)
    DAngle := - DoublePi * SampleInx / SampleNum ;
    for HarmInx := 0 to SampleTabSize do
    begin (* 2 *)
      Angle := DAngle  * HarmInx ;
      SampleTable [ SampleInx ] := SampleTable [ SampleInx ] +
                                   DFTTable [ HarmInx ] . Re * cos ( Angle ) +
                                   DFTTable [ HarmInx ] . Im * sin ( Angle ) ;
    end (* 2 *) ;
    SampleTable [ SampleInx ] := SampleTable [ SampleInx ] / SampleTabSize ;
  end (* 1 *) ;
end (* IDFT *) ;


Nie ma jak własne testy

Można się samemu zagłębić w zagadnienie. Została wygenerowana seria próbek w następujący sposób (zawiera dwie częstotliwości):

Kod: Zaznacz cały

procedure TDFTForm.InputSampleTable ( SampleInterv : double ) ;

var
  Loop                  : word ;
  Time                  : double ;
begin (* TDFTForm.InputSampleTable *)
  for Loop := 0 to SampleNum - 1 do
  begin (* 1 *)
    Time := SampleInterv * Loop ;
    TimeTable [ Loop ] := Time ;
    InpSampleTable [ Loop ] := 10.0 + 5.0 * sin ( DoublePi * 15.0 * Time ) +
                                      7.0 * cos ( DoublePi * 25.0 * Time ) ;
  end (* 1 *) ;
  TotalTime := SampleNum * SampleInterv ;
end (* TDFTForm.InputSampleTable *) ;

gdzie parametr wywołania funkcji określa okres próbkowania (użyta wartość jest 0.002).
Generuje to następujące dane:

Kod: Zaznacz cały

Liczba probek= 512
( . . )
Okres calego sygnalu: 1024,000 ms
*****************************************************************************
Okres probkowania: 2,000000 ms
Czestotliwosc probkowania: 500,000 Hz
Czestotliwosc podstawowa (pierwsza harmoniczna): 0,977 Hz
Pulsacja podstawowa (pierwsza harmoniczna): 6,135923 rad/s

Skoro zastosowana w danych częstotliwość próbkowania wynosi 500Hz, to częstotliwość Nyquista ma wartość 250Hz (jest dwa razy mniejsza).
Po wygenerowaniu danych, program obliczył transformatę „tam” oraz „z powrotem”, porównał wartości jakie zostały wygenerowane z wartościami wynikającymi z podróży „tam” i „z powrotem”. Nie jest dziwne, że nieznacznie się różnią (reprezentacja liczb w komputerze również nie jest ciągła i ma swoją „ziarnistość”).
Na wejściu jest sygnał (ilustracja 2 w funkcji numeru próbki lub ilustracja 3 w funkcji czasu):
sig70.png

sig71.png

Kolejna różnica w stosunku do szeregu trygonometrycznego jest taka, że N próbek sygnału daje N próbek częstotliwości (w szeregu trygonometrycznym N próbek generowało N/2 harmonicznych). Występuje tu swoista redundancja wyników, część danych jest identyczna i właściwie można wyrachować połowę, a drugą skopiować. Dopiero spojrzenie całościowe pozwala ujrzeć symetrię w wynikach, gdyż część rzeczywista transformaty ma charakter parzysty oraz część urojona transformaty ma charakter nieparzysty. Wyrachowane wyniki prezentuje ilustracja 4. Istotna jest właściwa interpretacja tego co wychodzi. Jeżeli próbki (jako tablica w programie) są indeksowane od 0 do N-1 (patrz reprezentacja typów), to dla wyników programu dla części rzeczywistej transformaty jest (w przypadku liczby próbek=512):
  • pod indeksem 0 ukrywa się składowa stała,
  • pod indeksami 1 .. N/2-1 (dla N=512 wychodzi → 255 ) ukrywa się rzeczywiste widmo podstawowe,
  • pod indeksem N/2 (dla N=512 → 256) ukrywa się granica Nyquista (wartości zerowe: czytaj bliskie zeru),
  • pod indeksami N/2+1 .. N ukrywa się widmo powielone (można je traktować jako widmo dla częstotliwości ujemnych).
sig72.png

Z wydruków wygenerowanych przez program wynika, że

Kod: Zaznacz cały

Harm. 14 .re=55,5921916630 .im=-258,6466245991  Czest. =13,672 Hz
Harm. 15 .re=59,6666444390 .im=-1012,0509247227  Czest. =14,648 Hz
Harm. 16 .re=64,6764661718 .im=587,6381113817  Czest. =15,625 Hz
Harm. 17 .re=70,9411580662 .im=236,0272549218  Czest. =16,602 Hz
Harm. 18 .re=78,9490702103 .im=150,5277874873  Czest. =17,578 Hz
Harm. 19 .re=89,4852511184 .im=111,8331820397  Czest. =18,555 Hz
Harm. 20 .re=103,8962551572 .im=89,6928381388  Czest. =19,531 Hz
Harm. 21 .re=124,7003767125 .im=75,3094712585  Czest. =20,508 Hz
Harm. 22 .re=157,2149564675 .im=65,1846223126  Czest. =21,484 Hz
Harm. 23 .re=214,9390286604 .im=57,6509242013  Czest. =22,461 Hz
Harm. 24 .re=345,1209218922 .im=51,8125053015  Czest. =23,438 Hz
Harm. 25 .re=910,0019285673 .im=47,1448851844  Czest. =24,414 Hz
Harm. 26 .re=-1350,5951458007 .im=43,3204757202  Czest. =25,391 Hz
Harm. 27 .re=-382,0569497582 .im=40,1240654503  Czest. =26,367 Hz

energia prążka 15Hz oraz 25Hz rozlała się po okolicy i zasiliła sąsiednie prążki, gdyż próbkowanie nie trafiło w częstotliwość 15Hz (było 14.648Hz lub 15.625Hz) oraz w częstotliwość 25Hz (było 24.414Hz lub 25.391Hz).

Kod: Zaznacz cały

Harm. 254 .re=-5,5787259318 .im=0,0560186278  Czest. =248,047 Hz
Harm. 255 .re=-5,5786038331 .im=0,0280083080  Czest. =249,023 Hz
Harm. 256 .re=-5,5785631385 .im=0,0000001539  Czest. =250,000 Hz
Harm. 257 .re=-5,5786038352 .im=-0,0280080002  Czest. =250,977 Hz
Harm. 258 .re=-5,5787259360 .im=-0,0560183200  Czest. =251,953 Hz

Wyraźnie widać częstotliwość Nuquista (harmoniczna 256 stanowi oś symetrii: część rzeczywista harmonicznej 255 jest równa harmonicznej 257 oraz zachodzi podobne zjawisko między harmoniczną 254 i 258, oczywiście należy wziąć „lupę” numeryczną i ograniczyć się przykładowo do 7 cyfr po przecinku).
Po drugiej stronie częstotliwości Nyquista (250Hz) wystąpiły następujące prążki dla harmonicznej 486/487 i 496/497.

Kod: Zaznacz cały

Harm. 485 .re=-382,0569525196 .im=-40,1240822263  Czest. =473,633 Hz
Harm. 486 .re=-1350,5951174269 .im=-43,3205377508  Czest. =474,609 Hz
Harm. 487 .re=910,0019632004 .im=-47,1448413477  Czest. =475,586 Hz
Harm. 488 .re=345,1209292216 .im=-51,8124877931  Czest. =476,563 Hz
Harm. 489 .re=214,9390322480 .im=-57,6509126488  Czest. =477,539 Hz
Harm. 490 .re=157,2149588271 .im=-65,1846132795  Czest. =478,516 Hz
Harm. 491 .re=124,7003786520 .im=-75,3094634924  Czest. =479,492 Hz
Harm. 492 .re=103,8962571497 .im=-89,6928309664  Czest. =480,469 Hz
Harm. 493 .re=89,4852536543 .im=-111,8331749191  Czest. =481,445 Hz
Harm. 494 .re=78,9490741259 .im=-150,5277796696  Czest. =482,422 Hz
Harm. 495 .re=70,9411655781 .im=-236,0272444121  Czest. =483,398 Hz
Harm. 496 .re=64,6764897364 .im=-587,6380824086  Czest. =484,375 Hz
Harm. 497 .re=59,6665927333 .im=1012,0509472941  Czest. =485,352 Hz
Harm. 498 .re=55,5921747719 .im=258,6466245501  Czest. =486,328 Hz

Częstotliwość jednego to 474,609Hz-500Hz=-25.391Hz (co koreluje z harmoniczną 26). Podobnie jest w drugim skupieniu: 484,375 Hz-500Hz=-15.625Hz.
Przebieg „krzywej” (bo de facto to również ciąg próbek) wykazuje oś symetrii dla częstotliwości Nyquista. Nawet jeżeli na pierwszy rzut oka tego nie widać, to jednak ona zachodzi. Składowa stała czyli harmoniczna numer 0 powieli się jako harmoniczna numer 512 → następna od 511 (przy założeniu, że tablica jest indeksowana jako 0..N-1 → 0..511). Ideologicznie prawą połowę za częstotliwością Nyquista można traktować jako widmo częstotliwości ujemnych. W piśmie obrazkowym to wygląda następująco (ilustracja 5):
sig76.png

Teraz już wyraźnie widać, że widmo części rzeczywistej wykazuje symetrię parzystą.
W przypadku wyników części urojonej, wynik pokazuje ilustracja 6. Również tu istnieje określona interpretacja wyników:
  • pod indeksem 0 ukrywa się składowa stała, ma wartość zero gdyż wejściowa funkcja czasu jest czysto-rzeczywista,
  • pod indeksami 1 .. N/2-1 (dla N=512 wychodzi → 255 ) ukrywa się urojone widmo podstawowe,
  • pod indeksem N/2 (dla N=512 → 256) ukrywa się granica Nyquista (wartości zerowe: czytaj bliskie zeru),
  • pod indeksami N/2+1 .. N ukrywa się widmo powielone (można je traktować jako widmo dla częstotliwości ujemnych).
sig73.png

Również tu można dostrzec odpowiednie korelacje oraz dokonać „symetryzacji”, której wynik pokazuje ilustracja 7. Część urojona widma wykazuje symetrię nieparzystą.
sig77.png

Zastosowanie wyliczonych danych jako transformata odwrotna daje ponownie sygnał wejściowy, ilustracja 8 (z niewielkim błędem: różnicę między spróbkowanym sygnałem a wynikiem transformaty „tam” i „z powrotem” pokazuje ilustracja 9 – nie przekracza małego ułamka procentowego wartości sygnału oryginalnego).
sig74.png

sig75.png

Dla tropicieli:
DFT.7z
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse


Wróć do „Podstawy elektroniki - teoria i praktyka”

Kto jest online

Użytkownicy przeglądający to forum: Obecnie na forum nie ma żadnego zarejestrowanego użytkownika i 3 gości