Cyfrowe przetwarzanie sygnałów

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

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » wtorek 01 lut 2022, 22:23

Widok z okna

Weźmy do analizy dokładnie ten sam przykład okienkowany w sposób symetryczny. Poprzednio przedział wynikający z okienkowania był niesymetryczny, gdyż dotyczył (0,T). Nowy przedział, czyli funkcja okienkująca niech będzie wyglądała jak na ilustracji 1, czyli w obszarze zainteresowania znajduje się przedział (-T/2 .. T/2).
sig78.png

Nawet zadbałem, by próbkowana funkcja była dokładnie ta sama jak poprzednio. Została przesunięte w dziedzinie czasu w taki sposób, by pokryła się dla dodatnich argumentów czasu w poprzednim przebiegu. Pokazuje to ilustracja 2.
sig79.png

Spowodowało to pewną ingerencję w sposób generacji ciągu próbek. Fragment kodu jest następujący (argument czasu został przesunięty o połowę długości czasu okienkowania):

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 - SampleNumDiv2 ) ;
    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 ; (* Okres calego sprobkowanego sygnalu *)
end (* TDFTForm.InputSampleTable *) ;

Oczywiście to ma spore implikacje w algorytmie. Transformata Fouriera obecnie jest obrabiana w innym przedziale czasowym a to ma swoje odbicie na egzystencję w innej przestrzeni. Zmienia się początek nawijania.

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 - Pi ;
      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 *) ;

Nawijanie „tam” musie mieć ekwiwalentne rozwiązanie w nawijaniu „z powrotem”, inaczej nie trafimy do domu.

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 + Pi ;
      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 *) ;

Rezultat działania programu powinien spowodować pewną refleksję. Dane wejściowe dotyczą tej samej funkcji natomiast wyniki w znacznym stopniu się różnią. Przede wszystkim ma to swoje odbicie w zawartości harmonicznych. Część rzeczywistą pokazuje poniższa tabelka i ilustracja 3.

Kod: Zaznacz cały

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
Harm. 0 .re=-5082,4909935847 .im=-0,0000020849  Czest. =0,000 Hz
Harm. 1 .re=37,5737730002 .im=-3,1477393408  Czest. =0,977 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
( . . . )
Harm. 483 .re=-154,4966348165 .im=-35,0690825601  Czest. =471,680 Hz
Harm. 484 .re=-220,7931480689 .im=-37,4083567146  Czest. =472,656 Hz
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

sig80.png

Widmo to wykazuje jedynie harmoniczną o częstotliwości 25Hz (podobnie jak poprzednio energia rozlała się i zasiliła sąsiednie prążki widma → próbkowanie nie trafiło precyzyjnie w częstotliwość sygnału). Patrząc na formułę sygnału, z której był generowany ciąg próbek:10+5 sin ( 2 π f1 t ) + 7 cos ( 2 π f2 t), gdzie f1=15Hz i F2=25Hz, nie powinno dziwić, że widmo funkcji parzystej (cos) przebiło się przez te wszystkie „zawirowania na wałku” natomiast widmo funkcji nieparzystej (sin) nie było w stanie dotrzeć.
Oczywiście wystąpiła harmoniczna symetryczna w stosunku do częstotliwości Nyquista.
Wybrane elementy widma części urojonej pokazuje poniższa tabelka i całe widmo pokazuje ilustracja 4.

Kod: Zaznacz cały

Harm. 13 .re=52,2355272581 .im=-143,3330134370  Czest. =12,695 Hz
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. 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
Harm. 499 .re=52,2355152801 .im=143,3330137584  Czest. =487,305 Hz

sig81.png

Wskazuje ono jedynie harmoniczną o częstotliwości 15Hz (podobnie jak poprzednio energia rozlała się i zasiliła sąsiednie prążki widma). Analogicznie jak w części rzeczywistej, tu na wierzch wydobyły się jedynie echa funkcji nieparzystych (o częstotliwości 15Hz). Oczywiście wystąpiła harmoniczna po drugiej stronie częstotliwości Nyquista.
Odwrotna transformata zrobiła to co do niej należało, odtworzony sygnał pokazuje ilustracja 5.
sig82.png

Błąd odtworzenia funkcji (wynik transformaty odwrotnej) jest na podobnym poziomie (ilustracja 6).
sig83.png

Patrząc na historię zaistniałych zdarzeń, można zadać sobie pytanie: co się stało, że wyszło jak wyszło?
Generalnie zagadnienie wykazuje dużą wrażliwość na zaistniałe warunki. Wybranie jakiejś funkcji okienkującej wpływa na wynik końcowy: wystąpienie określonej przyczyny generuje ściśle określone skutki. Nastąpił przeciek energii nawet pomiędzy częścią rzeczywistą a częścią urojoną widma. Tu również występuje swoista symetria: określone skutki są następstwem ściśle określonych przyczyn.

Program dla tropicieli:
DFT.7z
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.
Ostatnio zmieniony środa 02 lut 2022, 15:11 przez gaweł, łącznie zmieniany 1 raz.

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

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

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » środa 02 lut 2022, 02:26

Liczby zespolone w C

