Cyfrowe przetwarzanie sygnałów
Cyfrowe przetwarzanie sygnałów
w Auxerre, zm. 16 maja 1830 r. w Paryżu) – francuski
matematyk, który jest sprawcą
całego zamieszania.
Wstęp
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ć.
Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse
Re: Cyfrowe przetwarzanie sygnałów
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.
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.
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:
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.
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ę:
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ą):
Niech rzeczywisty sygnał będzie aproksymowany następującym szeregiem trygonometrycznym:
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 2π. 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):
Identycznie jest z funkcją cos.
Wszystkie rozkminy rachunkowe sprowadzają się do przedziału [-T/2 , T/2].
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 kω odpowiada określonej harmonicznej (jest całkowitoliczbową krotnością harmonicznej podstawowej ω). Mamy do obliczenia:
Warto tu zauważyć, że powyższą sumę można obliczyć „po kawałku”. Na początek zostaje obliczona:
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 kω.
Kolejnym fragmentem jest:
Z punktu widzenia rachunkowego koniecznością staje się rozpatrzenie dwóch wariantów, gdzie k≠n oraz k=n. Pierwszy przypadek w postaci:
gdzie jest k≠n, co oznacza „wyfiltrowanie” innej niż k-tej harmonicznej z całego zestawu harmonicznych:
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:
Kolejny fragment całki
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:
drugi przypadek:
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:
lub
Jest to kwestia punktu widzenia: czy brany jest pod uwagę okres T funkcji okresowej, czy okres (wynoszący 2π) 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).
Gdzie poszczególne elementy to:
Ostatecznie, skoro
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:
A poszczególne jego współczynniki wiadomo jak obliczyć.
Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse
Re: Cyfrowe przetwarzanie sygnałów
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.
Przypomniały mi się szkolne lata. Ostatni raz liczyłem całki w ubiegłym wieku.
X is work. Y is play. Z is keep your mouth shut."A. Einstein
Re: Cyfrowe przetwarzanie sygnałów
Re: Cyfrowe przetwarzanie sygnałów
Zegar pisze:Ostatni raz liczyłem całki w ubiegłym wieku.
No ja to robią ciągle.
To nie jest jeszcze przesadna matma, to jeszcze w miarę proste rzeczy
Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse
Re: Cyfrowe przetwarzanie sygnałów
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.
Można „wyciąć” odpowiedni fragment i przemieścić go w inne miejsce, to uzyskuje się ponownie ten sam przebieg (ilustracja 2).
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
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.
Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse
Re: Cyfrowe przetwarzanie sygnałów
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łę:
można zapisać w postaci (gdzie b0=0, choć w sumie to nieistotne):
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:
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.
Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse
Re: Cyfrowe przetwarzanie sygnałów
Weźmy dobrze znaną z świecie elektroniki funkcję piły. Niech jej szczegóły przedstawia ilustracja 1:
Dla prostoty rachunków, niech jej okres wynosi 2π sekund, czyli pulsacja podstawowa ω=1 rad/s → f=0.3183 hz. Składowa stała a0 to:
pozostałe współczynniki an szeregu trygonometrycznego to:
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:
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):
Funkcję piły można aproksymować następującym szeregiem trygonometrycznym:
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ą.
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.
Z parametrów funkcji wynika, że T=628 ms, f=1.59Hz, ω=10 rad/s.
Trochę rachunków. Składowa stała:
Współczynniki an są obliczone z formuły (A=5 jest amplitudą):
Współczynniki bn, to:
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ść
i finalnie szereg ten można zapisać w postaci (przy założonych parametrach):
Funkcja fali prostokątnej aproksymowana kilkoma początkowymi harmonicznym ma przebieg pokazany na ilustracji 4.
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)
- sygnał fali trójkątnej (ilustracja 6)
- sygnał piły (ilustracja 7)
Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse
Re: Cyfrowe przetwarzanie sygnałów
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.
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:
W chwili t0,t1 i t2 zachodzą związki:
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:
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
Wróćmy do formuły pozwalającej na obliczenie współczynników rozwinięcia funkcji w szereg trygonometryczny:
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
w ogólnym przypadku (całkowania jakiejś funkcji), należy utworzyć pola prostokątów i je zsumować.
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).
W przypadku drugiej harmonicznej jest to iloczyn wartości próbki z dwuokresową funkcją (ilustracja 4).
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).
Ł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.
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).
Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse
Re: Cyfrowe przetwarzanie sygnałów
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).
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)
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).
Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse
Re: Cyfrowe przetwarzanie sygnałów
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.
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
lub fazy (jako kąta utworzonego przez składowe harmonicznych, ilustracja 3).
W przypadku badanej funkcji, widmo fazowe pokazuje ilustracja 4.
Finalnie wyjściowa funkcja składa się z następujących elementów (ilustracja 5).
Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse
Re: Cyfrowe przetwarzanie sygnałów
Pozdrawiam! J23
BOB TAYLOR, PARC
Re: Cyfrowe przetwarzanie sygnałów
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
Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse
Re: Cyfrowe przetwarzanie sygnałów
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).
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
wykres fazy owych harmonicznych pokazuje ilustracja 3.
Pamiętając o tożsamości trygonometrycznej
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)
Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse
Re: Cyfrowe przetwarzanie sygnałów
gaweł pisze:Dla tropicieli, implementacja w Lazarus:
Trudno coś wytropić...
X is work. Y is play. Z is keep your mouth shut."A. Einstein
Re: Cyfrowe przetwarzanie sygnałów
Poniekąd temat na czasie
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
Re: Cyfrowe przetwarzanie sygnałów
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ę).
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 .
Dla tropicieli aktualny soft:
Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse
Re: Cyfrowe przetwarzanie sygnałów
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.
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).
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.
Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse
Re: Cyfrowe przetwarzanie sygnałów
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).
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:
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ć.
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ąć.
Na efekty nie trzeba długo czekać (ilustracja 3), bo w świecie DSP wszystko jest możliwe.
Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse
Re: Cyfrowe przetwarzanie sygnałów
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).
Próbkowany sygnał w funkcji czasu (ilustracja 2).
Rozkład harmonicznych w funkcji częstotliwości harmonicznych (ilustracja 3).
Dla tropicieli:
Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse
Re: Cyfrowe przetwarzanie sygnałów
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
Re: Cyfrowe przetwarzanie sygnałów
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.
Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse
Re: Cyfrowe przetwarzanie sygnałów
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?
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).
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).
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 ejω. Ile jakiejś ω zawiera się w f(t) określa wyrażenie:
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.
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.
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:
to mnożenie czegoś przez j realizuje na płaszczyźnie zespolonej obrót o kąt 90°, pokazuje to ilustracja 5.
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 ω.
Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse
Re: Cyfrowe przetwarzanie sygnałów
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.
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.
Odmiennie sprawa wygląda jeżeli wartości współczynników lekko będą się różnić od wyżej wymienionego.
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.
Niewielkie „zejście z kursu” daje wyraźny sygnał, że coś jest nie tak. Pokazuje to ilustracja 5.
Podobnie dla ω=3.0 (ilustracja 6 i 7).
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)?
Kolejne rysunki pokazują, jak subtelne jest to zjawisko.
Program do rolowania:
Prawdziwe słowa nie są przyjemne. Przyjemne słowa nie są prawdziwe.
Lao Tse
Re: Cyfrowe przetwarzanie sygnałów
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.
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:
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:
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
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):
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).
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):
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).
Również tu można dostrzec odpowiednie korelacje oraz dokonać „symetryzacji”, której wynik pokazuje ilustracja 7. Część urojona widma wykazuje symetrię nieparzystą.
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).
Dla tropicieli:
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 1 gość