rootuss

Juliusz

Nov 4th, 2019
221
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. % Autor opracowania: Juliusz Tarnowski
  2. % Proszę o nierozpowszechnianie bez zgody autora.
  3. % W przypadku pytań lub uwag proszę o wiadomość na messengerze.
  4. % Opracowanie zawiera materiał niezbędny do zaliczenia pierwszego kolokwium z laboratoriów PPS.
  5. % Jest to jednocześnie skrypt przystosowany do uruchomienia w Matlabie.
  6. % W Octave przed uruchomieniem należy w linii poleceń wpisać: pkg load signal, aby załadować odpowiednią bibliotekę.
  7. % Jeśli pojawi się błąd, należy ją doinstalować.
  8.  
  9. %% 1) Komendy ogólne
  10. clear all;            % Służy do wyczyszczenia workspace'a, czyli wszystkich zdefiniowanych zmiennych
  11.                       % Należy używać go na początku nowego skryptu, lub na początku "rozdziału" skryptu
  12. close all;            % Służy do zamknięcia wszystkich okienek – inaczej trzeba zamykać je ręcznie.
  13.  
  14.  
  15. %% 2) Definiowanie osi czasu.
  16. % Parametry:
  17. t_max = 5;            % Czas trwania sygnału. Wyrażony w sekundach.
  18. fp = 2000;            % Częstotliwość próbkowania, ilość próbek na sekundę. Wyrażona w hercach.
  19. N = 10001;            % Ilość próbek w całym sygnale. Należy pamiętać o próbce "numer 0".
  20.  
  21. % Wystarczy zdefiniować dwie spośród powyższych zmiennych. Trzecią można wyznaczyć:
  22. t_max = (N-1)/fp;     % Czas trwania sygnału.
  23. fp = (N-1)/t_max;     % Częstotliwość próbkowania.
  24. N = (t_max * fp) + 1; % Liczba próbek.
  25.  
  26.  
  27. %% 3) Metody generowania osi.
  28. % Na PPSach definiowane są dwie podstawowe osie: czasu oraz częstotliwości.
  29. % Pojawiają się też inne, np. oś przesunięcia (wyrażona w sekundach lub ilości próbek).
  30. % Przy tworzeniu każdego sygnału, widma czy charakterystyki należy zastanowić się, czy została zdefiniowana odpowiednia oś.
  31. % Należy wybrać taką metodę, jaka jest dla Ciebie najwygodniejsza, ale warto pamiętać o wszystkich.
  32.  
  33. t = 0 : 1/fp : (N-1)/fp;   % Podstawowa metoda generowania osi czasu.
  34. %   ^   ^^^^   ^^^^^^^^^   % Komentarze czytaj od dołu do góry. ;)
  35. %   |    |          \ ... do ostatniej chwili w sygnale.
  36. %   |    |     ^^^^  ^^^
  37. %   |    |      |     \ (Częstotliwość próbkowania)
  38. %   |    |       \ (Ostatni indeks próbki - jedna została zużyta na "chwilę zero".)
  39. %   |     \ ... co jeden okres próbkowania...
  40. %    \ Tablica wartości od zera...
  41.  
  42. t = 0 : 1/fp : t_max;      % To samo, co powyżej, ale z wyliczonym maksymalnym czasem.
  43.  
  44. t = (0:N-1) ./ fp;         % Najbardziej zwięzła forma generowania osi czasu. Niewydajna obliczeniowo.
  45. %   ^^^^^^^ ^^ ^^
  46. %      |    |   \ ... przez częstotliwość próbkowania.
  47. %      |     \ ... podzielony... (Kropka oznacza, że wykonujemy operację na próbkach wektora, a nie na całym wektorze.)
  48. %       \ Wektor z liczbami naturalnymi od 0 do N-1 wartości...
  49.  
  50. t = linspace(0,t_max,N);   % Funkcja linespace() najczęściej jest wykorzystywana do generowania osi częstotliwości.
  51. %   ^^^^^^^^ ^ ^^^^^ ^
  52. %      |     |   |   \ Ilość próbek.
  53. %      |     |    \ Wartość maksymalna.
  54. %      |      \ Wartość początkowa. W tym przypadku to chwila "0", ale może to być np. chwila rozpoczęcia się fonemu, lub przesunięcie.
  55. %       \ Nazwa funkcji. Generuje dokładnie N równomiernie rozłożonych wartości, między 0 a t_max.
  56.  
  57.  
  58. %% 4) Najważniejsze rodzaje sygnałów
  59. % Definiując parametry sygnału warto dodawać indeks (np. "_s"), aby później nie było kolizji oznaczeń.
  60.  
  61. %% 4.1) Szum biały, gaussowki, normalny.
  62. % Podstawowy rodzaj szumu. W razie wątpliwości użyj właśnie jego.
  63. W_1 = 3; % Wariancja/sigma. Pełni rolę pseudoamplitudy. Zgodnie z prawem trzech sigm 99.8% wartości będzie leżało nie dalej, niż w odległości 3*W_1 od wartości średniej.
  64. sr_1 = 9; % Wartość średnia sygnału.
  65.  
  66. s_1 = W_1*randn(1,N) + sr_1;
  67. %     ^^^ ^^^^^ ^^^    ^^^^
  68. %      |    |    |       \ Wartość średnia.
  69. %      |    |    \ Argumenty funkcji: 1 wiersz, N kolumn. Pomijając pierwszy argument wygenerujemy macierz rozmiaru N x N.
  70. %      |     \ Funkcja generująca losowe wartości o rozkładzie normalnym. Posiada literkę "n".
  71. %       \ Wariancja.
  72.  
  73. %% 4.2) Szum jednorodny.
  74. % Dodatkowy rodzaj szumu, służący najczęściej do symulowania błędu kwantyzacji.
  75. A_2 = 3; % Pseudoamplituda. Odległość między skrajnymi wartościami szumu.
  76. sr_2 = 6; % Wartość średnia.
  77.  
  78. s_2 = A_2*(rand(1,N)-1/2) + sr_2; %
  79. %     ^^^  ^^^^ ^^^ ^^^^    ^^^^
  80. %      |     |   |   |        \ Wartość średnia.
  81. %      |     |   |    \ Unormowanie wartości średniej szumu.
  82. %      |     |   |      (Standardowo rand() ma wartość oczekiwaną równą 1/2, a wygodniej jest, przesuniemy ją do 0.)
  83. %      |     |   \ Argumenty funkcji: 1 wiersz, N kolumn.
  84. %      |      \ Funkcja generująca losowe wartości o rozkładzie jednorodnym [0,1]. Nie posiada literki "n".
  85. %       \ Pseudoamplituda.
  86.  
  87. %% 4.3) Sinusoida.
  88. % Parametry sinusoid:
  89. A_3 = 5; % Amplituda, czyli wychylenie od zera.
  90. f_3 = 24; % Częstotliwość sinusoidy. Musi być mniejsza od połowy częstotliwości próbkowania, inaczej pojawi się aliasing.
  91. fi_3 = pi/4; % Przesunięcie fazowe. Wartość wyrażona w radianach, czyli np. [0,2*pi] lub [-pi,pi]
  92.  
  93. s_3 = A_3*cos(t*2*pi*f_3+fi_3);
  94. %     ^^^ ^^^ ^ ^^^^ ^^^ ^^^^
  95. %      |   |  |  |    |    \ Przesunięcie fazowe.
  96. %      |   |  |  |     \ Częstotliwość.
  97. %      |   |  |   \ Przemnożenie osi czasu przez 2*pi, żeby unormować częstotliwość (móc wyrażać ją w wygodnych dla człowieka hercach, a nie radianach).
  98. %      |   |   \ Oś czasu.
  99. %      |    \ Jeśli prowadzący nie zwraca na to uwagi, lepiej używać funkcji cos(), a nie sin().
  100. %      |      Wynika to z własności transformat oraz ze wzoru Eulera: e^(i*x) = cos(x) + i*sin(x)
  101. %       \ Amplituda.
  102.  
  103. % Sinusoida o zmiennej częstotliwości.
  104. % Ten sygnał na PPS jest najczęściej wykorzystywany do testowania filtrów.
  105. % Zachęcam posłuchać, jak brzmi ten sygnał na stronie: https://en.wikipedia.org/wiki/Chirp
  106. % Parametry:
  107. f0 = 1;                      % Częstotliwość dla czasu 0s.
  108. t1 = t_max;                  % Koniec osi czasu.
  109. f1 = 50;                     % Częstotliwość dla czasu t1.
  110.  
  111. s_4 = chirp(t, f0, t1, f1);  % Sinusoida o liniowo zmieniającej się częstotliwości.
  112. %     ^^^^^ ^  ^^  ^^  ^^
  113. %       |   |  |   |    \ Częstotliwość w zadanej chwili.
  114. %       |   |  |    \ Wybrana w dowolny sposób chwila. Najwygodniej podać czas trwania sygnału (ostatni moment).
  115. %       |   |   \ Częstotliwość w chwili 0.
  116. %       |    \ Oś czasu. Jako jedyna jest wektorem, wszystkie pozostałe zmienne są liczbami.
  117. %        \ Nazwa funkcji. Przyjmuje 4 argumenty.
  118.  
  119.  
  120. %% 5) Wyświetlanie sygnałów
  121. % Funkcje "porządkujące":
  122. figure(1); % Wybór okienka.
  123. subplot(2,2,1); % Wybór komórki.
  124. %       ^ ^ ^
  125. %       | | \ ... pierwsza komórka.
  126. %       |  \ ... dwie kolumny...
  127. %        \ Dwa wiersze...
  128.  
  129. % Funkcje "wyświetlające":
  130. plot(t,s_1,'g--');               % Standardowy wykres, generowany linią ciągłą.
  131. %    ^ ^^^ ^^^^^
  132. %    |  |   \ Ozdobniki, czyli kolor i kształt wykresu.
  133. %    |   \ Wartości sygnału. x oraz t muszą być takiego samego wymiaru.
  134. %     \ Oś czasu.
  135.  
  136. % Dostępne kolory: 'r', 'g', 'b', 'c', 'm', 'y', 'k'.
  137. % Dostępne ozgobniki: '-', '--', '.', 'o', '*' i inne.
  138.  
  139. subplot(222);
  140. plot(t,s_3,'b', t,s_4,'r');    % Funkcja plot() przewiduje możliwość dodania większej ilości wykresów.
  141.  
  142. subplot(223);
  143. stem(t,s_2,'ko');              % Wykres słupkowy, najczęściej wykorzystywany przy wizualizacji widma.
  144.  
  145. subplot(224);
  146. nbits = 50;                    % Ilość słupków histogramu.
  147. hist(s_1,nbits);               % Histogram, czyli wizualizacja zbioru wartości sygnału.
  148. %    ^^^ ^^^^^
  149. %     |    \ Drugi argument jest opcjonalny. Jeśli go nie podamy, Matlab automatycznie dobierze odpowiednią wartość.
  150. %      \ Analizowany sygnał.
  151.  
  152. % Komendy specjalne:
  153. % Poniższe komendy należy podać poniżej funkcji plot(). Wywołanie jej po subplot(), bez wywołania plot() nic nie zmieni.
  154. hold on;  % Komenda pozwalająca na wyświetlanie wielu wykresów na jednej komórce.
  155. hold off; % Antagonista powyższej komendy, używany znacznie rzadziej.
  156. grid on;  % Aktywowanie siatki na wykresie. Ją również można dezaktywować.
  157.  
  158. % Funkcje generujące opis:
  159. title('Tytuł wykresu, napisany pogrubioną czcionką, na górze.')
  160. xlabel('opis osi x (zwykle to czas / częstotliwość)')
  161. ylabel('opis osi y (zwykle to sygnał / widmo)')
  162.  
  163.  
  164. %% 6) Wczytywanie nagrania
  165. % Wczytanie nagrania wykorzystuje funkcję wavread(). Alternatywą dla niej jest audioread. Działają niemal identycznie.
  166. [x, fpx] = wavread('mbi04wyczy.wav'); % Wczytanie sygnału i jego częstotliwości próbkowania.
  167. %^  ^^^    ^^^^^^^ ^^^^^^^^^^^^^^^^
  168. %|   |        |            \ Nazwa pliku (musi być w tym samym folderze co skrypt).
  169. %|   |         \ Funkcja pozwalająca na wczytanie pliku. Można zastąpić ją funkcją audioread().
  170. %|    \ Zmienna, do której przypiszemy częstotliwość próbkowania nagrania z pliku.
  171. % \ Zmienna, do której przypiszemy sygnał z pliku.
  172. Nx = length(x);                       % Ustalanie ilości próbek nagrania.
  173. tx_max = (Nx-1)/fpx;                  % Czas trwania sygnału.
  174. tx = 0 : 1/fpx : (Nx-1)/fpx;          % Generowanie osi czasu.
  175.  
  176.  
  177. %% 7) Decymacja
  178. % Decymacja to ograniczenie rozmiaru sygnału, poprzez usunięcie nadmiarowych próbek.
  179. % Można to zrealizować na dwa sposoby: ręcznie, oraz z wykorzystaniem funkcji decimate().
  180. % Decymacja ręczna nie jest kompletna, ponieważ pozwala na występowanie aliasingu.
  181. % Funkcja decimate() posiada specjalny filtr antyaliasingowy.
  182.  
  183. dr = 6;                               % Rząd decymacji. Mogą pojawić się takie implementacje, gdzie musi być to potęga dwójki.
  184. y_1 = x(1:dr:Nx);                     % Decymacja ręczna. W praktyce nie wykorzystywana.
  185. %     ^ ^ ^^ ^^
  186. %     | | |   \ (Kończymy na ostatniej.)
  187. %     | |  \ ... co dr-tą próbkę.
  188. %     |  \ (Zaczynamy od pierwszej.)
  189. %      \ Oryginalny sygnał z którego pozostawiamy...
  190.  
  191. y_2 = decimate(x, dr);                % Decymacja sygnału z użyciem odpowiedniej funkcji.
  192. %              ^  ^^
  193. %              |   \ Rząd decymacji.
  194. %               \ Sygnał do decymacji.
  195.  
  196. % Po decymacji należy wygenerować nową oś czasu:
  197. Ny = length(y_2);                     % Ustalanie ilości próbek nagrania.
  198. fpy = fpx / dr;                       % Nowa częstotliwość.
  199. ty = 0 : 1/fpy : (Ny-1)/fpy;          % Generowanie osi czasu.
  200.  
  201.  
  202. %% 8) Analiza częstotliwości
  203. % Do analizy częstotliwości używamy funkcji fft().
  204. % Wymaga ona podania jako argumentu liczby próbek widma, będących potęgą dwójki.
  205. % Wynik będzie symetryczny, dlatego wyświetlać będziemy tylko połowę widma.
  206. % Transformata Fouriera generuje wartości zespolone, więc przed narysowaniem wykresu należy przeliczyć je na liczby rzeczywiste.
  207.  
  208. Nfx = 2 ^ nextpow2(Nx);           % Liczba próbek widma.
  209. %     ^   ^^^^^^^^ ^^
  210. %     |      |      \ Ilość próbek sygnału.
  211. %     |       \ Funkcja zwracająca wykładnik (ilość cyfr w systemie binarnym) pierwszej potęgi dwójki, która jest większa bądź równa zadanej liczbie.
  212. %     |         (Na przykład dla wartości od 17 do 32 funkcja zwróci wartość 5, ponieważ 32=2^5.)
  213. %      \ Dwójka, którą podniesiemy do potęgi wygenerowanej przez funkcję.
  214.  
  215. Nfx21 = Nfx/2 + 1;                % Zmienna pomocnicza, służąca wyświetlaniu tylko połowy widma (druga połowa zawiera takie same informacje).
  216.  
  217. fx = linspace(0, fpx/2, Nfx21);   % Oś częstotliwości.
  218. %    ^^^^^^^^ ^  ^^^^^  ^^^^^
  219. %       |     |    |      \ Ilość próbek, jakie chcemy wyświetlać.
  220. %       |     |     \ Częstotliwość próbkowania (cz. Nyquista), czyli maksymalna częstotliwość, jaka interesuje nas w analizowanym widmie.
  221. %       |      \ Zero, czyli minimalna interesująca nas częstotliwość.
  222. %        \ Funkcja zwracająca tablicę wartości: Nf21 równie oddalonych próbek, których wartości zaczynają od 0, a kończą się na fp/2.
  223.  
  224. vx = fft(x, Nfx);          % Wyznaczenie widma zespolonego.
  225. %    ^^^ ^  ^^^
  226. %     |  |   \ ... oraz liczbę próbek widma, będącą potęgą dwójki.
  227. %     |   \ ... sygnał...
  228. %      \ Szybka transformata Fouriera, która jako argument przyjmuje...
  229.  
  230. avx = abs(vx);              % Widmo amplitudowe, czyli moduł transformaty Fouriera.
  231. lavx = 20*log10(avx);       % Widmo amplitudowe w skali decybelowej / logarytmicznej.
  232. %      ^^ ^^^^^ ^^^
  233. %      |    |    \ Widmo, które przeliczamy na skalę decybelową.
  234. %      |     \ ... na dekadę.
  235. %       \ Dwadzieścia decybeli...
  236. fvx = angle(vx);            % Widmo fazowe, czyli argument wartości transformaty Fouriera. Przyjmuje wartości ze skali [-pi, pi].
  237. ufvx = unwrap(fvx);         % "Rozwinięte" widmo fazowe. Funkcja unwrap() nie pozwala na gwałtowne skoki wartości
  238.                             % (Zamiast przeskoku na granicy wartości, będziemy kontynuować wykres.
  239.                             % ... to trochę tak, jakbyśmy na zegarku zamiast wyświetlać po 23:59 godzinę 00:00, wyświetlilibyśmy 24:00.)
  240.  
  241. % Wizualizacja analizy widma:
  242. figure(2);
  243. subplot(311);
  244. plot(tx,x);
  245. title('Przebieg sygnału')
  246. xlabel('Czas [s]');
  247. ylabel('Wartość sygnału');
  248.  
  249. subplot(323);
  250. plot(fx,avx(1:Nfx21));
  251. title('Widmo amplitudowe')
  252. xlabel('Częstotliwość [Hz]');
  253.  
  254. subplot(324);
  255. plot(fx,fvx(1:Nfx21));
  256. title('Widmo fazowe')
  257. xlabel('Częstotliwość [Hz]');
  258.  
  259. subplot(325);
  260. plot(fx,lavx(1:Nfx21));
  261. title('Widmo amplitudowe w skali decybelowej')
  262. xlabel('Częstotliwość [Hz]');
  263.  
  264. subplot(326);
  265. plot(fx,ufvx(1:Nfx21));
  266. title('Rozprostowane widmo fazowe')
  267. xlabel('Częstotliwość [Hz]');
  268.  
  269.  
  270. %% 9) Wycinanie fragmentu nagrania.
  271. % Wycięcie sylaby "wy".
  272. t0_wy = 0.11;                     % Chwila w której zaczyna się sylaba. Odczytujemy ją z wykresu.
  273. t1_wy = 0.26;                     % Chwila w której kończy się sylaba.
  274. pr0_wy = round(fpx*t0_wy);        % Numer próbki w której zaczyna się sylaba.
  275. %   ^^   ^^^^^ ^^^ ^^^^^
  276. %   |      |    |    \ ... przemnożona przez chwilę w której sylaba się zaczyna.
  277. %   |      |     \ Częstotliwość próbkowania...
  278. %   |       \ Zaokrąglenie. W przypadku nagrania nie ma znaczenia, czy w górę, czy w dół.
  279. %   |         Dla sygnałów o niższej częstotliwości próbkowania warto użyć funkcji ceil() lub floor().
  280. %    \ Polecam pamiętać o indeksach. Lub zmienić język programowania, np. na Pythona z NumPy. :P
  281. pr1_wy = round(fpx*t1_wy);        % Numer próbki w której kończy się sylaba.
  282. Nwy = pr0_wy - pr1_wy;            % Długość sylaby w próbkach.
  283. wy = x(pr0_wy:pr1_wy);            % Wycięcie sylaby z nagrania.
  284. twy = t0_wy : 1/fpx : t1_wy;      % Oś czasu dla sylaby.
  285. %     ^^^^^   ^^^^^   ^^^^^
  286. %       |       |       \ ... kończymy na chwili końcowej.
  287. %       |        \ ... (nie zmieniając czêstotliwości próbkowania, bo sylaba jest fragmentem nagrania)...
  288. %        \ Tym razem nie zaczynamy od zera, a od chwili początkowej sylaby...
  289.  
  290. t0_czy = 0.36;                    % Chwila w której zaczyna się sylaba. Odczytujemy ją z wykresu.
  291. t1_czy = 0.5;                     % Chwila w której kończy się sylaba.
  292. pr0_czy = round(fpx*t0_czy);      % Numer próbki w której zaczyna się sylaba.
  293. pr1_czy = round(fpx*t1_czy);      % Numer próbki w której kończy się sylaba.
  294. Nczy = pr0_czy - pr1_czy;         % Długość sylaby w próbkach.
  295. czy = x(pr0_czy:pr1_czy);         % Wycięcie sylaby z nagrania.
  296. tczy = t0_czy : 1/fpx : t1_czy;   % Oś czasu dla sylaby.
  297.  
  298. figure(3);
  299. subplot(411);
  300. plot(tx,x);
  301. title('Słowo "wyczyść"')
  302. xlabel('czas [s]');
  303.  
  304. subplot(423);
  305. plot(twy,wy);
  306. title('Sylaba "wy"')
  307. xlabel('czas [s]');
  308.  
  309. subplot(424);
  310. plot(tczy,czy);
  311. title('Sylaba "czy"')
  312. xlabel('czas [s]');
  313.  
  314.  
  315. %% 10) Korelacja, autokorelacja
  316. % Aby obliczyć korelację, należy podać ilość próbek, o jaką ma zostać przesunięty sygnał.
  317. % Może być zadana na dwa sposoby:
  318. kmax = 230;                      % Jako liczba naturalna.
  319.  
  320. tk_max = 0.15;                   % Lub jako czas, w tym przypadku 0.15 sekundy.
  321. kmax = tk_max * fpx;             % Wyznaczenie ilości próbek.
  322. %      ^^^^^^   ^^^
  323. %        |       \ ... przemnożone przez ilość próbek na sekundę.
  324. %         \ Maksymalne opóźnienie...
  325.  
  326. r = xcorr(wy,czy,kmax);          % Obliczanie autokorelacji.
  327. %   ^^^^^ ^^ ^^^ ^^^^
  328. %     |   |   |     \ ... ilość próbek, definiujące maksymalne opóźnienie.
  329. %     |    \ ... dwukrotnie ten sam sygnał, oraz...
  330. %      \ Funkcja do wyznaczania autokorelacji, która przyjmuje jako argumenty...
  331.  
  332. pr = -kmax : kmax;               % Oś przesunięcia (w próbkach).
  333. tr = -kmax/fp : 1/fp : kmax/fp;  % Oś opóźnienia/przesunięcia (w sekundach).
  334. %    ^^^^^^^^   ^^^^   ^^^^^^^
  335. %       |        |        \ Maksymalne opóźnienie sygnału.
  336. %       |         \ Czas próbkowania, wielkość co którą będą pojawiać się próbki czasu.
  337. %        \ Ujemne opóźnienie sygnału.
  338.  
  339. figure(3);
  340. subplot(413);                          % Wybór komórki dla trzeciego wykresu.
  341. plot(tr, r);                           % Wykres autokorelacji.
  342. xlabel('Przesunicie [s]');             % Opis poziomej osi.
  343. title('Korelacja sylab "wy" i "czy"'); % Tytuł.
  344.  
  345. rwy = xcorr(wy,wy,kmax);               % Autokorelacja sylaby 'wy'.
  346. rczy = xcorr(czy,czy,kmax);            % Autokorelacja sylaby 'czy'.
  347.  
  348. subplot(427);                          % Wybór komórki.
  349. plot(tr, rwy);                         % Wykres autokorelacji.
  350. xlabel('Przesunięcie [s]');            % Opis poziomej osi.
  351. title('Autokorelacja sylaby "wy"');    % Tytuł.
  352.  
  353. subplot(428);                          % Wybór komórki.
  354. plot(tr, rczy);                        % Wykres autokorelacji.
  355. xlabel('Przesunięcie [s]');            % Opis poziomej osi.
  356. title('Autokorelacja sylaby "czy"');   % Tytuł.
  357.  
  358.  
  359. %% 11) Definiowanie filtra i filtracja
  360. % Jest wiele rodzajów filtrów cyfrowych, tutaj omówione są filtry o skończonej odpowiedzi impulsowej.
  361. % Filtr FIR (finite impulse response) wymaga:
  362. % - Zdefiniowania długości odpowiedzi impulsowej.
  363. % - Określenia zakresu filtracji w częstotliwości unormowanej.
  364. %   (Częstotliwość unormowana jest zdefiniowana jako procent częstotliwości Nyquista.)
  365. % - Określenia typu filtracji.
  366.  
  367. M = 401; % Długość odpowiedzi filtra. Powinna być to liczba nieparzysta.
  368.  
  369. fg = 3000;                   % Częstotliwość graniczna, wyrażona w hercach.
  370. fgn = fg / (fpx/2);          % Częstotliwość graniczna unormowana.
  371. %     ^^    ^^^^^
  372. %     |       \ ... podzielona przez połowę częstotliwości próbkowania (częstotliwość Nyquista).
  373. %      \ Częstotliwość graniczna...
  374. hl = fir1(M-1, fgn, 'low');                     % Odpowiedź impulsowa filtra dolnoprzepustowego
  375. %    ^^^^ ^^^  ^^^  ^^^^^^
  376. %     |    |    |      \ Typ filtra. Są cztery typy, pozostałe zdefiniowano poniżej.
  377. %     |    |     \ Częstotliwość graniczna
  378. %     |     \ Długość odpowiedzi impulsowej. Dla zachowania spójności musi być pomniejszona o 1.
  379. %      \ Funkcja generująca filtr.
  380. hh = fir1(M-1, fgn, 'high');                    % Filtr górnoprzepustowy.
  381.  
  382. % Możliwe jest definiowanie filtrów, które modyfikują zakres widma.
  383. fg_min = 5000;                                  % Dolny zakres częstotliwości granicznej.
  384. fg_max = 7500;                                  % Górny zakres częstotliwości granicznej.
  385. fgn_min = fg_min / (fpx/2);                     % Unormowana częstotliwość graniczna
  386. fgn_max = fg_max / (fpx/2);                     % Jak wyżej.
  387. hb = fir1(M-1, [fgn_min, fgn_max], 'bandpass'); % Filtr pasmowoprzepustowy.
  388. hs = fir1(M-1, [fgn_min, fgn_max], 'stop');     % Filtr pasmowozaporowy.
  389.  
  390. ph = 0:(M-1);                                   % Oś próbek dla wykresu odpowiedzi impulsowej.
  391. th = 0 : 1/fp : (M-1)/fp;                       % Os czasu.
  392.  
  393. % Aby zobaczyć działanie filtra odkomentuj wybraną linijkę.
  394. h = hl;                         % 'low', filtr dolnoprzepustowy
  395. % h = hh;                         % 'high', filtr górnoprzepustowy
  396. % h = hb;                         % 'bandpass', filtr pasmowoprzepustowy
  397. % h = hs;                         % 'stop', filtr pasmowozaporowy
  398.  
  399. y = filter(h,1,x);                % Filtracja. Sygnał przefiltrowany jest zapisany w zmiennej y.
  400. %   ^^^^^^ ^ ^ ^
  401. %     |    | | \ Filtrowany sygnał.
  402. %     |    |  \ Jedynka. Wstawimy tu inną wartość, gdy będziemy korzystać z filtrów IIR (nieskończona odpowiedź impulsowa).
  403. %     |     \ Odpowiedź impulsowa.
  404. %      \ Funkcja realizująca filtrację.
  405.  
  406. vy = fft (y,Nfx);                 % Widmo zespolone sygnału po filtracji.
  407. avy = abs(vy);                    % Widmo amplitudowe sygnału po filtracji.
  408.  
  409. % Filtry nie są sygnałami, tylko systemami, jednak możemy je analizować w dziedzinie częstotliwości.
  410. % "Widmo" odpowiedzi impulsowej nazywamy transmitancją (transfer function, funkcja transmisji).
  411. vh = fft (h,Nfx);                 % Transmitancja.
  412. avh = abs(vh);                    % Moduł transmitancji.
  413.  
  414. figure(4);
  415. subplot(321);
  416. plot(tx,x);
  417. xlabel('Czas [s]');
  418. title('Sygnał przed filtracją');
  419.  
  420. subplot(322);
  421. plot(fx,avx(1:Nfx21));
  422. xlabel('Częstotliwość [Hz]');
  423. title('Moduł widma sygnału');
  424.  
  425. subplot(323);
  426. plot(th, h);
  427. xlabel('Czas [s]');
  428. title('Odpowiedź impulsowa');
  429.  
  430. subplot(324);
  431. plot(fx,avh(1:Nfx21));
  432. xlabel('Częstotliwość [Hz]');
  433. title('Modul transmitancji');
  434.  
  435. subplot(325);
  436. plot(tx, y);
  437. xlabel('Czas [s]');
  438. title('Sygnał po filtracji');
  439.  
  440. subplot(326);
  441. plot(fx,avy(1:Nfx21));
  442. xlabel('Częstotliwość [Hz]');
  443. title('Moduł widma sygnału po filtracji');
RAW Paste Data