Różni ludzie mają swoje różne sympatie i preferencje. By sympatycy programowania w C nie czuli się pokrzywdzeni w stosunku do sympatyków języka PASCAL (a podejrzewam, że ich jest większość), wspieram również tych drugich i nie pozostawiam samych sobie. Tak się składa, że w języku C typ liczb zespolonych jest standardem (w przeciwieństwie do PASCAL'a). Oznacza to, że istnieje typ reprezentujący te liczby (chyba od wersji C99). Dokładniej, to istnieją dwa typy reprezentujące liczby zespolone: z pojedynczą precyzją oraz z podwójną precyzją.

Kod: Zaznacz cały

float complex z1 ;
double complex z2 ;

By posługiwać się liczbami zespolonymi należy do programu dołączyć odpowiednie #include (a właściwie to dwa). Są to

Kod: Zaznacz cały

#include <complex.h>
#include <tgmath.h>

Drugi z nich zawiera funkcje matematyczne.
Pierwszy z nich (jak wskazuje sama nazwa) zawiera definicje typu oraz podstawowe funkcje operujące na liczba zespolonych. W tym include jest nawet zdefiniowana stała będąca jednostką urojoną, tak, że można używać jej w instrukcjach podstawienia. Stała ma identyfikator I (dużą literą). Można w programie zapisać :

Kod: Zaznacz cały

double complex Z ;
Z = 1.0 + 2.0 * I ; //Z=1+2i
printf ( "Z= %.8f + i*%.8f\n" , creal ( Z ) , cimag ( Z ) ) ;
printf ( "Z= %.8f + i*%.8f\n" , __real__ Z , __imag__ Z ) ;

W instrukcjach operujących na wybranych częściach liczby zespolonej można posługiwać się funkcją creal(z) lub cimag(z), gdzie z jest typu zespolonego, z tym, że to pozwoli wziąć odpowiedni fragment liczby zespolonej (nie da się tam włożyć). Można również użyć „zaklęcia” __real__ przed zmienną zespoloną do użycia części rzeczywistej oraz __imag__ do operowania na części urojonej. Tego zaklęcia można użyć „po obu stronach” operatora podstawienia. Legalne jest:

Kod: Zaznacz cały

__real__ Z = __real__ Z + ...coś.... ;
__imag__ Z = __imag__ Z + ...coś.... ;

Również można używać klasycznych operatorów arytmetycznych w wyrażeniach zespolonych. Fragment programu:

Kod: Zaznacz cały

double complex z1,z2,z3;

z1 = 1.0 + I ;
z2 = 1.0 - I ;
z3 = z1 * z2 ;
printf ( "Z1= %.8f + i*%.8f\n" , creal ( z1 ) , cimag ( z1 ) ) ;
printf ( "Z2= %.8f + i*%.8f\n" , creal ( z2 ) , cimag ( z2 ) ) ;
printf ( "Z3= %.8f + i*%.8f\n" , creal ( z3 ) , cimag ( z3 ) ) ;

da następujący wynik (zgodny z arytmetyką zespoloną):

Kod: Zaznacz cały

Z1= 1.00000000 + i*1.00000000
Z2= 1.00000000 + i*-1.00000000
Z3= 2.00000000 + i*0.00000000

Poniżej są dwa „pakunki” programów implementujących dyskretną transformatę Fouriera w wersji okienka nieparzystego i parzystego (także, sympatycy języka PASCAL i języka C powinni być usatysfakcjonowani jak i również inni poligloci, znajomość wielu języków często jest przydatna). Programy są o identycznej funkcjonalności jak wyżej omówione programy w języku PASCAL.

Pakunek nieparzysty:
dft_np.7z

Pakunek parzysty:
dft_pa.7z
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

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

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

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » piątek 04 lut 2022, 18:34

Przekłamania i fałszowanie rzeczywistości


Wcześniej mowa była o szeregu trygonometrycznym. On również przetwarza funkcję czasu (oryginał) na szereg trygonometryczny, w którym mowa jest o harmonicznych. W pewnym sensie można mówić tu o widmie (aczkolwiek samo widmo jest pojęciem zespolonym). Jeżeli „przetworzyć” do postaci szeregu jakąś funkcję uzyskuje się ciąg czegoś, co można koncepcyjnie przyrównać do widma. Jest jedynie kwestią punktu widzenia, że w szeregu trygonometrycznym uzyskuje się amplitudy fal sin oraz cos natomiast w transformacie Fouriera wynikłe liczby nie są wprost amplitudami, jednak można stworzyć coś w rodzaju wspólnej płaszczyzny porozumienia i odnieść to do wspólnego mianownika. Rozpatrzmy program, który wygenerował ciąg próbek:

Kod: Zaznacz cały

void InputSampleTable ( real           SamplInt ,
                        real *         TotalTime ,
                        real           TimeTable [ ] ,
                        real           SampleTable [ ] ,
                        unsigned short Samples )
{
  unsigned short Loop ;   
  real Time ;
  real Angle ;
  /*----------------------------------------------------------*/   
  for ( Loop = 0 ; Loop <= Samples ; Loop ++ )
  {
    TimeTable [ Loop ] = SamplInt * ( real ) ( Loop - SampleNumDiv2 ) ;
    Angle = DoublePi * ( real ) ( Loop - SampleNumDiv2 ) / ( real) ( Samples ) ;
    SampleTable [ Loop ] = 10.0 + 5.0 * sin ( 15.0 * Angle ) + 7.0 * cos ( 25.0 * Angle ) ;
  } /* for */ ;
  * TotalTime = SampleNum * SamplInt ;
} /* InputSampleTable */

Wygenerowany sygnał pokazuje ilustracja 1.
sig90.png

Dla tego ciągu zostały obliczone współczynniki szeregu trygonometrycznego oraz w oparciu o uzyskane współczynniki została dokonana operacja odwrotna z jednoczesnym określeniem różnicy (między wygenerowanymi wartościami jako ciąg wejściowy próbek a po przejściu „tam” i „z powrotem”) oraz określeniem błędu wyrażonego w procentach w stosunku do danych oryginalnych. Sygnał wynikowy pokazuje ilustracja 2 (generalnie pokrywa się z sygnałem wejściowym z wyjątkiem samego początku oraz samego końca) oraz stopień błędów wyrażony w procentach ilustracja 3.
sig91.png

sig92.png

Tabelka z wynikami jest następująca:

Kod: Zaznacz cały

Testy
Probka [ 0 ] = 3.000000, obl. wart.=6.005859 [rozn.=3.005859 :100.195312%]
Probka [ 1 ] = 2.411658, obl. wart.=2.405799 [rozn.=-0.005859 :-0.242962%]
Probka [ 2 ] = 2.477431, obl. wart.=2.483291 [rozn.=0.005859 :0.236511%]
Probka [ 3 ] = 3.138474, obl. wart.=3.132615 [rozn.=-0.005859 :-0.186695%]
( . . . )
Probka [ 509 ] = 8.384371, obl. wart.=8.378512 [rozn.=-0.005859 :-0.069885%]
Probka [ 510 ] = 6.076382, obl. wart.=6.082241 [rozn.=0.005859 :0.096429%]
Probka [ 511 ] = 4.242057, obl. wart.=4.236198 [rozn.=-0.005859 :-0.138127%]
Probka [ 512 ] = 3.000000, obl. wart.=6.005859 [rozn.=3.005859 :100.195311%]

Dające się zauważyć piki w środku wynikają z ziarnistości liczb rzeczywistych w komputerze, co pokazuje tabelka:

Kod: Zaznacz cały

Probka [ 79 ] = 1.030414, obl. wart.=1.024555 [rozn.=-0.005859 :-0.568643%]
Probka [ 80 ] = 0.022365, obl. wart.=0.028224 [rozn.=0.005859 :26.199272%]
Probka [ 81 ] = -0.301668, obl. wart.=-0.307527 [rozn.=-0.005859 :1.942327%]
Probka [ 82 ] = 0.123067, obl. wart.=0.128927 [rozn.=0.005859 :4.761114%]
Probka [ 83 ] = 1.298601, obl. wart.=1.292742 [rozn.=-0.005859 :-0.451207%]
( . . . )
Probka [ 144 ] = -1.769423, obl. wart.=-1.763564 [rozn.=0.005859 :-0.331146%]
Probka [ 145 ] = -1.132114, obl. wart.=-1.137974 [rozn.=-0.005859 :0.517560%]
Probka [ 146 ] = 0.246828, obl. wart.=0.252687 [rozn.=0.005859 :2.373870%]
Probka [ 147 ] = 2.242991, obl. wart.=2.237132 [rozn.=-0.005859 :-0.261230%]
Probka [ 148 ] = 4.684244, obl. wart.=4.690103 [rozn.=0.005859 :0.125087%]
( . . . )
Probka [ 242 ] = 4.488323, obl. wart.=4.494182 [rozn.=0.005859 :0.130547%]
Probka [ 243 ] = 1.960085, obl. wart.=1.954225 [rozn.=-0.005859 :-0.298935%]
Probka [ 244 ] = -0.020138, obl. wart.=-0.014279 [rozn.=0.005859 :-29.096216%]
Probka [ 245 ] = -1.303952, obl. wart.=-1.309811 [rozn.=-0.005859 :0.449355%]
Probka [ 246 ] = -1.799914, obl. wart.=-1.794054 [rozn.=0.005859 :-0.325536%]

I tu zachodzi dziwna sytuacja: na samych krańcach wynik przetworzony „tam” i „z powrotem” różni się o 100% w stosunku do sygnału wejściowego. Nie zależy to od gęstości próbkowania ani od precyzji rachunków → wystąpi zawsze. Wcześniej, gdy mowa była o rozwijaniu funkcji w szereg trygonometryczny takie zjawisko nie zaistniało. Jest to kwestia doboru pewnych cech występujących w sygnale. Nie będę ukrywał, że wyszło tak, bo taki był mój zamiar i robiłem to tendencyjnie.
Warto zadać sobie pytanie: a dlaczego wyszło jak wyszło?
W tym przypadku dotknęło nas zjawisko nakładania się widm zwane aliasingiem. Spróbkowany sygnał o n+1 próbkach generuje harmoniczną zerową oraz n harmonicznych fal sin oraz cos. Całe zagadnienie jest okresowe: po każdym okresie następuje powtórka wszystkich zjawisk. Dochodzą do tego nawarstwiające się problemy z „przeszłości” i „przyszłości”. Każda harmoniczna o numerze 256 jest obciążona zaszłością z widma sąsiadującego. Próbę wyjaśnienia problemów pokazuje ilustracja 4.
sig95.png

Harmoniczna 256 jest wspólna dla sąsiadujących ze sobą zestawów, toteż na wynik odtworzenia sygnału mają wpływ niejako obce elementy, które dodają się w 100%.
Dlaczego ten efekt nie wystąpił wcześniej? Odpowiedź jest banalnie prosta: wystąpił, tylko nie był widoczny, gdyż dodanie do zera 100% wartości nadal jest zerem. W jakimś sensie jest to sposób na pozbycie się takich problemów, ale nie zawsze jest to możliwe.
Można nawet sobie pokombinować, by efekt ten zminimalizować. Przykładowo wygenerowałem sobie dane, które po zamianie na szereg trygonometryczny i podróży z powrotem dały następujące wyniki:

Kod: Zaznacz cały

Testy
Probka [ 0 ] = 2.623226, obl. wart.=6.169453 [rozn.=3.546227 :135.185704%]
Probka [ 1 ] = 2.364933, obl. wart.=2.359235 [rozn.=-0.005698 :-0.240920%]
Probka [ 2 ] = 2.736134, obl. wart.=2.741832 [rozn.=0.005698 :0.208235%]
Probka [ 3 ] = 3.651474, obl. wart.=3.645776 [rozn.=-0.005698 :-0.156035%]
Probka [ 4 ] = 4.979406, obl. wart.=4.985104 [rozn.=0.005698 :0.114423%]
( . . . )
Probka [ 509 ] = 9.647593, obl. wart.=9.641895 [rozn.=-0.005698 :-0.059057%]
Probka [ 510 ] = 7.169316, obl. wart.=7.175014 [rozn.=0.005698 :0.079472%]
Probka [ 511 ] = 5.085967, obl. wart.=5.080270 [rozn.=-0.005698 :-0.112026%]
Probka [ 512 ] = 3.540529, obl. wart.=6.169453 [rozn.=2.628924 :74.252280%]

Jak widać, nic w przyrodzie nie ginie i „oszczędności” w jednym miejscu skutkują „zwiększonymi wydatkami” w drugim miejscu (a per saldo, średnia arytmetyczna , wychodzi praktycznie na to samo).
Zjawiska tego nie obserwuje się przy transformacie Fouriera. Dokładnie ta sama funkcja spróbkowana w identyczny sposób daje po przejściu „tam” (transformata DFT) i „z powrotem” dokładnie to samo (ilustracja 5).
sig93.png

Te „dokładnie” można traktować za rzeczywiście dokładnie, gdyż różnica wyrażona w procentach jest praktycznie zerowa, co pokazuje ilustracja 6.
sig94.png

A dlaczego tak się dzieje?
Już pisałem o tym (tak między wierszami). Po przejściu danych przez transformatę DFT mamy (dla tablicy indeksowanej od 0 do N-1, konkretnie dla n=512):
  • pod indeksem 0 ukrywa się składowa stała,
  • pod indeksami 1 .. N/2-1 (czyli od 1 do 255) ukrywa się widmo podstawowe,
  • pod indeksem N/2 (czyli dla 256) ukrywa się granica Nyquista,
  • pod indeksami N/2+1 .. N (czyli od 257 do 511) ukrywa się widmo powielone.
Mamy od granicy Nyquista w jedną stronę 256 prążków, w drugą stronę 255 prążków, które stanowią cały okres próbkowania. Liczby w obie strony nie są równe, więc nie spowodują nałożenia się na siebie widm powielonych. Podobna koncepcja występuje w kodach liczb całkowitych ze znakiem (czyli U2), liczb dodatnich zapisanych na 8 bitach jest 127, liczb ujemnych jest 128 i jedno zero, co daje łącznie 256 kombinacji.
Użyte softy:
fseries.7z
dft.7z
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

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

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

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » sobota 05 lut 2022, 22:59

Dążenie do szczegółów

Może się tak zdarzyć, że widmo sygnału jest dosyć rozmyte (zbyt wiele energii prążka wyciekło do okolicy), natomiast nas interesują szczegóły. A jak wiadomo, to właśnie szczegóły często są tak decydujące, że mogą determinować nasze działania.
Został spróbkowany jakiś sygnał (analizując tekst programu, to nawet wiadomo jaki: s(t)=10+5 sin(2 π f1 t)+7 cos(2 π f2 t), gdzie f1=16.5Hz oraz f2=28.3Hz). Częstotliwość próbkowania wynosi 500Hz, więc nawet nieuzbrojonym okiem widać, ze próbki widma nie koniecznie musiały trafić w częstotliwości rzeczywistego sygnału. Ma to swoje odbicie w widmie (poniższa tabelka i ilustracja 1).

Kod: Zaznacz cały

Okres probkowania: 2.000 ms
Okres calego sygnalu: 512.000 ms
Czestotliwosc probkowania: 500.000 Hz
Czestotliwosc Nyquista: 250.000 Hz
Czestotliwosc podstawowa (pierwsza harmoniczna): 1.953 Hz
Pulsacja podstawowa (pierwsza harmoniczna): 12.27184630 rad/s

Harm. 0 .re = -2593.9969522045  .im= 0.0000000000  Ampl= 10.1328005945  Czest. = 0.000 Hz
Harm. 1 .re = -34.1852440119  .im= 5.7331620934  Ampl= 0.2708020377  Czest. = 1.953 Hz
Harm. 2 .re = -34.7610914871  .im= 11.9751317894  Ampl= 0.2872342163  Czest. = 3.906 Hz
Harm. 3 .re = -35.7590484045  .im= 19.3978774885  Ampl= 0.3178244830  Czest. = 5.859 Hz
Harm. 4 .re = -37.2426374849  .im= 29.1234550149  Ampl= 0.3693577552  Czest. = 7.813 Hz
Harm. 5 .re = -39.3151353984  .im= 43.4506459440  Ampl= 0.4577910687  Czest. = 9.766 Hz
Harm. 6 .re = -42.1394460165  .im= 68.3193116290  Ampl= 0.6271088082  Czest. = 11.719 Hz
Harm. 7 .re = -45.9750147220  .im= 125.9461629096  Ampl= 1.0474618787  Czest. = 13.672 Hz
Harm. 8 .re = -51.2497627024  .im= 436.6153319353  Ampl= 3.4344756427  Czest. = 15.625 Hz
Harm. 9 .re = -58.7113818161  .im= -375.4794445292  Ampl= 2.9690772807  Czest. = 17.578 Hz
Harm. 10 .re = -69.7814725513  .im= -140.2075503521  Ampl= 1.2235385467  Czest. = 19.531 Hz
Harm. 11 .re = -87.5191783826  .im= -88.8769332194  Ampl= 0.9744889186  Czest. = 21.484 Hz
Harm. 12 .re = -119.9078368523  .im= -66.1765909526  Ampl= 1.0699768671  Czest. = 23.438 Hz
Harm. 13 .re = -196.3807840874  .im= -53.2661529027  Ampl= 1.5896603358  Czest. = 25.391 Hz
Harm. 14 .re = -586.8710577856  .im= -44.8748453945  Ampl= 4.5983142598  Czest. = 27.344 Hz
Harm. 15 .re = 554.1758916672  .im= -38.9469321288  Ampl= 4.3401779829  Czest. = 29.297 Hz
Harm. 16 .re = 184.7301531357  .im= -34.5136099492  Ampl= 1.4681768070  Czest. = 31.250 Hz
Harm. 17 .re = 109.8493656352  .im= -31.0578509001  Ampl= 0.8918396745  Czest. = 33.203 Hz
Harm. 18 .re = 77.7819753174  .im= -28.2781530718  Ampl= 0.6465847791  Czest. = 35.156 Hz

Amplitudą harmonicznej jest przeskalowana (podzielona przez długość połowy okresu całego spróbkowanego sygnału z wyjątkiem składowej stałej, gdzie należy podzielić przez całą długości okresu próbkowania) wartość bezwzględna zespolonego elementu widmowego.
Kawałek programu do rachunków może być następujący:

Kod: Zaznacz cały

  for ( Loop = 0 ; Loop < SampleNum ; Loop ++ )
    MagnTable [ Loop ] = sqrt ( __real__ DFTRes [ Loop ] * __real__ DFTRes [ Loop ]  +
                                __imag__ DFTRes [ Loop ] * __imag__ DFTRes [ Loop ] ) /
                               ( double ) SampleNumDiv2 ;
  MagnTable [ 0 ] /= 2.0 ;

sig96.png

By dokładniej określić wartość częstotliwości można „górkę” amplitudy aproksymować kawałkiem funkcji kwadratowej. Taka funkcja jest wystarczająco dobrze poznana z punktu widzenia analitycznego. Mając znane trzy różne punkty należące do paraboli można wyznaczyć wszystkie jej parametry. Jeżeli te punkty będą: P1 (x1,y1), P2 (x2,y2), P3 (x3,y3) to można określić równanie paraboli przechodzące przez zadane punkty. Ogólne równanie ma postać: f(x)=a x2 +b x +c, gdzie a≠0. Jeżeli zadane punkty należą do paraboli, to ich współrzędne muszą spełniać to równanie.
part18_form01.png

W powyższym układzie równań, poszukiwanymi niewiadomymi są współczynniki równania kwadratowego a, b, c. Równanie można rozwiązać metodą wyznaczników. Wyznacznik układu jest następujący:
part18_form02.png

Kolejne wyznaczniki:
part18_form03.png

part18_form04.png

part18_form05.png

współczynniki oblicza się dzieląc odpowiednie wyznaczniki:
part18_form06.png

Mając znane współczynniki, bez problemów można określić położenie wierzchołka paraboli (ilustracja 2).
sig97.png

Dokładnie tej technologii używa program do uszczegółowienia zarówno częstotliwości jak i wartości amplitudy, przykładową procedurę pokazuje poniższy kawałek programu:

Kod: Zaznacz cały

static void ComputeFrequency ( double H1 ,
                               double M1 ,
                               double H2 ,
                               double M2 ,
                               double H3 ,
                               double M3 )
{
  double DetW ;
  double DetWA ;
  double DetWB ;
  double DetWC ;
  double delta ;
  double A ;
  double B ;
  double C ;
  double f ;
  double Mag ;
  /*----------------------------------------------------------------------*/
  fprintf ( OutFile , "P1 = %.8f  %.8f \n" , H1 , M1 ) ;
  fprintf ( OutFile , "P2 = %.8f  %.8f \n" , H2 , M2 ) ;
  fprintf ( OutFile , "P3 = %.8f  %.8f \n" , H3 , M3 ) ;
  DetW = H1*H1*H2+H3*H3*H1+H2*H2*H3-H3*H3*H2-H1*H1*H3-H2*H2*H1 ;
  DetWA = H2*M1+H1*M3+H3*M2-H2*M3-H3*M1-H1*M2 ;
  DetWB = H1*H1*M2+H3*H3*M1+H2*H2*M3-H3*H3*M2-H1*H1*M3-H2*H2*M1 ;
  DetWC = H1*H1*H2*M3+H3*H3*H1*M2+H2*H2*H3*M1-
          H3*H3*H2*M1-H1*H1*H3*M2-H2*H2*H1*M3 ;
  A = DetWA / DetW ;
  B = DetWB / DetW ;
  C = DetWC / DetW ;
  delta = B * B - 4.0 * A * C ;
  f = - B / ( 2.0 * A ) ;
  Mag = - delta / ( 4.0 * A ) ;
  fprintf ( OutFile , "Czest. = %.8f\n" , f ) ;
  fprintf ( OutFile , "Ampl. = %.8f\n" , Mag ) ;
} /* ComputeFrequency */

Pozostaje jeszcze kwestia lokalizacji potencjalnych miejsc na widmie amplitudowym sygnału. Realizuje to odpowiednie procedura poszukująca swoistych ekstremów na charakterystyce. Przykładowy kod jest następujący:

Kod: Zaznacz cały

void LocateFrequency ( double         FreqTab [ ] ,
                       double         MagnitTab [ ] ,
                       unsigned short TabSize )
{
  unsigned short Loop ;
  double Average ;
  SampleArrayType DerivativeTable ;
  /*----------------------------------------------------------------------*/
  Average = 0.0 ;
  for ( Loop = 1 ; Loop < TabSize ; Loop ++ )
  {
    Average += MagnitTab [ Loop ] ;
    DerivativeTable [ Loop ] = 0.0 ;
  } /* for */ ;
  Average /= ( double ) TabSize ;
  for ( Loop = 2 ; Loop < TabSize ; Loop ++ )
  {
    if ( MagnitTab [ Loop ] >= Average )
      DerivativeTable [ Loop ] = MagnitTab [ Loop ] - MagnitTab [ Loop - 1 ] ;
  } /* for */ ;
  fprintf ( OutFile , "Wart. sred. = %.8f\n" , Average ) ;
  for ( Loop = 1 ; Loop < TabSize ; Loop ++ )
  {
    if ( ( DerivativeTable[Loop]>0.0 ) && ( DerivativeTable[Loop+1 ]<=0.0 ) )
    {
      fprintf ( OutFile , "\n\nExtremum dla n=%d\n" , Loop ) ;
      ComputeFrequency ( FreqTab [ Loop - 1 ] , MagnitTab [ Loop - 1 ] ,
                         FreqTab [ Loop ] , MagnitTab [ Loop ] ,
                         FreqTab [ Loop + 1 ] , MagnitTab [ Loop + 1 ] ) ;
    } /* if */ ;
  } /* for */ ;
} /* LocateFrequency */

Przede wszystkim poszukiwanie odpowiednich lokacji zaczyna się z pominięciem elementu indeksowanego jako zero (gdyż tam ukrywa się składowa stała o częstotliwości zero → jest poza obszarem zainteresowania) i kończy się na częstotliwości Nyquista. Algorytm w pierwszej kolejności oblicza wartość średnią sygnału by odsiać i pominąć w przetwarzaniu szum. Wszystkie wartości amplitudy poniżej wartości średniej są odcięte i nie są brane pod uwagę. Rozwiązanie problemu pokazuje ilustracja 3 i sprowadza się do określenia przyrostów w danych przedziałach (punktach). W oparciu o tę bardzo przybliżoną pochodną sygnału poszukiwane są miejsca, gdzie pochodna zmienia znak z dodatniej na ujemną (ewentualnie na zero, co sprowadza się do zejścia z ekstremum na sygnał o wartości poniżej średniej).
sig98.png

Ten algorytm zapuszczony na przykładowe dane znalazł wszystkie elementy co miał znaleźć, wynik pokazuje tabelka:

Kod: Zaznacz cały

Wart. sred. = 0.28843682


Extremum dla n=8
P1 = 13.67187500  1.04746188
P2 = 15.62500000  3.43447564
P3 = 17.57812500  2.96907728
Czest. = 16.28289145
Ampl. = 3.59629508


Extremum dla n=14
P1 = 25.39062500  1.58966034
P2 = 27.34375000  4.59831426
P3 = 29.29687500  4.34017798
Czest. = 28.16597984
Ampl. = 4.88779364

Wyniki może nie koniecznie powalają swoją precyzją, ale takie są zawsze lepsze od żadnych. Zawsze zachodzi poprawienie jakości uzyskanych danych.

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

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

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

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » poniedziałek 14 lut 2022, 00:34

Szybka transformata Fouriera (FFT)

Nadszedł czas na bardziej zaawansowane działania.

Szybka transformata Fouriera (FFT) koncepcyjnie nie różni się od dyskretnej transformaty Fouriera (DFT), a właściwie stanowi jej szczególny przypadek. Ta szczególność dotyczy liczby próbek, która w istotny sposób rzutuje na optymalność rachunkową.
Rozpatrzmy ciąg próbek sygnału w ilości 2 sztuk. Realizując algorytm DFT dla tych próbek (s[n]) kolejne harmoniczne (H[n]) wyliczane są w oparciu o następującą formułę:
part19_form01.png

Jeżeli przez R(α)=e-j α oznaczymy operator obrotu w przestrzeni zespolonej o kąt α, to:
part19_form02.png

Ponieważ e-j 0=cos(0)-j sin(0)=1 oraz e-j π=cos(π)-j sin(π)=-1 to powyższe można zapisać jako
part19_form03.png

co graficznie pokazuje ilustracja 1.
sig100.png



Rozpatrzmy ciąg próbek sygnału w ilości 4 sztuk. Realizując algorytm DFT dla tych próbek kolejne harmoniczne wyliczane są w oparciu o następującą formułę:
part19_form04.png

Biorąc pod uwagę redukcję kątów (obrót o kąt 2 k π + α=α) powyższą formułę można przedstawić jako:
part19_form05.png

Podobnie jak wyżej można już dokonać pewnych uproszczeń:
part19_form06.png

Powyższe rachunki w piśmie obrazkowym pokazuje ilustracja 2.
sig101.png

Przyglądając się uważnie tej formule, daje się dostrzec, że idea rozwiązania jest fraktalna, niejako rekurencyjna. Jednak ta prostota jest okupiona koniecznością ułożenia próbek sygnału w innej kolejności (o czym jest dalej).
Podobnie jak poprzednio, mamy dwie dwupunktowe próbki (s[0] i s[2] oraz s[1] i s[3]), które są przetworzone jako próbki dwupunktowe. Wyniki tego przetworzenia są ponownie przetworzone w koncepcyjnie identyczny sposób.


Rozpatrzmy ciąg próbek sygnału w ilości 8 sztuk. Realizując algorytm DFT dla tych próbek kolejne harmoniczne wyliczane są w oparciu o następującą formułę:
part19_form07.png

Co daje po redukcjach wielokrotności kątów następujące zależności:
part19_form08.png

Działania pokazuje ilustracja 3.
sig102.png

Niektóre operacje mogą budzić wątpliwości, są to jednak tylko pozory. Często występująca tożsamość to:
part19_form09.png

czasami występująca:
part19_form10.png

part19_form11.png

Wyjaśnienia wymaga jeszcze sposób przeindeksowania tablicy próbek. Zachodzi tu inwersja kolejności bitów w bitowym rozwinięciu numeru indeksu. Przykładowo w zakresie 3 bitów (próbki numerowane od 0 do 7) inwersja bitów pokazuje ilustracja 4.
sig103.png

Biorą pod uwagę uwagę dwie próbki, rachunki są mało skomplikowane → występuje obrót o 0 lub π radianów (czyli dodać lub odjąć). Fraktalne złożenie dwóch zestawów dwupunktowych próbek wymaga podwojonej ich liczby. Kolejne fraktalne złożenie dwóch zestawów próbek 4-ro punktowych wymaga kolejnego podwojenia liczby próbek. Jak łatwo się domyślić, liczba próbek musi być „liczbą bożą” (cytując dialog z filmu „Sami swoi”, dwa to liczba boża – i każda jej potęga). Przy spełnieniu takich wymagań, zachodzi swoista symetria rachunków i zamiast wielokrotnie wyliczać wartości funkcji trygonometrycznych można znacząco zmniejszyć ich liczbę. Warto pamiętać, że komputery nie mają tablic matematycznych i funkcje trygonometryczne rozwijają w szereg Maclaurina, a to generuje zapotrzebowanie na dużą moc obliczeniową, gdyż (liczba branych pod uwagę wyrazów rzutuje na precyzję rachunków):
part19_form12.png

Tak z grubsza wygląda algorytm zaproponowany przez Cooleya-Tukeya. Co ciekawe, nie byli oni pierwszymi, którzy wpadli na ten pomysł (swoją drogą trzeba mieć wyobraźnię, by dostrzec pewne detale). Algorytm ten, włączając jego rekurencyjne działanie, został wynaleziony około roku 1805 przez Carla Friedricha Gaussa, ale jego praca nie była powszechnie uznawana (opublikowana dopiero pośmiertnie). Różne ograniczone formy zostały ponownie odkryte kilka razy, zarówno w XIX w., jak i na początku XX wieku. Transformacje FFT stały się popularne po tym jak Cooley (IBM) i Tukey (Princeton) opublikowali w 1965 roku artykuł o ponownym wynalezieniu algorytmu i opisaniu metody na jego użycie na komputerze. Nie bez znaczenia jest również fakt, że dopiero pojawienie się odpowiedniego wyposażenia technicznego (komputery) spowodowało wręcz eksplozję zaiteresowania tym algorytmem. Mimo że Gauss opisał ten sam algorytm, został on zrealizowany dopiero kilka lat po artykule Cooleya i Tukeya. Powstało wiele realizacji i rozwiązań programowych i temat można uważać nadal za otwarty, gdyż z pewnością nie zostało powiedziane jeszcze ostatnie słowo w temacie.
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

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

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

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » środa 16 lut 2022, 15:14

Motylki w brzuchu, czyli
implementacja szybkiej transformaty Fouriera (FFT)


Algorytm liczenia szybkiej transformaty Fouriera jest w gruncie rzeczy algorytmem rekurencyjnym. Jak wiadomo, każdy algorytm rekurencyjny można zrealizować w sposób iteracyjny. Jednak nie spotkałem się z rekurencyjną jego realizacją i zapewne wynika to z jakichś powodów, ale nie wnikałem w tą problematykę.
Cechą charakterystyczną algorytmu liczenia FFT są „motylki”. Jak to często bywa, chyba ktoś dla żartu tak nazwał określony fragment algorytmu i tak już zostało. Jako motylek należy rozumieć fragment kodu:

Kod: Zaznacz cały

AddPart = OutTable [ AddPartIndex ] ;
SubPart = OutTable [ SubPartIndex ] * Rotation ;
OutTable [ AddPartIndex ] = AddPart + SubPart ;
OutTable [ SubPartIndex ] = AddPart - SubPart;

W piśmie obrazkowym odpowiada to zawsze parze próbek w określonej strefie obrotów jak jest pokazane na ilustracji 1.
sig105.png

W przypadku 8 próbek możliwymi obrotami są:
  • o 0 oraz π radianów w pierwszej strefie obrotów
  • o 0,π/2, π i 3π/2 (będące w rzeczywistości obrotem -π/2) radianów w drugiej strefie obrotów,
  • o 0,π/4, π/2, 3π/4, π, 5π/4, 3π/2, 7π/4 w kolejnej strefie obrotów (kąty obrotów powyżej π radianów są traktowane jako obrót o kąt dopełnienia do radianów w stronę przeciwną).
Pokazuje to ilustracja 2.
sig106.png

Implementacja algorytmu FFT przejścia z funkcji czasu do przestrzeni zespolonej (czyli podróży tam) ma ciekawą dodatkową własność, można użyć jej do realizacji przejścia z powrotem. Sam algorytm różni się w jednym szczególe (złośliwi mogą twierdzić, że w dwóch). Przede wszystkim, co jest elementem istotnym, różni się kierunkiem nawijania na koło zespolone. Ten szczegół jest ujęty w parametrach wywołania funkcji realizującej obliczenia.
Drugim szczegółem jest to, że algorytm FFT „nie skaluje” wyliczonych próbek. W przypadku podróży tam, dane odpowiadają widmu (nie są to amplitudy harmonicznych). By uzyskać amplitudy, należy uzyskane dane przeliczyć. Analogicznie przy podróży z powrotem, należy uzyskane dane podzielić przez długość okresu próbkowania (czyli liczbę próbek). Jeżeli odbędziemy podróż tam w towarzystwie określonych próbek, to podróż z powrotem (z podzieleniem wyniku przez liczbę próbek) da dokładnie to samo (o czym można się przekonać samemu).
Skoro podróż tam bazuje (jako dane wejściowe) na tablicy próbek funkcji czasu i daje próbki w przestrzeni zespolonej, to podróż z powrotem bazuje (jako dane wejściowe) na próbkach widma i daje próbki w funkcji czasu. Ta dziwna symetria oznacza, że próbki w funkcji czasu należy umieścić w tablicy zespolonej w częściach rzeczywistych zerując części urojone. Operacja „z powrotem” również da wyniki jako dane zespolone, lecz część urojona będzie zawierać zera.
To oznacza, że dane wejściowe jak i wyjściowe są reprezentowane przez ten sam typ danych. Ma to swoje odbicie w oprogramowaniu:

Kod: Zaznacz cały

#define Pi                  ( ( double ) 3.141592653589793 )
#define DoublePi            ( ( double ) 6.283185307079586 )


#define MaxSampleTableSize  14
#define SampleNumDiv2       ( 256 )
#define SampleNum           ( 2 * SampleNumDiv2 )

typedef double complex SampleArrayType [ SampleNum ] ;
typedef double RealSampleArrayType [ SampleNum ] ;

static SampleArrayType InpSampleTable ;
static SampleArrayType FFTTable ;
static SampleArrayType IFFTTable ;

Podróż tam to:

Kod: Zaznacz cały

void FFT ( double complex SpectrumTable [ ] ,
           double complex SampleTable [ ] ,
           unsigned short SampleTabSize )
{
  /*----------------------------------------------------------------*/
  FourierFransformation ( 1.0, SpectrumTable, SampleTable, SampleTabSize) ;
} /* FFT */

podróż z powrotem to:

Kod: Zaznacz cały

void IFFT ( double complex SampleTable [ ] ,
            double complex SpectrumTable [ ] ,
            unsigned short SampleTabSize )
{
  unsigned short Loop ;
  /*----------------------------------------------------------------------*/
  FourierFransformation( -1.0, SampleTable, SpectrumTable, SampleTabSize );
  for ( Loop = 0 ; Loop < SampleTabSize ; Loop ++ )
    SampleTable [ Loop ] = SampleTable [ Loop ] / ( double ) SampleTabSize ;
} /* IFFT */

zawiera wspólną funkcję do obróbki danych. Jednym z parametrów jest kierunek nawijania (1.0 lub -1.0).
Funkcja do przetwarzania danych to:

Kod: Zaznacz cały

void FourierFransformation ( double         Signum ,
                             double complex OutTable [ ] ,
                             double complex InpTable [ ] ,
                             unsigned short TableLength )
{
  double complex Tmp ;
  double complex AddPart ;
  double complex SubPart ;
  double complex TwiddleFactor ;
  double complex Rotation ;
  unsigned short NumberOfBits ;
  unsigned short StageIndex ;
  unsigned short Loop ;
  unsigned short RIndex ;
  unsigned short IntIndex ;
  unsigned short Index ;
  unsigned short AddPartIndex ;
  unsigned short SubPartIndex ;
  double Period ;
  double Angle ;
  /*----------------------------------------------------------------------*/
  NumberOfBits = NumberOfBitsNeeded ( TableLength ) ;
  if ( NumberOfBits == 0 )
    return ;
  Period = Signum * DoublePi ;
  for ( Loop = 0 ; Loop < TableLength ; Loop++ )
  {
    RIndex = ReverseBits ( Loop , NumberOfBits ) ;
    OutTable [ RIndex ] = InpTable [ Loop ] ;
  } /* for */ ;
  StageIndex = 2 ;
  while ( StageIndex <= TableLength )
  {
    Angle = Period / ( double ) StageIndex ;
    __real__ TwiddleFactor = cos ( Angle ) ;
    __imag__ TwiddleFactor = sin ( Angle ) ;
    Rotation =  ( 1.0F +  0.0iF ) ;
    for ( Index = 0 ; Index < StageIndex / 2 ; Index ++ )
    {
      for ( IntIndex = 0 ; IntIndex < TableLength ; IntIndex += StageIndex )
      {
        AddPartIndex = Index + IntIndex ;
        SubPartIndex = AddPartIndex + StageIndex / 2 ;
        AddPart = OutTable [ AddPartIndex ] ;
        SubPart = OutTable [ SubPartIndex ] * Rotation ;
        OutTable [ AddPartIndex ] = AddPart + SubPart ;
        OutTable [ SubPartIndex ] = AddPart - SubPart;
      } /* for */ ;
      Rotation = Rotation * TwiddleFactor ;
    } /* for */ ;
    StageIndex = 2 * StageIndex ;
  } /* while */ ;
} /* FourierFransformation */

Funkcja posiłkuje się pomocniczą funkcja do przeindeksowania tablic:

Kod: Zaznacz cały

static unsigned short ReverseBits ( unsigned short Index ,
                                    unsigned short NumBits )
{
  unsigned short Loop ;
  unsigned short Reverse ;
  /*----------------------------------------------------------------------*/
  Reverse = 0 ;
  for ( Loop = 0 ; Loop < NumBits ; Loop ++ )
  {
    Reverse = ( Reverse << 1 ) | ( Index & 0x01 ) ;
    Index = Index >> 1 ;
  } /* for */ ;
  return ( Reverse ) ;
} /* ReverseBits */

Prosty test to:

Kod: Zaznacz cały

int main ( int argc , char** argv )
{
  double SamplInt ;
  double TotalTime ;
  unsigned short Loop ;
  double SignalFreq ;
  double Freq ;
  double NyquistFreq ;
  /*----------------------------------------------------------------------*/
  OutFile = fopen ( "fft.txt" , "w" ) ;
  fprintf ( OutFile , "* * * Program obliczający transformate FFT * * *\n\nGeneracja sygnalu:\n" ) ;
  SamplInt = 0.002 ;
  InputSampleTable ( SamplInt , & TotalTime , TimeTable , InpSampleTable , SampleNum ) ;
  for ( Loop = 0 ; Loop < SampleNum ; Loop ++ )
  {
    fprintf ( OutFile , "Probka [ %d ] = %.8f\n" , Loop , __real__ InpSampleTable [ Loop ] ) ;
  } /* for */ ;
  SignalFreq = 1.0 / TotalTime ;
  NyquistFreq = 0.5 / SamplInt ;
  fprintf ( OutFile , "***********************************************************************\n\n" ) ;
  fprintf ( OutFile , "Okres probkowania: %.3f ms\n" , 1000.0 * SamplInt ) ;
  fprintf ( OutFile , "Okres calego sygnalu: %.3f ms\n" , 1000.0 * TotalTime ) ;
  fprintf ( OutFile , "Czestotliwosc probkowania: %.3f Hz\n" , 1.0 / SamplInt ) ;
  fprintf ( OutFile , "Czestotliwosc Nyquista: %.3f Hz\n" , NyquistFreq ) ;
  fprintf ( OutFile , "Czestotliwosc podstawowa (pierwsza harmoniczna): %.3f Hz\n" ,
            SignalFreq ) ;
  fprintf ( OutFile , "Pulsacja podstawowa (pierwsza harmoniczna): %.8f rad/s\n\n" ,
            DoublePi * SignalFreq ) ;
  FFT ( FFTTable , InpSampleTable , SampleNum ) ;
  for ( Loop = 0 ; Loop < SampleNum ; Loop ++ )
  {
    Freq = Loop * SignalFreq ;
    FreqTable [ Loop ] = Freq ;
    if ( Loop > SampleNumDiv2 )
      Freq = Freq - ( 1.0 / SamplInt ) ;
    MagnTable [ Loop ] = sqrt ( __real__ FFTTable [ Loop ] * __real__ FFTTable [ Loop ]  +
                                __imag__ FFTTable [ Loop ] * __imag__ FFTTable [ Loop ] ) /
                                ( double ) SampleNumDiv2 ;
    if ( Loop == 0 )
      MagnTable [ 0 ] /= 2.0 ;
    fprintf ( OutFile , "Harm. %d .re = %.10f  .im= %.10f  Ampl= %.10f  Czest. = %.3f Hz\n" , Loop , creal ( FFTTable [ Loop ] ) ,
                         cimag ( FFTTable [ Loop ] ) , MagnTable [ Loop ] , Freq ) ;
  } /* for */ ;
  IFFT ( IFFTTable , FFTTable , SampleNum ) ;
  for ( Loop = 0 ; Loop < SampleNum ; Loop ++ )
  {
    fprintf ( OutFile , "IFFT. %d .re = %.10f  .im= %.10f\n" , Loop ,
              creal ( IFFTTable [ Loop ] ) ,
          cimag ( IFFTTable [ Loop ] ) ) ;
  } /* for */ ;
  fprintf ( OutFile , "\n\nPorownanie probek\n\n" ) ;
  for ( Loop = 0 ; Loop < SampleNum ; Loop ++ )
  {
    fprintf ( OutFile , "prob. %d = %.10f [org.= %.10f, %.10f] roznica= %.10f, %.10f\n" ,
              Loop , __real__ IFFTTable [ Loop ] ,
          __real__ InpSampleTable [ Loop ] , __imag__ InpSampleTable [ Loop ] ,
             __real__ IFFTTable [ Loop ] - __real__ InpSampleTable [ Loop ] ,
             __imag__ IFFTTable [ Loop ] - __imag__ InpSampleTable [ Loop ] ) ;
  } /* for */ ;
  fclose ( OutFile ) ;
} /* main */

Uruchomiony program daje następujące wyniki:

Kod: Zaznacz cały

Okres probkowania: 2.000 ms
Okres calego sygnalu: 1024.000 ms
Czestotliwosc probkowania: 500.000 Hz
Czestotliwosc Nyquista: 250.000 Hz
Czestotliwosc podstawowa (pierwsza harmoniczna): 0.977 Hz
Pulsacja podstawowa (pierwsza harmoniczna): 6.13592315 rad/s

Harm. 0 .re = 5122.0457286895  .im= 0.0000000000  Ampl= 10.0039955638  Czest. = 0.000 Hz
Harm. 1 .re = 2.0489809989  .im= -0.2590904800  Ampl= 0.0080675658  Czest. = 0.977 Hz
Harm. 2 .re = 2.0588347310  .im= -0.5243270019  Ampl= 0.0082990295  Czest. = 1.953 Hz
Harm. 3 .re = 2.0755894067  .im= -0.8023542164  Ampl= 0.0086924759  Czest. = 2.930 Hz
Harm. 4 .re = 2.0997764969  .im= -1.1009014918  Ampl= 0.0092612281  Czest. = 3.906 Hz
Harm. 5 .re = 2.1322151735  .im= -1.4295714206  Ampl= 0.0100277447  Czest. = 4.883 Hz
Harm. 6 .re = 2.1741080435  .im= -1.8010140422  Ampl= 0.0110280829  Czest. = 5.859 Hz
Harm. 7 .re = 2.2272033940  .im= -2.2328175791  Ampl= 0.0123191936  Czest. = 6.836 Hz
Harm. 8 .re = 2.2940773487  .im= -2.7507669879  Ampl= 0.0139915255  Czest. = 7.813 Hz
Harm. 9 .re = 2.3786502676  .im= -3.3948605622  Ampl= 0.0161923629  Czest. = 8.789 Hz
Harm. 10 .re = 2.4872027530  .im= -4.2313298994  Ampl= 0.0191726177  Czest. = 9.766 Hz
Harm. 11 .re = 2.6305793444  .im= -5.3791227871  Ampl= 0.0233902224  Czest. = 10.742 Hz
Harm. 12 .re = 2.8296443492  .im= -7.0763319081  Ampl= 0.0297699719  Czest. = 11.719 Hz
Harm. 13 .re = 3.1315836918  .im= -9.8805871901  Ampl= 0.0404882049  Czest. = 12.695 Hz
Harm. 14 .re = 3.6753624829  .im= -15.4783708027  Ampl= 0.0621435455  Czest. = 13.672 Hz
Harm. 15 .re = 5.1589609325  .im= -32.4777012918  Ampl= 0.1284565997  Czest. = 14.648 Hz
Harm. 16 .re = -100.1183841374  .im= 1275.4552418006  Ampl= 4.9975729044  Czest. = 15.625 Hz
Harm. 17 .re = 0.2181536980  .im= 32.7422543536  Ampl= 0.1279022699  Czest. = 16.602 Hz
Harm. 18 .re = 1.7107868847  .im= 16.9996179574  Ampl= 0.0667401763  Czest. = 17.578 Hz
Harm. 19 .re = 2.4351963931  .im= 11.6381625941  Ampl= 0.0464461190  Czest. = 18.555 Hz
Harm. 20 .re = 3.0488426993  .im= 8.9200735631  Ampl= 0.0368231466  Czest. = 19.531 Hz
Harm. 21 .re = 3.7399756753  .im= 7.2645269296  Ampl= 0.0319168999  Czest. = 20.508 Hz
Harm. 22 .re = 4.6732568567  .im= 6.1369752974  Ampl= 0.0301317996  Czest. = 21.484 Hz
Harm. 23 .re = 6.1510284697  .im= 5.3001795109  Ampl= 0.0317169830  Czest. = 22.461 Hz
Harm. 24 .re = 9.0316361811  .im= 4.6137815972  Ampl= 0.0396166615  Czest. = 23.438 Hz
Harm. 25 .re = 17.5335590214  .im= 3.8867503605  Ampl= 0.0701530876  Czest. = 24.414 Hz
Harm. 26 .re = 1791.4388675185  .im= -50.0091362543  Ampl= 7.0005341787  Czest. = 25.391 Hz
Harm. 27 .re = -16.8180108582  .im= 4.2293530170  Ampl= 0.0677408306  Czest. = 26.367 Hz
Harm. 28 .re = -8.0840059685  .im= 3.7018777604  Ampl= 0.0347316045  Czest. = 27.344 Hz
Harm. 29 .re = -5.1858143742  .im= 3.3889336491  Ampl= 0.0241990665  Czest. = 28.320 Hz
Harm. 30 .re = -3.7374896825  .im= 3.1500906275  Ampl= 0.0190934927  Czest. = 29.297 Hz
Harm. 31 .re = -2.8682550837  .im= 2.9531163270  Ampl= 0.0160811271  Czest. = 30.273 Hz
Harm. 32 .re = -2.2885908441  .im= 2.7847944584  Ampl= 0.0140802450  Czest. = 31.250 Hz
Harm. 33 .re = -1.8745300378  .im= 2.6379496283  Ampl= 0.0126411954  Czest. = 32.227 Hz
Harm. 34 .re = -1.5640926186  .im= 2.5080332723  Ampl= 0.0115460032  Czest. = 33.203 Hz
Harm. 35 .re = -1.3228288976  .im= 2.3918857875  Ampl= 0.0106769996  Czest. = 34.180 Hz
Harm. 36 .re = -1.1300521273  .im= 2.2871822047  Ampl= 0.0099653179  Czest. = 35.156 Hz
Harm. 37 .re = -0.9725841363  .im= 2.1921441243  Ampl= 0.0093680115  Czest. = 36.133 Hz

Nie będę ukrywał, że parametry częstotliwościowe sygnału specjalnie spreparowałem. Wygenerowany sygnał „trafia próbkami” w częstotliwości, toteż na widmie widać, że prążki są ostre a nie rozmyte.

W realizacji pomocna była dosyć szeroka bibliografia, z której przede wszystkim wymienię:
  • Les Thede, „Practical Analog And Digital Filter Design”

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

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

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

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » wtorek 22 lut 2022, 14:03

Jest stronka poświęcona wolnym (na licencji Creative Common) książkom i innym publikacjom poświęcona problematyce DSP: https://www.dsprelated.com/freebooks.php
Sporo jest na temat przetwarzania dźwięku. Można sprawić, że procesory przemówią normalnym głosem. Również tej tematyce jest poświęcona popularna u nas książka "Sygnały cyfrowe, przetwarzanie i zastosowania" pod redakcją Alana V. Oppenheima.
Jak dokładniej ją obmacać (znaczy tą stronkę :) ), to można znaleźć znacznie więcej z nowoczesnych rozwiązań jak elektronika, mikrokontrolery, systemy wbudowane czy robotyka.

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

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

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » sobota 05 mar 2022, 01:29

Funkcje specjalne

W cyfrowym przetwarzaniu danych jest kilka „dziwnych” funkcji, które mają ogromne znaczenia w temacie. Do takich należy między innymi funkcja δ(t) – delta Diraca - realizuje niezwykłe działania i ma przedziwne własności. Każda funkcja matematyczna w jakimś stopniu daje się wyobrazić (w jaki sposób działa). W przypadku tej funkcji jest to trochę utrudnione. Jakkolwiek opis słowny jest w miarę prosty, to jej realizacja fizyczna już do prostych nie należy. Funkcja ta jest zdefiniowana następująco:
jeżeli argument funkcji jest różny od zera, to wartość funkcji wynosi zero, jeżeli argument funkcji ma wartość zero, to wartość funkcji wynosi plus nieskończoność.
Dodatkowo, funkcja ta ma tą własność, że
part21_form01.png

czyli iloczyn szerokości przedziału, gdzie funkcja ma wartość różną od zera (czyli wynoszącą zero, gdyż funkcja ma wartość zerową wszędzie z wyjątkiem argumentu zero) i wartości funkcji (wynoszącą plus nieskończoność) wynosi jeden. Patrząc na to wszystko z punktu widzenia logiki, nic tu nie trzyma się kupy, jakby to była jakaś kosmiczna pomyłka. Jednak nie jest to pomyłka.
Wykresem funkcji δ(t) jest wąska szpilka do plus nieskończoności.
Z funkcją tą wiąże się ciekawa historia. Otóż niejaki Dirac zajmując się swoją pracą rozwiązując równanie różniczkowe z zakresu fizyki kwantowej wprowadził właśnie tak zdefiniowaną funkcję. Okazała się bardzo przydatna. Nawet sama fizyka kwantowa wymyka się jakimkolwiek wyobrażeniom z fizyki klasycznej, tam wszystko jest inne. Jego praca była tryumfem rozumu, logicznego myślenia i dziwnej wyobraźni. Jego rozwiązanie przewidywało rzecz niepojętą: elektron ma ładunek elektryczny dodatni. Po raz pierwszy teoretycznie została przewidziana rzecz niezwykła, nierealna i niemożliwa (no bo jak elektron może mieć ładunek dodatni). Kilka lat później został eksperymentalnie odkryty. To jedynie utwierdziło przekonanie, że te zjawiska jakkolwiek wymykające się normalnemu wyobrażeniu, opisują całkiem prawdziwy świat.
Dlaczego wspominam o funkcji δ(t)?
A no właśnie transformata Fouriera (ta analityczna) dla funkcji sin i cos (dwóch najbardziej użytecznych) ma właśnie postać funkcji delty Diraca. No więc:
F ( sin ( a t) ) = j π [ δ ( ω + a ) - δ ( ω - a ) ]
oraz
F ( cos ( a t) ) = π [ δ ( ω + a ) + δ ( ω - a ) ]
Dowód matematyczny dla tych zależności jest wystarczająco dziwny, bym go tutaj przytaczał a zainteresowani zawsze mogą sięgnąć do źródeł.
Obrazkowo to wygląda następująco [dla sin(a t)]
sig110.png

oraz [dla cos (a t)]
sig111.png

Przyglądając się powyższym związkom, zrozumiałe staje się, że widmo zespolone składa się „z cieniutkich szpilek”. Dodatkowo można dostrzec jeszcze, że w przypadku funkcji sin, widmo jest nieparzyste i czysto urojone. Funkcja cos daje widmo parzyste oraz czysto rzeczywiste. Przy okazji wypłynął kolejny dziwny element: częstotliwość ujemna.
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

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

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

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » niedziela 06 mar 2022, 19:38

Najważniejsza

Ostatnio wspomniałem o funkcji δ(t). Wśród wielu różnych funkcji mających zastosowanie w DSP (można tu przykładowo wspomnieć o funkcji Heaviside), jest ona niezmiernie ważna (może nawet i najważniejsza, choć w gruncie rzeczy wszystkie są ważne). Jej wyjątkowość polega na tym, że transformata Fouriera dla funkcji δ(t) wynosi 1 (jeden). Jeżeli jakiś układ zostanie pobudzony impulsem będącym funkcją Diraca, to układ zareaguje na to jakąś odpowiedzią. Znajomość tej odpowiedzi jest tu decydująca, gdyż pozwala uzyskać odpowiedź na każdy inny sygnał wymuszający. Oczywiście teoria teorią a życie toczy się swoją drogą. W rzeczywistości uzyskanie pobudzenia impulsem o amplitudzie sięgającej nieskończoności jest praktycznie nieosiągalne. Dodatkowo, impuls ten musi trwać nieskończenie krótko. Pomimo, że taki impuls ma skończoną energią, jego uzyskanie w naszych warunkach jest niedostępne.
Można trochę „poluzować” swoje wymagania w stosunku do funkcji δ(t). Po pierwsze, niech czas trwania tego impulsu będzie wystarczająco krótki w stosunku do czasu reakcji układu na wymuszenie. To już wiele upraszcza. Drugie założenie, to niech jego amplituda wynosi 1 (jeden). W ten sposób funkcja δ(t) staje się funkcją do dystrybucji sygnału. Na rysunku 1 element δ(0) tworzy próbkę sygnału dla chwili t=0 (dokładniej dla numeru próbki wynoszącej zero) i pomnożony przez jego wartość daje sygnał dla t=0 (n=0): a0 δ(0) . Analogicznie próbka numer 1 to a1 δ(1), próbka numer n to an δ(n). Pokazuje to rysunek 1.
sig115.png

Dla dowolnego sygnału można zrealizować dekompozycję, czyli rozmontować na elementy składowe i potraktować odpowiedź układu na dowolny sygnał jako ileś odpowiedzi na impuls i połączyć je wszystkie do kupy. Niezbędna jest tu jedynie znajomość odpowiedzi impulsowej układu.
Do uzyskania odpowiedzi używana jest operacje splotu. Splot, jest to operacja matematyczna, która wiąże ze sobą dwa elementy dając w wyniku trzeci. Te dwa elementy to sygnał wymuszający i charakterystyka układu jako jego odpowiedź impulsowa. Wynikiem jest ich wspólny produkt: odpowiedź układu.
W matematyce splot oznacza się znakiem „*” (generuje to trochę nieporozumień, gdyż w językach programowania znak „*” jest używany do oznaczenia operacji mnożenia). No więc splotowa * to
part22_form01.png

Powyższa operacja jest operacją ciągłą (ma continuum czasowe). Oczywiście w sensie próbek koniecznością staje się przejście do świata dyskretyzacji. I zrobimy to w kolejnej części.
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

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

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

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » niedziela 06 mar 2022, 22:04

Tematyka związana z Cyfrowym przetwarzaniem sygnałów siłą rzeczy „dotyka” zagadnień numerycznych. To już samo w sobie wykracza poza ramy samego DSP i każdy z tych elementów może żyć własnym życiem. Jedynym elementem ich scalającym jest pewnego rodzaju podobieństwo zastosowania. Z tego powodu zdecydowałem się na utworzenie niezależnego wątku poświęconego szeroko pojętym zagadnieniom numerycznym <tutaj> (rozwiązywania określonych problemów matematycznych posiłkując się liczbami zmiennoprzecinkowymi i językiem programowania → w tym konkretnym przypadku językiem C).
W sumie wszystko opiera się o wiedzę uniwersalną, wypracowaną przez wysiłek wielu ludzi, chociaż w paru przypadkach kierowałem się własną koncepcją oraz potrzebą. Z pewnością niektóre rozwiązania można zrealizować inaczej (nie znaczy, że lepiej lub gorzej → po prostu inaczej).

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

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

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » poniedziałek 14 mar 2022, 14:08

Robimy splot

Operacja splotu pozwala określić sygnał wyjściowy. Jeżeli dany układ zostaje pobudzony sygnałem z zewnątrz, to, jak łatwo się domyślić, zostaje wygenerowany sygnał wyjściowy. Tu pierwszym pytaniem jakie się ciśnie, jest: co z tego wyjdzie?
Jak było wspomniane, jeżeli jest znana odpowiedź układu na pobudzenie impulsowe, to możliwe jest znalezienie odpowiedzi na każde pobudzenie. Sygnałem wejściowym jest spróbkowany jakiś sygnał, czyli ciąg próbek odzwierciedlający mierzony sygnał. To może być sygnał uzyskany z przetwornika ADC znajdującego się w mikrokontrolerze lub innego sygnału uzyskanego w jakikolwiek sposób. Do określenia odpowiedzi układu niezbędna jest znajomość odpowiedzi impulsowej. Na razie nie będziemy wnikać, w jaki sposób jest ona pozyskana. Odpowiedź impulsowa jest również ciągiem próbek. Jak pisałem wcześniej, do uzyskania rezultatu zastosowana jest operacja splotu. Sama matematyka mówi o splocie jako operacji ciągłej. W zastosowaniach przy sygnałach próbkowanych ta operacja ulega dyskretyzacji i operacja całkowania zamienia się na operację sumowania. W podręcznikach dotyczących cyfrowego przetwarzania sygnałów można znaleźć informacje pozwalające na implementację lub wręcz gotowe rozwiązania operacji splotu dla sygnałów spróbkowanych.
Operacja splotu dotyczy ciągu próbek. Z jednej strony są to próbki sygnału wejściowego z drugiej jest spróbkowana odpowiedź impulsowa układu (coś co charakteryzuje układ). Ciąg próbek sygnału zawiera określoną liczbę próbek, odpowiedź impulsowa również jest ciągiem próbek o określonej wielkości. Powstały sygnał jest sygnałem o długości będącej sumą wielkości sygnału wejściowego oraz odpowiedzi impulsowej. Przykładową operację splotu przedstawia poniży kawałek kodu:

Kod: Zaznacz cały

void Convolution ( double OutSignal [ ] ,
                   double InpSignal [ ] ,
                   short  InpSignalSize ,
                   double FiltSignal [ ] ,
                   short  FiltSignalSize )
{
  unsigned short SignalLoop ;
  unsigned short OutEndLoop ;
  unsigned short FiltLoop ;
  /*-----------------------------------------------------------------------*/
  OutEndLoop = InpSignalSize + FiltSignalSize ;
  for ( SignalLoop = 0 ; SignalLoop < OutEndLoop ; SignalLoop ++ )
    OutSignal [ SignalLoop ] = 0.0 ;
  for ( SignalLoop = 0 ; SignalLoop < InpSignalSize ; SignalLoop ++ )
  {
    for ( FiltLoop = 0 ; FiltLoop < FiltSignalSize ; FiltLoop ++ )
    {
      OutSignal [ SignalLoop + FiltLoop ] = OutSignal [ SignalLoop + FiltLoop ] +
                           InpSignal [ SignalLoop ] * FiltSignal [ FiltLoop ] ;
    } /* for */ ;
  } /* for */ ;
} /* Convolution */

Parametr OutSignal jest tablicą, która musi pomieścić wyliczony splot dwóch sygnałów, InpSignal jest sygnałem wejściowym (w programie tablicą liczb zmiennoprzecinkowych), InpSignalSize określa liczbę próbek wchodzących w skład sygnału wejściowego, FiltSignal jest tablicą zawierającą odpowiedź impulsową układu, której wielkość określa parametr FiltSignalSize. Jako przykład zastosowania splotu przedstawiony zostaje filtr sygnału. W ogólnym przypadku za pomocą odpowiedzi impulsowej możliwe jest uzyskanie każdego rezultatu.
Na początek rozpatrzmy sygnał wolnozmienny (o niskiej częstotliwości) zawierający zakłócenia (sygnał zaszumiony). Fragment sygnału (jakieś pół okresu) pokazuje rysunek 1.
sig120.png

Sygnał ten przechodzi przez filtr, którego odpowiedź impulsową pokazuje rysunek 2. Odpowiada to klasycznemu układowi całkującemu, a to jak wiadomo realizuje funkcję filtru dolnoprzepustowego, przepuszcza jedynie te sygnały, które zawierają odpowiednią częstotliwość.
Odpowiedź impulsowa filtru jest „wbita” w soft:

Kod: Zaznacz cały

void SetFilter2 ( double           Table [ ] ,
                  unsigned short * TableSize )
{
  unsigned short Index ;
  /*------------------------------------------------------------------------*/   
  Index = 0 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  Table [ Index ++ ] = 0.1625 ;
  * TableSize = Index ;
} /* SetFilter2 */

sig121.png

Odpowiedź impulsowa tego filtru odpowiada 32 próbkom o stałej wartości.
Sygnał to z kolei następujący kawałek softu:

Kod: Zaznacz cały

void InputSampleTable ( double         SampleTable [ ] ,
                        unsigned short Samples )
{
  unsigned short Loop ;   
  double Time ;
  /*-----------------------------------------------------------------------*/   
  for ( Loop = 0 ; Loop < Samples ; Loop ++ )
  {
    Time = 0.01 * ( double ) Loop ;
    SampleTable [ Loop ] = 10.0 + 5.0 * sin ( DoublePi * 0.6 * Time ) +
                                  3.0 * cos ( DoublePi * 45.4 * Time ) ;
  } /* for */ ;
} /* InputSampleTable */

By dostrzec pewne korelacje między sygnałem wejściowym i wyjściowym, oba są przedstawione w jednym układzie współrzędnych (rysunek 3).
sig122.png

Łatwo zauważyć, że sygnał wyjściowy ma znacząco większą amplitudę. Jest to działanie celowe, gdyż chcąc przedstawić oba przebiegi w jednym układzie współrzędnych i by one nie nakładały się, sygnał wyjściowy jest znacząco wzmocniony. By wzmocnienie wynosiło 0 dB, próbki odpowiedzi impulsowej powinny być na poziomie ok. 0.03, gdyż pole owego prostokąta odpowiada wzmocnieniu (32 próbki [od 0 do 31] pomnożone przez 0.03 daje około 1, ale wtedy przebiegi będą się nakładać). Zniekształcenia na samym początki i na samym końcu wynikają z operacji splotu, splot musi się „wypełnić” 32-oma próbkami sygnału wejściowego by działać poprawnie i zanim tego nie zrobi wszystko się kaszani. Podobnie jest na samym końcu, skończyły się próbki sygnału wejściowego, gdyż ten przykład filtru jest filtrem o skończonej odpowiedzi (w literaturze często nazywany filtrem FIR).
Niejako uzysk z filtru o prostokątnej charakterystyce odpowiedzi impulsowej jest całkiem, całkiem ale można uzyskać znacznie więcej. Niech odpowiedź filtru będzie jak na rysunku 4, jest to również filtr dolnoprzepustowy.
sig123.png

Jego odpowiedź jest następująca:

Kod: Zaznacz cały

void SetFilter1 ( double           Table [ ] ,
                  unsigned short * TableSize )
{
  unsigned short Index ;
  /*-----------------------------------------------------------------------*/
  Index = 0 ;
  Table [ Index ++ ] = 0.000000 ;
  Table [ Index ++ ] = 0.012889 ;
  Table [ Index ++ ] = 0.024889 ;
  Table [ Index ++ ] = 0.036000 ;
  Table [ Index ++ ] = 0.046222 ;
  Table [ Index ++ ] = 0.055555 ;
  Table [ Index ++ ] = 0.063999 ;
  Table [ Index ++ ] = 0.071555 ;
  Table [ Index ++ ] = 0.078221 ;
  Table [ Index ++ ] = 0.083999 ;
  Table [ Index ++ ] = 0.088888 ;
  Table [ Index ++ ] = 0.092888 ;
  Table [ Index ++ ] = 0.095999 ;
  Table [ Index ++ ] = 0.098221 ;
  Table [ Index ++ ] = 0.099555 ;
  Table [ Index ++ ] = 0.099999 ;
  Table [ Index ++ ] = 0.099555 ;
  Table [ Index ++ ] = 0.098221 ;
  Table [ Index ++ ] = 0.095999 ;
  Table [ Index ++ ] = 0.092888 ;
  Table [ Index ++ ] = 0.088888 ;
  Table [ Index ++ ] = 0.083999 ;
  Table [ Index ++ ] = 0.078221 ;
  Table [ Index ++ ] = 0.071555 ;
  Table [ Index ++ ] = 0.063999 ;
  Table [ Index ++ ] = 0.055555 ;
  Table [ Index ++ ] = 0.046222 ;
  Table [ Index ++ ] = 0.036000 ;
  Table [ Index ++ ] = 0.024889 ;
  Table [ Index ++ ] = 0.012889 ;
  Table [ Index ++ ] = 0.000000 ;
  * TableSize = Index ;
} /* SetFilter1 */

Dokładnie ten sam sygnał z identycznymi zakłóceniami poddany przetworzeniu przez układ o charakterystyce jak na rysunku 4 daje następujący wynik (rysunek 5).
sig124.png

Algorytm splotu oczyścił sygnał z wszelkich śmieci. Tu również zastosowane jest „przeskalowanie” by przebiegi nie nakładały się na siebie i było widać wykonaną przez filtr akcję.

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

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

Zealota
Posty: 10
Rejestracja: środa 08 maja 2019, 08:49

Re: Cyfrowe przetwarzanie sygnałów

Postautor: Zealota » wtorek 15 mar 2022, 20:54

Doskonała seria. Dzięki za opublikowanie.
Czy planujesz może coś o Q-constant transform?

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

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » niedziela 20 mar 2022, 02:16

Zealota pisze:Doskonała seria. Dzięki za opublikowanie.
Czy planujesz może coś o Q-constant transform?

Raczej nie planowałem pisać o przekształceniach z grupy wavelet, ale... może się kiedyś przydarzy.

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

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

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » piątek 25 mar 2022, 23:56

Szum i zakłócenia

Otaczająca rzeczywistość potrafi znacząco wpływać na sygnał. Nie trzeba nikogo przekonywać o prawdziwości tego stwierdzenia. Wystarczy sobie uświadomić (przypomnieć) choćby przydźwięk zasilania w urządzeniach audio. Źle wykonany (lub uszkodzony) zasilacz potrafi przenieść częstotliwość sieci na głośniki, a to już daje się stwierdzić własnymi zmysłami.
Zatem, w dziedzinie przetwarzania sygnałów, przydaje się coś co nazywamy generatorem liczb losowych. Przykład takiego generatora zawarłem <tutaj>.
Zobaczmy, co taki generator może. W programie wygenerowany ciąg losowy o wielkości 1024 próbek (jak raz dobry kawałek dla algorytmu FFT). Co on sobą reprezentuje?

Kod: Zaznacz cały

#define SampleNumDiv2       ( 512 )
#define SampleNum           ( 2 * SampleNumDiv2 )

(...)

typedef double complex SampleArrayType [ SampleNum ] ;
typedef double RealSampleArrayType [ SampleNum ] ;
static long idum ;

(...)

void InputSampleTable ( double         SamplInt ,
                        double *       TotalTime ,
                        double         TimeTable [ ] ,
                        double complex SampleTable [ ] ,
                        unsigned short Samples )
{
  unsigned short Loop ;   
  double Time ;
  /*--------------------------------------------------------------*/   
  for ( Loop = 0 ; Loop < Samples ; Loop ++ )
  {
    Time = SamplInt * ( double ) Loop ;
    TimeTable [ Loop ] = Time ;
    __real__ SampleTable [ Loop ] = Random1 ( & idum ) ;
    __imag__ SampleTable [ Loop ] = 0.0 ;
  } /* for */ ;
  * TotalTime = SampleNum * SamplInt ;
} /* InputSampleTable */


int main ( int argc , char** argv )
{
  double SamplInt ;

(...)

  SamplInt = 0.002 ;
  InputSampleTable ( SamplInt , & TotalTime , TimeTable ,
                     InpSampleTable , SampleNum ) ;

(...)

} /* main */

Ten kawałek softu oznacza, że co 0.002s do ciągu próbek (tablica InputSampleTable) dodawana jest losowa liczba. Taki ciąg ma długość SampleNum (=1024) sztuk. Wygenerowany sygnał ma następującą postać (właściwie to za każdym uruchomieniem powinien być inny, ale tego już nie badałem), rys 1. Jak można się tego spodziewać, wartości sygnału zawierają się w przedziale 0 .. 1.
sig132.png

Ważne pytanie, jakie należy zadać to: jakie jest widmo owego sygnału. Z odpowiedzią pośpieszy nam program.
Podstawowe informacje to:

Kod: Zaznacz cały

Okres probkowania: 2.000 ms
Okres calego sygnalu: 2048.000 ms
Czestotliwosc probkowania: 500.000 Hz
Czestotliwosc Nyquista: 250.000 Hz
Czestotliwosc podstawowa (pierwsza harmoniczna): 0.488 Hz
Pulsacja podstawowa (pierwsza harmoniczna): 3.06796158 rad/s

Z wydruku parametrów widma wyłaniają się ciekawe i ważne informacje. Po pierwsze:

Kod: Zaznacz cały

Harm. 0 .re = 512.2770795622  .im= 0.0000000000  Ampl= 0.5002705855

Amplituda zerowej harmonicznej to praktycznie 0,5. Wartość średnia funkcji przyjmującej wartości z zakresu 0 .. 1 wynosi 0,5, czyli „energia” sygnału jest właściwa, posiada chwilowe nadwyżki jak i chwilowe deficyty ale w sumie wszystko się równoważy.
Po drugie:

Kod: Zaznacz cały

Harm. 1 .re = 12.0981914367  .im= -2.2592057165  Ampl= 0.0240377440  Czest. = 0.488 Hz
Harm. 2 .re = -3.0232728487  .im= 3.3259582678  Ampl= 0.0087786781  Czest. = 0.977 Hz
Harm. 3 .re = 4.3317005201  .im= -3.7091115447  Ampl= 0.0111381460  Czest. = 1.465 Hz
Harm. 4 .re = 3.3308105715  .im= -3.6944050414  Ampl= 0.0097152858  Czest. = 1.953 Hz
Harm. 5 .re = 9.9106874861  .im= -8.0465006323  Ampl= 0.0249333752  Czest. = 2.441 Hz
Harm. 6 .re = 3.9768902366  .im= 2.4105340165  Ampl= 0.0090828356  Czest. = 2.930 Hz
Harm. 7 .re = -5.8611331299  .im= 12.5055613105  Ampl= 0.0269744838  Czest. = 3.418 Hz
Harm. 8 .re = -14.0153314555  .im= -11.6661197222  Ampl= 0.0356159113  Czest. = 3.906 Hz
Harm. 9 .re = 2.9338389371  .im= -1.4089765327  Ampl= 0.0063567020  Czest. = 4.395 Hz
Harm. 10 .re = 0.3480894650  .im= 7.1857995590  Ampl= 0.0140512218  Czest. = 4.883 Hz
Harm. 11 .re = 1.3588932537  .im= 0.8047107487  Ampl= 0.0030845467  Czest. = 5.371 Hz
Harm. 12 .re = 8.8169593886  .im= 1.0383533335  Ampl= 0.0173396311  Czest. = 5.859 Hz
Harm. 13 .re = -0.8290988902  .im= 1.7418951255  Ampl= 0.0037678629  Czest. = 6.348 Hz
Harm. 14 .re = -15.0503938630  .im= -2.4398981985  Ampl= 0.0297790695  Czest. = 6.836 Hz
Harm. 15 .re = 2.2791569589  .im= 0.1596729726  Ampl= 0.0044623892  Czest. = 7.324 Hz
Harm. 16 .re = -0.6430333279  .im= -3.3486345842  Ampl= 0.0066597970  Czest. = 7.813 Hz
Harm. 17 .re = 2.8051158794  .im= -6.2943514063  Ampl= 0.0134592187  Czest. = 8.301 Hz
Harm. 18 .re = 0.3250837642  .im= 3.8034228175  Ampl= 0.0074556450  Czest. = 8.789 Hz
Harm. 19 .re = 4.6159652506  .im= 11.0650987099  Ampl= 0.0234166203  Czest. = 9.277 Hz
Harm. 20 .re = 2.1242021107  .im= -2.9578497400  Ampl= 0.0071124622  Czest. = 9.766 Hz
Harm. 21 .re = -4.1391963778  .im= 8.9289748504  Ampl= 0.0192221179  Czest. = 10.254 Hz
Harm. 22 .re = -5.6822747947  .im= 1.6381423961  Ampl= 0.0115501804  Czest. = 10.742 Hz
Harm. 23 .re = -2.8813410297  .im= 1.0130701938  Ampl= 0.0059653302  Czest. = 11.230 Hz
Harm. 24 .re = -11.7397253132  .im= -8.9286692029  Ampl= 0.0288072553  Czest. = 11.719 Hz

Każda harmoniczna ma wartość niezerową (małe kilkadziesiąt wartości w powyższej ramce, by poznać więcej należy sięgnąć do archiwum → paczka na końcu).
Odnosząc wygenerowane w programie wartości do pisma obrazkowego, oba widma (rzeczywiste oraz urojone) są następujące (rys 2 i 3).
sig131.png

sig130.png

Jak się zagłębić w lekturę wygenerowanych wyników, to można wyciągnąć parę istotnych elementów. Przykładowo:

Kod: Zaznacz cały

Harm. 103 .re = -0.4790450039  .im= 3.7839146216  Ampl= 0.0074494487  Czest. = 50.293 Hz
Harm. 104 .re = -0.9806061082  .im= 0.1494102364  Ampl= 0.0019373501  Czest. = 50.781 Hz

Harm. 320 .re = -2.4781356587  .im= -2.0299731674  Ampl= 0.0062566942  Czest. = 156.250 Hz

Jest to kwestia interpretacji. Czy należy uznać wartości rzędu kilku tysięcznych za wartość zerową (czyli dostatecznie bliską zera)? Takich miejsc jest więcej.
Jeżeli tak, to widmo te ma kilka dziur i mówiąc kolokwialnie jest trochę szczerbate.
Dla tropicieli.
fft.7z
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

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

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

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » poniedziałek 28 mar 2022, 23:56

Zastosowanie odpowiedzi impulsowej w medycynie:
Siedzi facet z bolącym zębem na fotelu u dentysty.
Sympatyczna pani doktór chcąc ustalić szczegóły związane
z leczeniem pacjenta, uderza lekko posiadanym
w ręku narzędziem w ząb.
I to jest funkcja delta (impuls Diraca).
W odpowiedzi facet o mało co... i przypir...ił lekarce w …
(ale się powstrzymał). I to jest odpowiedź
układu na impuls Diraca.


Odpowiedź impulsowa

Rozpatrzmy układ elektryczny jak na rysunku 1. Jak wiadomo, pełni on rolę filtru dolnoprzepustowego.
sig138.png

Teraz na jego wejściu przyłożone jest wymuszenie w postaci impulsu Diraca. W obwodzie popłynie prąd i takie tam... Można zapisać następujące równanie:
part25_form01.png

ponieważ
part25_form02.png

to:
part25_form03.png

Stosując do równania obustronnie transformatę Laplace'a, mamy (transformata Laplace'a funkcji delta Diraca wynosi 1).
part25_form04.png

Stosując transformatę odwrotną do powyższej zależności powstaje:
part25_form05.png

No to należy ją spróbkować. Kawałek programu generuje odpowiednie ciągi próbek:

Kod: Zaznacz cały

void SetFilter ( double           Table [ ] ,
                 unsigned short * TableSize )
{
  unsigned short Loop ;
  /*-----------------------------------------------------------------*/   
  * TableSize = 45 ;
  for ( Loop = 0 ; Loop < * TableSize ; Loop ++ )
  {
    Table [ Loop ] = 0.2 * exp ( - 0.1 * Loop ) ;
  } /* for */ ;
} /* SetFilter */

Tu oczywiście nie trzymałem się jakichś konkretnych wartości parametrów zastosowanych komponentów (bardziej tu chodzi o zrozumienie filozofii), ale można to zrobić.
Podobnie został wygenerowany w miarę wolnozmienny przebieg sin, na który został nałożony sygnał zakłócający generowany losowo (o ile ta losowość jest losowa, ale z pewnością jest elementem niepożądanym), który należy usunąć.

Kod: Zaznacz cały

void InputSampleTable ( double         SampleTable [ ] ,
                        unsigned short Samples )
{
  unsigned short Loop ;   
  double Time ;
  /*-----------------------------------------------------------------*/   
  for ( Loop = 0 ; Loop < Samples ; Loop ++ )
  {
    Time = 0.01 * ( double ) Loop ;
    SampleTable [ Loop ] = 10.0 + 5.0 * sin ( DoublePi * 0.6 * Time ) +
                           3.0 * Random1 ( & idum ) ;
  } /* for */ ;
} /* InputSampleTable */

Oba te sygnały zostały splecione (obliczony splot sygnałów), co odpowiada w rzeczywistości przepuszczeniu sygnału analogowego przez filtr RC.
Sygnał, odpowiedź impulsowa i wynik splotu pokazują rysunki 2, 3 i 4.
sig135.png

sig136.png

sig137.png

No muszę przyznać, że wynik filtracji mnie nie powalił. Z drugiej strony to jest filtr jedynie z jednym biegunem, więc również nie należy się niczemu dziwić. Z punktu widzenia procesu DSP, to można dowolnie kształtować charakterystyki filtrów i w ten sposób uzyskać niejako cuda na kiju. To w końcu jedynie matematyka, która czasami (za przeproszeniem) jest oderwana od rzeczywistości. Jednak zanim, to wróćmy do normalności i popatrzmy, co jeszcze ma w zanadrzu prawdziwa elektrotechnika.
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

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

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

Re: Cyfrowe przetwarzanie sygnałów

Postautor: gaweł » wtorek 05 kwie 2022, 00:31

Różniczkowanie sygnału

Zastosowanie operacji splotu służy nie tylko do filtrowania sygnałów. Wszystko zależy od charakterystyki odpowiedzi impulsowej. To ona decyduje o efektach końcowych. Nie zawsze cechy tej odpowiedzi muszą odzwierciedlać określoną realizację układową. Tu panuje szeroka dowolność i można uzyskać każdą możliwość. Jako przykład niech będzie jakiś układ, którego odpowiedź impulsowa jest następująca:
  • próbka o numerze 0 ma wartość 1,
  • próbka o numerze 1 ma wartość -1.
Takie coś pozwala uzyskać na wyjściu pochodną sygnału wejściowego. Jak łatwo się zorientować, nie należy to do zbyt skomplikowanych opisów. W programie to wygląda następująco:

Kod: Zaznacz cały

void SetFilter ( double           Table [ ] ,
                 unsigned short * TableSize )
{
  unsigned short Index ;
  /*------------------------------------------------------------------------*/   
  Index = 0 ;
  Table [ Index ++ ] = 25.0 ;
  Table [ Index ++ ] = -25.0 ;
  * TableSize = Index ;
} /* SetFilter */

Dałem trochę więcej (25.0 zamiast 1.0) gdyż to jednocześnie wzmacnia sygnał. Oprócz zwiększenia amplitudy nic nie wnosi do sygnału. Sam sygnał w programie wygenerowałem możliwie najprostszy, by można było gołym, okiem zaobserwować zachodzące zdarzenia.

Kod: Zaznacz cały

void InputSampleTable ( double         SampleTable [ ] ,
                        unsigned short Samples )
{
  unsigned short Loop ;   
  double Time ;
  /*----------------------------------------------------------------------*/   
  for ( Loop = 0 ; Loop < Samples ; Loop ++ )
  {
    Time = 0.01 * ( double ) Loop ;
    SampleTable [ Loop ] = sin ( DoublePi * 0.6 * Time ) ;
  } /* for */ ;
} /* InputSampleTable */

Algorytm operacji splotu jest już wcześniej wystarczająco dokładnie przedstawiony. Umieszczając sygnał wejściowy oraz odpowiedź układu na ten sygnał w jednym układzie współrzędnych można się o tym przekonać naocznie.
sig140.png

Miejsca, gdzie jeden z sygnałów przechodzi przez zero, to drugi osiąga ekstremum, no nic innego jak tylko pochodna.
I nie należy się tu niczemu dziwić, nie ma w tym żadnej magii. Dwupunktowy sygnał odpowiedzi impulsowej oblicza różnicą pomiędzy dwoma kolejnymi próbkami sygnału wejściowego. Biorąc pod uwagę, że przyrost argumentu jest stały (wynoszący jeden, bo to kolejne próbki), to wynik będzie proporcjonalny do dyskretnej pochodnej sygnału (jakkolwiek niezwykle to brzmi).
Nie masz wymaganych uprawnień, aby zobaczyć pliki załączone do tego posta.

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


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

Kto jest online

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