1 MATLAB

2 Literatura B. Mrozek, Z. Mrozek, „Matlab 5.x, Simulink 2.x – poradnik użytkownika”, Wyd. PLJ Warszawa. J. Brzózka, L. Dorobczyński, „Programowanie w Matlabie”, Mikom, Warszawa 1998. A. Zalewski, R. Cegieła, „Obliczenia numeryczne i ich zastosowania”, Nakom 1997. D. Hanselman, B. Littliefield, „Mastering Matlab 6, A Comprehensive Tutorial and Reference, Prentice Hall, 2001. D. J. Higham, N. J. Higham, „Matlab Guide, SIAM, Philadelphia.

3

4 Scilab, Octave, i in. – darmowe wersje MATLAB’a

5 Jak powstawał Matlab Lata 70 – Univ. New Mexico (USA), Cleve Moler, MATrix LABoratory (za pomocą Fortranu), wspomaganie zajęć z algebry 1984 – nowa edycja (na bazie C) 1987 – Math Works Inc., Matlab 3.0, ulepszony interpreter, grafika 1992 – Matlab 4.0, system pod Windows, animacje, GUI, macierze rzadkie 1997 – Matlab 5.0, programowanie obiektowe, macierze wielowymiarowe, hipertekstowy HELP, nowe narzędzia ODE 2000 – Matlab 6.0, współpraca z Javą, PDE

6 CZYM JEST MATLAB? Matlab to pakiet przeznaczony do wykonywania obliczeń numerycznych oraz graficznej prezentacji wyników. Interpreter skryptowy. Dostępny jest na różnych platformach sprzętowych oraz systemowych (np.. Windows, Macintosh). Podstawową strukturą danych w Matlabie jest macierz.

7 Kiedy Matlab ? Interaktywny język wysokiego poziomu. Przejrzysty kod.System do obliczeń numerycznych operuje na tablicach (macierzach) danych Optymalizacja czasochłonnych operacji tablicowych Język skryptowy – interpretowany (możliwość kompilacji poprzez jęz. C). Bogate (darmowe) biblioteki – otwarty kod, dowolny system operacyjny). Możliwość łączenia z funkcjami w C/C++, klasy Javy

8 Środowisko pracy – okno poleceń

9 Simulink

10 PRACA Z PAKIETEM MATLABW trybie bezpośrednim – typowy tryb roboczy, umożliwiający prowadzenie dialogu pomiędzy użytkownikiem a pakietem na zasadzie: pytanie-odpowiedź. Użytkownik wpisuje polecenia bezpośrednio do okna poleceń W trybie pośrednim – umożliwiającym szybkie i efektywne wykonanie obliczeń i prezentację wyników za pomocą uruchomienia programu napisanego w języku pakietu Matlab, czyli tzw. Skryptu (zwanego również m-plikiem, np. moj_skrypt.m). >> moj_skrypt

11 POLECENIA Po wydaniu polecenia i naciśnięciu klawisza Enter Matlab natychmiast wyświetla jego wynik. Umieszczenie po poleceniu średnika spowoduje wykonanie obliczeń, ale bez zwracania wyniku. Polecenie powinno się mieścić w jednym wierszu. Kilka poleceń w jednym wierszu oddzielamy od siebie przecinkami lub średnikami.

12 Matlab (interpreter) jako kalkulator

13 Charakterystyka M-plikówSkrypty Nie posiadają argumentów wejściowych, ani wyjściowych Operują na danych w przestrzeni roboczej Wygodne, kiedy często należy wykonać sekwencję tych samych poleceń Funkcje Mogą posiadać argumenty wejściowe i wyjściowe Wewnętrzne zmienne są lokalne dla funkcji Wygodne do rozszerzenia MATLABA o własne aplikacje

14 M-pliki

15 Edytor Matlaba

16 Matlab - język numeryczny

17 POMOC SYSTEMOWA Uzyskanie informacji o funkcjach Matlaba:>>help nazwa_funkcji Help Desk-podręcznik opracowany w postaci stron HTML.

18 Help

19 LICZBY Stałopozycyjne-z opcjonalnym użyciem znaku + lub – oraz kropki dziesiętnej; Zmiennopozycyjnej - z użyciem znaku e lub E poprzedzającego wykładnik potęgi 10, np. 1e2=100; Do zapisu części urojonej liczb zespolonych używa się stałej i lub j. UWAGA: Domyślnie Matlab traktuje wszystkie liczby jako zespolone (ostrożnie z pierwiastkowaniem)

20 ZMIENNE Nazwa zmiennej musi rozpoczynać się literą i może składać się z dowolnej liczby liter, cyfr i znaków podkreślenia. Pakiet Matlab nie wymaga deklarowania zmiennych ani określenia ich rozmiaru (można, w uzasadnionych wypadkach - b. duże macierze). Aby sprawdzić wartość istniejącej już zmiennej, należy w wierszu poleceń wpisać jej nazwę. Np. >>A Matlab rozróżnia duże i małe litery. Standardowe polecenia pakietu pisane są zawsze małymi literami.

21 Deklaracja zmiennych

22 DEFINIOWANIE MACIERZY (WEKTORA)Elementy w wierszu macierzy muszą być oddzielane spacją lub przecinkami; A=[1 3 4;3 4 5]; B=[1,2;7,8]; C=[3:7]; D=[2:0.1:15]; Średnik lub znak nowego wiersza kończy wiersz macierzy i powoduje przejście do następnego; Cała lista elementów musi być ujęta w nawiasy kwadratowe.

23

24 ODWOŁANIA DO FRAGMENTÓW MACIERZYx(j:k) – elementy wektora wierszowego x o numerach od j do k A(i,:) – wszystkie elementy w wierszu i macierzy A A(i,j:l)- wszystkie elementy w wierszu i macierzy A o numerach od j do l A(i:k,j:l)-wszystkie elementy w kolumnach od j do l wierszy od i do l

25 A(x,j:l)-wszystkie elementy w kolumnach od j do l w wierszach macierzy A o numerach określonych przez elementy wektora x A(:,:) – cała dwuwymiarowa macierz A A(:)-cała macierz A w postaci wektora kolumnowego.

26 Przeszukiwanie macierzyfind (A>3) A>3

27

28  A=[17,0,-9;30,-2,38]  A(A>1)  ans = 17 30 38

29 WYŚWIETLANIE MACIERZY I ICH ROZMIARÓWdisp(A)-wyświetla zawartość macierzy A w oknie poleceń; size(A)-wyświetla rozmiar dwuwymiarowej macierzy A (liczbę wierszy i kolumn) w postaci dwuelementowego wektora wierszowego; [n m]=size(A)-przypisuje zmiennej n liczbę wierszy, a zmiennej m liczbę kolumn;

30 n=size(A,1)-przypisuje zmiennej n liczbę wierszy macierzy A;m=size(A,2)-przypisuje zmiennej m liczbę kolumn macierzy A; length(x)-zwraca długość wektora x lub dłuższy z wymiarów macierzy.

31 ARYTMETYKA MACIERZOWA I TABLICOWAA*B, B*A C1=B/A, C2=A\B A^2=A*A A’ A+B A-B A.*B=B.*A B./A=A.\B A.^2 A’

32 Wektory i macierze - podstawowy typ zmiennych Matlaba

33 FUNKCJE GENERUJĄCE I PRZEKSZTAŁCAJĄCE MACIERZEeye(n)-tworzy macierz jednostkową nxn; ones(n)-tworzy macierz nxn o elementach równych 1; zeros(n)-macierz zerowa nxn; rand(n)-macierz nxn wypełniona liczbami pseudolosowymi z przedziału <0,1> o rozkładzie jednostajnym; randn(n)-macierz nxn wypełniona liczbami pseudolosowymi o rozkładzie normalnym ze średnią 0 i wariancją 1.

34 A=diag(x)-macierz przekątniowa A ze składnikami wektora x na głównej przekątnej;x=diag(A)-utworzenie wektora x z elementów znajdujących się na głównej przekątnej macierzy A; inv(A)-utworzenie macierzy odwrotnej do A; repmat(A,n,m)-utworzenie macierzy przez powielenie podmacierzy A m razy w poziomie i n razy w pionie;

35 reshape(A,n,m)-utworzenie macierzy o n wierszach i m kolumnach z elementów branych kolejno kolumnami z macierzy A; rot90(A)-obrócenie macierzy A o 90 stopni w kierunku przeciwnym do wskazówek zegara; tril(A)-utworzenie z macierzy A macierzy trójkątnej dolnej; triu(A)- utworzenie z macierzy A macierzy trójkątnej górnej.

36 Zmiana kształtu macierzy size(A) ans = length(A) 2 Zmiana kształtu macierzy  C=reshape(A,2,3)  C =

37 Funkcje logiczne służąc do badania własności całych macierzyexist('nazwa') isempty(X) issparce(X) isstr(X) isglobal(X)

38 Funkcje logiczne badające własności elementów macierzyany(X) – czy którykolwiek element niezerowy? all(X) – czy wszystkie niezerowe isinf(X) finite(X) I=find(X)

39 MACIERZE WIELOWYMIAROWEMatlab dopuszcza definiowanie macierzy wielowymiarowych. Odwoływanie się do elementów takich macierzy wymaga liczby indeksów większej niż 2. Pierwszy indeks-wiersz macierzy (wymiar 1); Drugi indeks-kolumna macierzy (wymiar 2); Trzeci indeks-strona macierzy (wymiar 3) Czwarty indeks-książka macierzy (wym. 4); Piąty indeks-tom macierzy (wym.5); itd.;

40 METODY TWORZENIA TABLIC WIELOWYMIAROWYCHprzez indeksowanie; przez zastosowanie funkcji (ones, zeros, randn, repmat-tworzy tablice wielowymiarową wypełnioną jednakowymi wartościami); przez zastosowanie funkcji cat (konkatenacja, scalanie tablic); cat(dim,A,B) scala dwie macierze A i B zgodnie z podanym wymiarem dim; cat(2,A,B) oznacza to samo co [A, B]; cat(1,A,B) oznacza to samo co [A; B]; B=cat(dim,A,B,C,...) scala macierze A, B,C,... zgodnie z podanym wymiarem dim;

41 PRZYKŁAD Na rysunku widoczna jest macierz trójwymiarowa o rozmiarze 2x3x2 (2 wiersze i 3 kolumny na każdej stronie, 2 strony); >>D(:,:,1)=[1 3 0; 5 7 2] % str.1; >>D(:,:,2)=[4 7 8; 1 0 5] % str.2; kolumny strony wiersze

42 TYPY DANYCH Matlab dopuszcza użycie sześciu podstawowych typów danych:Double-liczby podwójnej precyzji; podstawowy typ danych dla zmiennych MATLAB-a (wszystkie obliczenia w Matlabie są prowadzone w trybie podwójnej precyzji dla zmiennych numerycznych i łańcuchowych); DOMYŚLNY 2. Char-znaki i łańcuchy znaków; łańcuch znakowy definiuje się za pomocą apostrofów i przechowywany jest w pamięci w postaci wektora liczb całkowitych reprezentujących kody ASCII poszczególnych znaków;

43 3. Sparse - dotyczy dwuwymiarowych macierzy rzadkich podwójnej precyzji; (macierz rzadka to taka macierz, w której zapamiętywane są tylko elementy niezerowe; redukuje to zapotrzebowanie pamięci); 4. Cell-typ komórkowy; elementy tablic komórkowych mogą zawierać inne tablice; 5. Struct-typ strukturalny; tablice strukturalne odwołują się do nazw pól, które mogą zawierać inne tablice; 6. Uint8-typ przeznaczony do efektywnego wykorzystania pamięci, integer ze znakiem lub bez (uint); możliwe są takie operacje, jak zmiana wymiarów lub kształtu tablicy, ale niedozwolone są żadne operacje matematyczne; Funkcje zamiany: i = int8(x) i = int16(x) i = int32(x) i = int64(x) uint8(x) Oprócz wymienionych typów istnieje typ UserObject, który jest typem definiowanym przez użytkownika.

44 FUNKCJE MATLABA Wbudowane-część jądra pakietu, do których użytkownicy nie mają dostępu (np.sqrt); Implementowane w m-plikach -przechowywane w ogólnie dostępnych plikach np. peaks, takie m-pliki użytkownicy mogą tworzyć sami;

45 LISTA ŚCIEŻEK Lista ścieżek to lista katalogów, do których Matlab ma dostęp. Polecenie „path” Jest zdefiniowana w pliku pathdef.m, znajdującym się w podkatalogu toolbox\local katalogu z Matlabem. Można ją wyświetlać lub zmieniać poleceniem path.

46 PODSTAWOWE FUNKCJE I STAŁE MATEMATYCZNEFunkcja opis sin(z), cos(z), tan(z), cot(z) Sinus, cosinus, tangens, cotangens; argument funkcji w radianach; asin(z), acos(z), atan(z), acot(z) Funkcje cyklometryczne; wynik w radianach; sinh(z), cosh(z), tanh(z), coth(z) Funkcje hiperboliczne; argument w radianach; sqrt(z) Pierwiastek z ; z<0 – wynik zespolony;

47 exp(z) log(z) log2(z) abs(z) angle(z) real(z), imag(z) conj(z) ezln(z); z<0 – wynik zespolony; log2(z) log2z; z<0 – j.w. abs(z) wartość bezwzględna lub moduł liczby zespolonej; angle(z) argument liczby zespolonej; real(z), imag(z) część rzeczywista i urojona liczby z conj(z) liczba zespolona sprzężona;

48 complex(x,y) ceil(z) floor(z) fix(z) round(z) rem(x,y); mod(x,y)utworzenie liczby zespolonej; ceil(z) zaokrąglenie liczby w górę; floor(z) zaokrąglenie liczby w dół; fix(z) zaokrąglenie liczby dodatniej w dół, ujemnej w górę; round(z) zaokrąglenie do najbliższej liczby całkowitej; rem(x,y); mod(x,y) reszta z dzielenia x przez y; sign(x) funkcja signum;

49 FUNKCE OPERUJĄCE NA WEKTORACH (na MACIERZACH podane wcześniej)max(x) największy element wektora x; min(x) najmniejszy element wektora x; sum(x) sumę elementów wektora x; prod(x) iloczyn elementów wektora x; mean(x) średnia arytmetyczna elementów wektora x; length(x) długość wektora

50 pierwiastek z liczby –1; Inf lub inf STAŁE MATEMATYCZNE stałe opis pi przybliżenie wartości eps względna dokładność zmiennoprzecinkowa; i lub j pierwiastek z liczby –1; Inf lub inf nieskończoność (ang. Infinity); jest rezultatem operacji, która przekracza zakres arytmetyki komputera, np.dzielenie przez 0; NaN lub nan nie liczba; jest wynikiem matematycznie niezdefiniowanych operacji;

51 Instrukcje sterujące M-pliki: skrypty i funkcje. Przykłady skryptów. Globalność zmiennych. Metody debuggowania. Funkcje Struktura funkcji Zmienne nargin, nargout, zmienne lokalne, zmienne globalne Subfunkcje i funkcje prywatne. Pseudokompilacje funkcji, usuwanie funkcji z pamięci Przykłady funkcji Kolejne funkcje własne Matlab’a

52 Przykład 1- life C:\Matlab\toolbox\matlab\demos

53 PODSTAWOWE FUNKCJE I STAŁE MATEMATYCZNEFunkcja opis sin(z), cos(z), tan(z), cot(z) Sinus, cosinus, tangens, cotangens; argument funkcji w radianach; asin(z), acos(z), atan(z), acot(z) Funkcje cyklometryczne; wynik w radianach; sinh(z), cosh(z), tanh(z), coth(z) Funkcje hiperboliczne; argument w radianach; sqrt(z) Pierwiastek z ; z<0 – wynik zespolony;

54 exp(z) log(z) log2(z) abs(z) angle(z) real(z), imag(z) conj(z) ezln(z); z<0 – wynik zespolony; log2(z) log2z; z<0 – j.w. abs(z) wartość bezwzględna lub moduł liczby zespolonej; angle(z) argument liczby zespolonej; real(z), imag(z) część rzeczywista i urojona liczby z conj(z) liczba zespolona sprzężona;

55 complex(x,y) ceil(z) floor(z) fix(z) round(z) rem(x,y); mod(x,y)utworzenie liczby zespolonej; ceil(z) zaokrąglenie liczby w górę; floor(z) zaokrąglenie liczby w dół; fix(z) zaokrąglenie liczby dodatniej w dół, ujemnej w górę; round(z) zaokrąglenie do najbliższej liczby całkowitej; rem(x,y); mod(x,y) reszta z dzielenia x przez y; sign(x) funkcja signum;

56 FUNKCE OPERUJĄCE NA WEKTORACH (na MACIERZACH podane wcześniej)max(x) największy element wektora x; min(x) najmniejszy element wektora x; sum(x) sumę elementów wektora x; prod(x) iloczyn elementów wektora x; mean(x) średnia arytmetyczna elementów wektora x; length(x) długość wektora

57 1.     Kolejne funkcje własne Matlab’atype – type mój %wyświetla zawartość skryptu moj disp – wyświetla zmienną bez jej nazwy, num2str – numeryczna -> string disp(z), disp(‘obliczenia zakończone a=’); disp (a), greckie litery nie działają. pause – zatrzymuje wykonywanie m-pliku do wciśnięcia klawisza, pause(n)- na n sekund, pause on, pause off, umożliwia lub uniemożliwia ponowne wykorzystanie tej funkcji input – podanie zmiennej input(‘podaj a, a=’) lub a= input(‘podaj a, a=’) error – komunikat błędu wewnątrz funkcji, przerwanie f-cji i przekazanie sterowania do okna poleceń ... error(‘za mało argumentów’)

58 menu – zwraca numer wyboru:K = menu('Choose a color','Red','Blue','Green') (uwaga: istnieje również funkcja uimenu w GUI) keyboard – przerywa wykonywanie funkcji i wraca do okna poleceń „>>K”. Po wpisaniu „return” wraca do funkcji. plot - wykresy: plot (x,y) - najprostsze zastosowanie, gdzie np. x=[…], y=funkcja([x]) clock – zwraca 6-elementowy wektor zawierający obecną datę i czas w formie dziesiętnej c=clock => c=[rok, miesiąc, dzień, godzina, minuty, sekundy] find – znajduje indeksy elementów spełniających podany warunek X=[ find(X~=0); Przykład: Usuwanie elementów zerowych X=[ ]; index=find(X~=0); nowy=X(index) =>[ ];

59 size – rozmiar macierzyall – czy wszystkie elementy wektora są niezerowe lub w kolumnie macierzy, isnan – czy liczba jest NaN => 0 lub 1 (0/0=>NaN) isinf– czy liczba jest Inf => 0 lub 1 finite– czy liczba jest skończona rem – reszta z dzielenia rem(6,4)=>2

60 roots – oblicza miejsca zerowe wielomianu.Algorytm w oparciu o obliczenie wartości własnych macierzy współczynników (patrz help) przykład: s3-6s2-72s-27 p=[1 –6 –72 –27]; r=roots(p) => ; –5.7345; –0.3884 prod – mnożenie elementów wektora lub elementów kolumn macierzy cross – cross(A x B) – iloczyn wektorowy dot – iloczyn skalarny eval , feval (przy łańcuchach znaków) !!!!! >> R=input(‘co robic’,'s'); eval(R) >> F feval(F,9.64)  foo(9.64)

61 try wyrażenie catch zrób_gdy_wyrażenie_niegramatyczne endINSTRUKCJE POMOCNICZE error(‘tekst’) – opuszczenie programu lasterr – zwraca ostatni komunikat o błędzie, często używana z funkcją eval try wyrażenie catch zrób_gdy_wyrażenie_niegramatyczne end - Obsługa błędu składni return- opuszczenie danego skryptu lub funkcji, a następnie powrót do wywołującego ją programu warning (‘tekst’)- wyświetla tekst podobnie jak disp, ale może tego uniknąć poprzez polecenie Warning off (Warning on). Warning backtrace – wyświetla tekst podając dodatkowo ścieżkę i nazwę pliku, skąd jest „warning”.

62 Instrukcje sterujące IF if n>0 …….. (instrukcje) elseif n = = 0IF if n>0 …….. (instrukcje) elseif n = = 0 else end

63 Relacje i wyrażenia logiczne Operatory porównania (przy instrukcji if ..): A= =B A~ =B A A<=B A>B A>=B

64 WHILE while prod(1:n)<10^100 n=n+1 end FOR for i=1:n % i=[1:n] x(i)=0 !!! przykład nieefektywnego programowania  zeros(1:n)

65 SWITCH: SWITCH zmienna CASE wartość zmiennej, polecenia CASE {wartość1, wartość2,…} ... OTHERWISE, END Uwaga: W przeciwieństwie do języka C Matlab wykonuje tylko pierwszy pasujący przypadek i nie sprawdza pozostałych. Dlatego nie jest używana instrukcja „break”

66 Przykład: w=4; s=0.5; operator=’*’; switch operator case ‘+’ w=w+s; case ‘{‘*’,’.*’} w=w*s; case ‘/’ w=w/s otherwise w=w-s; end

67    Skrypty Wywołanie skryptu poprzez nazwę. Zmienne skryptu są przechowywane w pamięci jako zmienne globalne. Debuggowanie – menu edytora Debug lub Breakpoints. Ewentualne wyświetlanie wartości zmiennych (brak średnika).

68 Funkcja Ciało funkcji: function b=oblicz(x)% pomoc (help) do uzyskania po poleceniu „help oblicz” % j.w (wolna linia) % komentarz programisty % chcemy uzyskać funkcję f(x)=sin(x2)/(x2-x) b=sin(x.^2)./(x.^2-x); % średnik zapobiega podwójnemu wyświetleniu wyniku wywołanie: >> oblicz(3) >> u=oblicz([1,3,8,7])

69 Uwagi: ·        a) Bardzo ważna jest nazwa pliku, a nie nazwa f-cji zadeklarowana wewnątrz pliku (nazwa funkcji do 31 znaków). ·       b) Funkcja lookfor rozpoznaje słowa tylko z I linii. ·      c)  Wszystkie zmienne w funkcji są lokalne (zapominane po wyjściu – sprawdzić na ćw.). ·       d) Nie można wewnątrz funkcji odwoływać się do zmiennych zewnętrznych. Z (c) i (d) wynika, że można dublować nazwy zmiennych e) Funkcję można wywoływać rekurencyjnie, ale patrz (c).

70 · Można zmienną zmienić na globalną·        Można zmienną zmienić na globalną. Należy wewnątrz funkcji zadeklarować global nazwa_zmiennej. Wówczas zmienna jest pamiętana w kolejnych wywołaniach funkcji, ale nie jest pamiętana w przestrzeni roboczej i innych funkcjach lub skryptach. Można to uzyskać deklarując global nazwa_zmiennej również tam. Wówczas jest to sama zmienna. ·        clear global – usuwa wszystkie zmienne globalne. ·        Funkcja może mieć kilka argumentów wejściowych (input arguments) ·        Ilość argumentów wejściowych może być zmienna ·        Funkcja może zwracać kilka argumentów wyjściowych (output arguments) przykład: [X,Y]=fun2(3), gdzie X i Y mogą być macierzami o innych wymiarach.

71  Funkcja ze zmienną globalną:function nasza1(x) global a; b=a+x;  function nasza2(x) b=a-x; >> global a; a=7; >>nasza1(4)*nasza2(3) =>11*4=44 Można teraz zmieniać zmienną a bezpośrednio z przestrzeni roboczej bez konieczności edycji obu funkcji.

72 Funkcja zwracająca 2 wartości:function [mean,stdev]=stat(x) [m,n]=size(x) if m= =1 m=n end mean=sum(x)/m stdev=sqrt(sum(x.^2)/(m-mean) wywołanie: >> [M,S]=stat(rand(10,1000))

73 Funkcje o zmiennej ilości argumentów wej/wyjściowychnargin (number of input arguments) nargout (number of output arguments)  Przykład: function y=wiele(x1,x2) if nargin= = 1 x2=10; end y=sin(x1*x2);  >> wiele(pi,2) =>sin(2 pi)=0 >> wiele(pi/20) =>sin(pi/2)=1  nargout dotyczy sposobu wywołania funkcji w innej funkcji lub oknie poleceń.

74 W przypadku kiedy w deklaracji funkcji nie chcemy określić maksymalnej liczby argumentów wejściowych lub wyjściowych można użyć funkcji: varargin varargaout Przykład: function wynik=zmienna2(varargin) a=1; z=varargin{2}; wynik=a+z; >> zmienna2(2,6,11) =>7

75 Przestrzeń robocza funkcji:Każda funkcja ma przydzieloną pamięć odseparowana od przestrzeni roboczej Matlab’a. Dlatego funkcje operuja na zmiennych lokalnych. Po wywołaniu funkcji Matlab przekształca funkcję w pseudokod i zapamiętuje go, co przyspiesz wykonywanie funkcji. Pseudokod można usunąć z pamięci komputera poleceniem: clear nazwa_funkcji clear functions – usuwa wszystkie funkcje clear all – usuwa wszystkie funkcje i zmienne

76 Pseudokompilacja funkcji Do przyspieszenia działania dużych programów z rozbudowanym GUI (lub ukrycia kodu). Uruchomienie funkcji z Matlab’a t.j m-plik, ale nie można po kompilacji zmieniać nazwy pliku zewnętrznie. pcode nazwa_funkcji => p-plik Powstaje plik nazwa_funkcji.p

77 Subfunkcje ·        M-pliki mogą zawierać kody więcej niż jednej funkcji. Funkcja zapisana jako pierwsza jest funkcją główna i jej nazwa pokrywa się z nazwą m-pliku. ·        Zmienne wszystkich funkcji z tego samego m-pliku są nadal lokalne dla tych funkcji. ·        Subfunkcje mogą być wywołwane tylko przez funkcje z tego samego m-pliku.

78 Przykład (Higham): function max_err=poly1err(n)% POLY1ERR(N) Błąd w interpolacji liniowej dowolnej funkcji (p=a*x+c; a=(f1-f0)/(1-0); c=f0) max_err=0; f0=f(0); f1=f(1); for x=linspace(0,1,n); p=x*f1+(1-x)*f0; err=abs(f(x)-p); max_err=max(max_err,err); end; % Subfunkcja function y=f(x) % Function to be interpolated, F(x) y=sin(x);

79 Ta sama funkcja z optymalnym wykorzystaniem operacji wektorowych: function max_err=poly2err(n) f0=f(0); f1=f(1); x=linspace(0,1,n); p=x*f1-(x-1)*f0; err=abs(f(x)-p); max_err=max(err); % Subfunkcja function y=f(x) % Function to be interpolated, F(x) y=sin(x);

80 Przykład – funkcja PEAKS

81 Przykład – funkcja PEAKSfunction [xz,y,z] = peaks(arg1,arg2); %PEAKS A sample function of two variables. % PEAKS is a function of two variables, obtained by translating and % scaling Gaussian distributions, which is useful for demonstrating % MESH, SURF, PCOLOR, CONTOUR, etc. % There are several variants of the calling sequence: % % Z = PEAKS; Z = PEAKS(N); Z = PEAKS(V); % Z = PEAKS(X,Y); PEAKS; PEAKS(N); % PEAKS(V); PEAKS(X,Y); [X,Y,Z] = PEAKS; % [X,Y,Z] = PEAKS(N); [X,Y,Z] = PEAKS(V);

82 function [xz,y,z] = peaks(arg1,arg2);if nargin = = 0 [x,y] = meshgrid(-3:dx:3); elseif nargin = = 1 if length(arg1) = = 1 [x,y] = meshgrid(-3:6/(arg1-1):3); % grafika3D else [x,y] = meshgrid(arg1,arg1); % grafika3D end x = arg1; y = arg2; z = 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) - 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ... - 1/3*exp(-(x+1).^2 - y.^2);

83 if nargout > 1 xz = x; elseif nargout = = 1 xz = z; else % Self demonstration disp('z = 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) ... ') disp(' *(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ... ') disp(' - 1/3*exp(-(x+1).^2 - y.^2) ') surf(x,y,z) % grafika3D axis([min(min(x)) max(max(x)) min(min(y)) max(max(y)) ... min(min(z)) max(max(z))]) xlabel('x'), ylabel('y'), title('Peaks') end

84 ŚRODOWISKO MATLABA Okno poleceń. Funkcje obsługujące okno poleceń. clc lit Powtórzenie ostatniego polecenia lub ostatniego polecenia zaczynającego się od „lit” clc wyczyszczenie okna poleceń i umieszczenie kursora w jego lewym górnym rogu; home umieszczenie wiersza poleceń i kursora w lewym górnym rogu okna poleceń; Ctrl+C przerwanie obliczeń

85 echo on/echo off more on/more off diary plik diary off/onwłącza/wyłącza wysyłanie na ekran treści wykonywanych poleceń; more on/more off włącza/wyłącza stronicowanie tekstów wysyłanych na ekran; diary plik polecenia i teksty (bez grafiki) wysyłane na ekran będą zapisywane w pliku o podanej nazwie; diary off/on przełącznik funkcji diary loose/compact zmiana interlinii w wyświetlanym tekście

86 Format Opis Wynik dla 1/23 short 0,0435 short e 4,3478e-002 longFormaty liczb. Do określenia sposobu wyświetlania liczb rzeczywistych w oknie służy funkcja format. Użycie funkcji nie ma wpływu na dokładność wykonywanych obliczeń, a tylko na widok liczby na ekranie. UWAGA: Przed każdym poleceniem dopisać format Format Opis Wynik dla 1/23 short 5-cyfr.liczba stałopozycyjna (format domyślny); 0,0435 short e 5-cyfr.liczba zmiennopozycyjna; 4,3478e-002 long 15-cyfr.liczba stałopozycyjna; 0,

87 15-cyfr.liczba zmiennopozycyjna; short g long e 15-cyfr.liczba zmiennopozycyjna; 4, e-002 short g 5 znaczących cyfr liczby stało- lub zmiennopozycyjnej; 0,043478 long g 15 znaczących cyfr liczby stało- lub zmiennopozycyjnej; 0, hex liczba szestnastkowa; 3fa642c8590b2164 + drukuje znak + dla liczb dodatnich,- dla ujemnych, spację dla zera; bank format walutowy 0,04 rat przybliża liczbę ułamkami małych liczb całkowitych; 1/23

88 pierwiastek z liczby –1; Inf lub inf STAŁE MATEMATYCZNE stałe opis pi przybliżenie wartości eps względna dokładność zmiennoprzecinkowa; i lub j pierwiastek z liczby –1; Inf lub inf nieskończoność (ang. Infinity); jest rezultatem operacji, która przekracza zakres arytmetyki komputera, np.dzielenie przez 0; NaN lub nan nie liczba; jest wynikiem matematycznie niezdefiniowanych operacji;

89 Reprezentacja zmiennych i format wyświetlania zmienne typu "double precision", domyślnie wyświetla 4 miejsca po przecinku.  Zmiana formatu wyświetlania – "format" a=13/3  a = 4.3333  format long a

90 format short – 4 miejsc (stało-przecinkowa) – domyślnyformat – powrót do domyślnego long – 15 (stało) short e – 5 (zmienno-przecinkowa) long e – 15 (zmienno) rat – w postaci ułamka format rat a  a = 13/3

91 format a  a = 4.3333 compact – pomijanie pustych linii przy wyświetlaniu loose – wprowadzenie pustych linii przy wyświetlaniu  UWAGA NA NAZWY: Zmienne specjalne i stałe  realmax  realmin – najmniejsza liczba, która można dodać do 1 I otrzymać większa od 1)  realmin

92 0/0  NaN  9/0  Inf ans beep pi eps  i lub j

93 nargin (ile we argumentów funkcji)nargout (ile wy wartości funkcji) Największa użytkowa dodatnia liczba całkowita bitmax varargin – zmienna ilość arg. we funkcji varargout – zmienna ilość wartości wy funkcji

94 Przykład 1- life C:\Matlab\toolbox\matlab\demos

95

96 n = [m 1:m-1]; e = [2:m 1]; s = [2:m 1]; w = [m 1:m-1]; while get(axHndl,'UserData')==play, % How many of eight neighbors are alive. N = X(n,:) + X(s,:) + X(:,e) + X(:,w) + ... X(n,e) + X(n,w) + X(s,e) + X(s,w);

97 Wykresy i grafika w MatlabieOperacje na plikach Wykresy 2D Wykresy 3D Animacje i filmy Tworzenie interfejsu użytkownika GUI Ręczne Wspomagane (GUIDE)

98 (standardowa obsługa zmiennych)PRACA NA PLIKACH (I) (standardowa obsługa zmiennych) save nazwa_pliku nazwa_zmiennej save nazwa_pliku nazwa_zmiennej –ascii load nazwa_pliku load nazwa_pliku zmienne (jeżeli format .mat)

99 zapisuje binarnie wszystkie zmienne w pliku matlab.mat; save plik zapisuje binarnie wszystkie zmienne w pliku o nazwie plik.mat; save plik lista zapisuje binarnie w pliku o nazwie plik.mat tylko zmienne wymienione; load wczytuje zmienne zapisane w pliku matlab.mat; load plik wczytuje zmienne zapisane w pliku plik.mat; load plik.rozsz wczytuje zmienne zapisane w pliku tekstowym o podanej nazwie i dowolnym rozszerzeniu; dane muszą tworzyć tablicę prostokątną; zostaną zapisane w macierzy o nazwie plik;

100 Przykład: who  Your variables are: A a b B ans x  save plik A save plik.dat A -ascii clear A A ??? Reference to a cleared variable A.

101 load plik.dat  who  Your variables are: B ans plik a b x  load plik.mat A a b x

102 OPERACJE NA PLIKACH (II)Przed zapisaniem lub odczytaniem danych należy otworzyć plik za pomocą funkcji fopen: id_pliku=fopen(nazwa_pliku,rodzaj_dostępu)

103 otwarcie pliku do odczytu ‘w’ Wartość argumentu opis ‘r’ otwarcie pliku do odczytu ‘w’ usunięcie zawartości istniejącego pliku lub otworzenie nowego i otwarcie go do zapisu ‘a’ otwarcie pliku w celu dopisywania elementów na jego końcu ‘r+’ otwarcie pliku do odczytu i zapisu ‘w+’ usunięcie zawartości istniejącego pliku lub utworzenie nowego i otwarcie go do odczytu i zapisu ‘a+’ otwarcie pliku w celu czytania lub dopisywania elementów na jego końcu

104 Funkcja fopen otwiera plik wskazany łańcuchem nazwa_pliku i zwraca unikatowy identyfikator pliku (zmienną id_pliku). Identyfikator ten powinien być używany we wszystkich operacjach wejścia i wyjścia wykonywanych na danym pliku. Jeśli operacja otwarcia pliku zakończy się sukcesem, zmienna id_pliku będzie nieujemną liczbą całkowitą, w przeciwnym wypadku przyjmie wartość –1.

105 Druga postać wywołania funkcji fopen jest następująca:[id_pliku,informacja]=fopen(nazwa_pliku,rodzaj_dostępu) Informacja jest łańcuchem znakowym, który może być pomocny w ustaleniu błędu. Jest on zwracany kiedy operacja otwarcia pliku zakończy się niepowodzeniem. Zamknięcie pliku o podanym identyfikatorze: status=fclose(id_pliku) Zamknięcie wszystkich otwartych plików: status=fclose(‘all’)

106 Zapisu elementów macierzy A w pliku binarnym określonym identyfikatorem id_pliku dokonujemy przy pomocy funkcji fwrite: liczba=fwrite(id_pliku,A,typ) Argument funkcji typ pozwala określić, na ilu bitach mają być zapisane dane i jak powinny być zinterpretowane. Wartością domyślną argumentu jest ‘uchar’.

107 Wartość argumentu typ f-cji fwrite Interpretacja‘uchar’ pojedynczy znak zapisany na 8 bitach bez znaku +/- ‘schar’ pojedynczy znak zapisany na 8 bitach, w tym jeden bit przeznaczony na +/- ‘int8’,’int16’,’int32’,’int64’ liczba całkowita ze znakiem zapisana odpowiednio na 8,16,32,64 bitach ‘’int8’,’int16’,’int32’,’int64’ liczba całkowita bez znaku

108 liczba zapisana w formacie zmiennopozycyjnym na 8 bitach‘single’ liczba zapisana w formacie zmiennopozycyjnym na 8 bitach ‘float32’ liczba zapisana w formacie zmiennopozycyjnym na 16 bitach ‘double’ liczba zapisana w formacie zmiennopozycyjnym na 32 bitach ‘float64’ liczba zapisana w formacie zmiennopozycyjnym na 64 bitach

109 Odczyt plików binarnychA=fread(id_pliku,rozmiar,typ) [A,liczba]=fread(id_pliku,rozmiar,typ) Funkcja fread wczytuje dane z pliku binarnego określonego przez identyfikator id_pliku i zapisuje je w macierzy A. Rozmiar określa liczbę argumentów, które powinny zostać wczytane z pliku.

110 Wartości argumentu rozmiar funkcji fread Opisodczytuje n elementów i zapisuje je w wektorze kolumnowym [m,n] odczytuje tyle argumentów, aby wypełniły kolumnami macierz o rozmiarze mxn; brakujące elementy są zastępowane zerami

111 SFORMATOWANE PLIKI TEKSTOWEFunkcja fprintf o wywołaniu: liczba=fprintf(id_pliku,format,A,...) Argument format jest łańcuchem znakowym określającym m.in..rodzaj konwersji, szerokość pola i liczbę cyfr znaczących każdej wymienionej w funkcji fprintf macierzy.

112 do zapisu liczb całkowitych %f Rodzaj konwersji Opis %d do zapisu liczb całkowitych %f do zapisu liczb rzeczywistych w formacie stałoprzecinkowym %e do zapisu liczb rzeczywistych w formacie zmiennoprzecinkowym %g automatyczny dobór krótszego formatu (%e lub %f) %c do zapisu pojedynczych znaków %s do zapisu łańcuchów znakowych

113 powrót karetki (przesunięcie kursora do początku wiersza)Znak specjalny opis \b cofnięcie o jeden znak \f nowa strona \n nowy wiersz \r powrót karetki (przesunięcie kursora do początku wiersza) \t znak tabulatora \” znak apostrofu \\ znak lewego ukośnika %% znak procentu

114 Odczyt danych z pliku tekstowego o identyfikatorze id_pliku dokonuje funkcja fscanf. Konwertuje ona dane w sposób określony w argumencie format i zapisuje w macierzy A: A=fscanf(id_pliku,format,rozmiar) [A,liczba]=fscanf(id_pliku,format,rozmiar) Rozmiar określa liczbę elementów, które powinny zostać wczytane z pliku. Może on przyjmować wartości identyczne jak w wypadku funkcji fread. Polecenie fscanf dopuszcza takie same rodzaje konwersji, jak funkcja fprintf.

115 Grafika

116 Wykresy plot plot(Y) plot(X1,Y1,...) plot(X1,Y1,LineSpec,...)plot(Y) plot(X1,Y1,...) plot(X1,Y1,LineSpec,...) plot(...,'Nazwa _cechy', Wartość_cechy,...) h = plot(...) x = -pi:pi/10:pi; y = tan(sin(x)) - sin(tan(x)); plot(x,y)

117 Oznaczenia kolorów: y yellow m magenta c cyan r red g green b blueOznaczenia kolorów: y yellow m magenta c cyan r red g green b blue w white k black Kształt punktów s square (prostokąt) d diamond (karo) v trójkąt w dół ^ trójkąt w góre < trójkąt w lewo > trójkąt w prawo p pentagram h hexagram Linia wykresu . punktowy - ciągła o kółka : kropkowana -. kropka-kreska -- przerywana x x-owana + plus * gwiazdki

118 plot(x,y,'-rs','LineWidth',2,... 'MarkerEdgeColor','k',... 'MarkerFaceColor','g',... 'MarkerSize',10) ;

119 Opis osi plot(x,sin(x)) xlabel('Oś X') ylabel('Oś Y')

120 plot(x,sin(x)); xlabel('-\pi \leq \Theta \leq \pi','FontSize',60) ylabel('sin(\Theta)','FontSize',60)

121 legend(h,'string1','string2',...) legend('string1','string2',...) Title text gtext legend legend(h,'string1','string2',...) legend('string1','string2',...) h = legend(...) x = -pi:pi/20:pi; plot(x,cos(x),'-ro',x,sin(x),'-.b') h = legend('cos','sin',2);

122 x = -pi:pi/20:pi; plot(x,cos(x),'-ro',x,sin(x),'-.b') h = legend('cos','sin',2); Uwaga: legendę można edytować I przesuwać. Podobnie jak wszystkie elementy wykresu !!

123 set x = -pi:.1:pi; y = sin(x); plot(x,y) set(gca, 'FontSize',30)x = -pi:.1:pi; y = sin(x); plot(x,y) set(gca, 'FontSize',30) set(gca,'XTick',-pi:pi/2:pi) set(gca,'XTickLabel',{'-pi','-pi/2','0','pi/2','pi'}) xlabel('-\pi \leq \Theta \leq \pi') ylabel('sin(\Theta)') title('Plot of sin(\Theta)') text(-pi/4,sin(-pi/4),'\leftarrow sin(-\pi\div4)',... 'HorizontalAlignment','left')

124 set(findobj(gca,'Type','line','Color',[0 0 1]),... 'Color','red',... h=plot(x,sin(x)) set(findobj(gca,'Type','line','Color',[0 0 1]),... 'Color','red',... 'LineWidth',2) h =

125 XLimMode, YLimMode, and ZlimModex=[0:.1:7]; plot(x,sin(x)) set(gca, 'Xlim',[0,1])

126 set(gca,'Xlabel',text('String','axis label')) XGrid, YGrid, ZGrid on|{off} włączenie/wyłączenie siatki set(gca,'XGrid','on') XScale, YScale, ZScale {linear} | log set(gca,'XScale','log') => semilogx(x,sin(x)) XTickLabel, YTickLabel, ZtickLabelstring set(gca,'XTickLabel',{'One';'Two';'Three';'Four'}) CameraPosition [x, y, z] axes coordinates set(gca,'CameraPosition',[1,3,4]) CLim - [cmin, cmax] - Color axis limits XColor, YColor, ZColor ColorSpec XDir, YDir, ZDir {normal} | reverse set(gca,'XDir','reverse')

127 set(gca) ALim ALimMode: [ {auto} | manual ] AmbientLightColor Box: [ on | {off} ] CameraPosition CameraPositionMode: [ {auto} | manual ] CameraTarget CameraTargetMode: [ {auto} | manual ] CameraUpVector CameraUpVectorMode: [ {auto} | manual ] CameraViewAngle CameraViewAngleMode: [ {auto} | manual ] CLim CLimMode: [ {auto} | manual ] Color ColorOrder DataAspectRatio DataAspectRatioMode: [ {auto} | manual ] DrawMode: [ {normal} | fast ] etc...

128 BackingStore: [ {on} | off ] set(gcf) Alphamap BackingStore: [ {on} | off ] CloseRequestFcn: string -or- function handle -or- cell array Color Colormap CurrentAxes CurrentCharacter CurrentObject CurrentPoint Dithermap DithermapMode: [ auto | {manual} ] DoubleBuffer: [ on | {off} ] ... etc

129 axis axis([xmin xmax ymin ymax])axis([xmin xmax ymin ymax]) axis([xmin xmax ymin ymax zmin zmax cmin cmax]) v = axis axis auto – sam dobiera wszystkie parametry axis manual – "zamraża" aktualne parametry do następnych wykresów axis tight – zakres dokł. z danymi axis fill – zostawiona przestrzeń axis ij – odwrócenie do numeracji "macierzowej" (wiersz, kolumna) axis xy – powrót do norm. axis equal axis image – spłaszczenie w pionie axis square – w kwadracie axis vis3d axis normal axis off - wyłączenie axis on [mode,visibility,direction] = axis('state')

130 subplot(2,1,1); subplot(2,1,2);Kilka wykresów na jednym układzie współrzędnych Hold on / hold off figure figure(numer) subplot subplot(m,n,p) subplot(h) subplot('Position',[left bottom width height]) h = subplot(...) in = [ ]; out = [ ]; subplot(2,1,1); plot(in) subplot(2,1,2); plot(out)

131 Inne typy wykresów semilogx semilogy x = 0:.1:10; semilogy(x,10.^x)

132 loglog x = logspace(-1,2); loglog(x,x.^7,'-s') grid on

133 Bar Y = round(rand(5,3)*10); subplot(2,2,1) bar(Y,'group')subplot(2,2,1) bar(Y,'group') title 'Group' subplot(2,2,2) bar(Y,'stack') title 'Stack' subplot(2,2,3) barh(Y,'stack') subplot(2,2,4) bar(Y,1.5) title 'Width = 1.5'

134 hist x = -2.9:0.1:2.9; y = randn(10000,1); hist(y,x)hist N = HIST(Y) – dzieli na 10 równomiernie rozłożonych przedziałów. Jeśli Y macierz – działa w dół kolumny N = HIST(Y,M), M - skalar, dzieli na M słupków. n = hist(Y,x), x-wektor [n,xout] = hist(...) – n-zliczenia i położenie słupków (xout). można potem bar(xout,n) aby narysować histogram. x = -2.9:0.1:2.9; y = randn(10000,1); hist(y,x)

135 h = findobj(gca,'Type','patch');set(h,'FaceColor','r','EdgeColor','w')

136 Cień line przykład – tworzenie cienia t = 0:pi/20:2*pi;t = 0:pi/20:2*pi; hline1 = plot(t,sin(t),'k'); hline2 = line(t+.06,sin(t),'LineWidth',4,'Color',[ ]); %wysunięcie linii na czoło set(gca,'Children',[hline1 hline2])

137 Wykresy trójwymiaroweplot3 plot3(X1,Y1,Z1,LineSpec,...) t = 0:pi/50:10*pi; plot3(sin(t),cos(t),t) grid on axis square

138 Mesh, meshgrid [X,Y] = meshgrid(-3:.125:3); Z = sin(X.^2)+cos(Y.^2);mesh(X,Y,Z); axis([ ])

139 meshc [X,Y] = meshgrid(-3:.125:3); Z = peaks(X,Y); meshc(X,Y,Z);[X,Y] = meshgrid(-3:.125:3); Z = peaks(X,Y); meshc(X,Y,Z); axis([ ])

140 meshz [X,Y] = meshgrid(3:.125:3); Z = peaks(X,Y); meshz(X,Y,Z)

141 contour [X,Y] = meshgrid(-2:.2:2,-2:.2:3); Z = X.*exp(-X.^2-Y.^2);[X,Y] = meshgrid(-2:.2:2,-2:.2:3); Z = X.*exp(-X.^2-Y.^2); contour(X,Y,Z,20);

142 colormap [X,Y] = meshgrid(-2:.2:2,-2:.2:3); Z = X.*exp(-X.^2-Y.^2);[X,Y] = meshgrid(-2:.2:2,-2:.2:3); Z = X.*exp(-X.^2-Y.^2); [C,h] = contour(X,Y,Z); clabel(C,h) colormap cool

143 contour3 [X,Y] = meshgrid([-2:.25:2]); Z = X.*exp(-X.^2-Y.^2);contour3(X,Y,Z,30)

144 surface figure; surface(peaks)

145 Eksport wykresów do plikuprint print -device -options filename -dwin % wysyła b/c do aktualnej drukarki -dwinc % wysyła w kolorze do aktualnej drukarki -dmeta % wysyła do schowka w formacie Meta -dbitmap % wysyła do schowka w formacie bitmapa -dsetup % okienko dialogowe drukowania bez drukowania -deps % Encapsulated PostScript -depsc % Encapsulated Color PostScript -deps2 % Encapsulated Level 2 PostScript -depsc2 % Encapsulated Level 2 Color PostScript

146 Animacje - getframe getframe F = getframe F = getframe(h)F = getframe F = getframe(h) F = getframe(h,rect) [X,Map] = getframe(...) movie movie(M) movie(M,n) movie(M,n,fps) movie(h,...) movie(h,M,n,fps,loc)

147 Przykład animacji plot(fft(eye(k+1))) plot(fft(eye(k+16))) plot(fft(eye(k+16))) for k = 1:16 plot(fft(eye(k+16))) axis equal M(k) = getframe; end movie(M,30)

148 Drawnow, imagesc % Przykład – potencjał czynnościowy sercaload napiecie clim=[-95,60] h=imagesc(napiecie(:,:,1),'CDataMapping','scaled','EraseMode','xor','Clipping','on',clim); colorbar [a,b,inc]=size(napiecie) for i=[2:inc] set(h,'CData', napiecie(:,:,i),'CDataMapping','scaled' ); drawnow; end

149 Potencjał czynnościowy serca

150 GUI

151 Obiekty graficzne uicontrol uicontrol('propertyname1',value1,'propertyname2,'value2,...) h=uicontrol('Style','pushbutton','String','przycisk') set(hh) BackgroundColor Callback: string -or- function handle -or- cell array CData Position String Style: [ {pushbutton} | togglebutton | radiobutton | checkbox | edit | text | slider | frame | listbox | popupmenu ]

152 h=uicontrol('Style','pushbutton','String','przycisk')

153 Callback 'Style','pushbutton', ... 'Units','normalized', ...stopHndl=uicontrol( ... 'Style','pushbutton', ... 'Units','normalized', ... 'Position',[xPos yPos-spacing btnLen btnWid], ... 'Enable','off', ... 'String','Stop', ... 'Callback',‘moja_funkcja(X,Y)');

154 uimenu uimenu('PropertyName',PropertyValue,...)uimenu(parent,'PropertyName',PropertyValue,...) handle = uimenu('PropertyName',PropertyValue,...) handle = uimenu(parent,'PropertyName',PropertyValue,...)

155 Przykład uimenu f = uimenu('Label','Workspace');uimenu(f,'Label','New Figure','Callback','figure'); uimenu(f,'Label','Save','Callback','save');

156 uigetfile / uiputfile uigetfile uigetfile('FilterSpec')uigetfile('FilterSpec','DialogTitle') uigetfile('FilterSpec','DialogTitle',x,y) [fname,pname] = uigetfile(...) [fname,pname] = uigetfile('*.m',’Podaj nazwe pliku)    uiputfile

157

158

159

160

161 Plik źródłowy generowany przez GUIUWAGA: typ generowanego pliku zmienia się w kolejnych wersjach MATLABA

162

163 % --- Executes on button press in pushbutton1. function pushbutton1_Callback(hObject, eventdata, handles) % hObject handle to pushbutton1 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) disp('napisz cos tam') plot([1 2 3], 'Parent',handles.axes1) % --- Executes during object creation, after setting all properties.

164

165 % --- Executes during object creation, after setting all properties.function edit1_CreateFcn(hObject, eventdata, handles) if ispc set(hObject,'BackgroundColor','white'); else set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor')); end function edit1_Callback(hObject, eventdata, handles) % hObject handle to edit1 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Hints: get(hObject,'String') returns contents of edit1 as text % str2double(get(hObject,'String')) returns contents of edit1 as a double

166 Grafika cd

167 inspect Alternatywa dla set(h) lub get(h) >> h=plot(x,y);>> inspect(h)

168

169 Statystyka danych

170 Kursor danych

171 Paleta

172 Wykresy 3D - cd Przykład: [X,Y] = meshgrid(-8:.5:8);R = sqrt(X.^2 + Y.^2) + eps; Z = sin(R)./R; surf(X,Y,Z,'FaceColor','interp',... 'EdgeColor','none',... 'FaceLighting','phong') daspect([5 5 1]) % osie axis tight view(-50,30) camlight left

173 Okna dialogowe - inputdlgAddOpts.Resize='on'; AddOpts.WindowStyle='normal'; AddOpts.Interpreter='tex'; textLabel={'opis1', 'opis2‘, ...} 'Nr of open neighbors to switch into open state (Noo) [0-4]',... newParam=inputdlg(textLabel, windowTitle, ile_linii,... startparameters,AddOpts);

174

175 FUNKCJE ZARZĄDZAJĄCE LISTĄ ŚCIEŻEKpath wyświetla listę ścieżek; path(path,s1) dodaje do listy ścieżek katalog określony łańcuchem s1; path(s1) zmienia listę ścieżek na składającą się tylko z katalogu określonego łańcuchem s1;

176 PRZESZUKIWANIE M-plikówtype nazwa_funkcji lookfor random path (path, "C:\…….") which psd C:\MATLABR12\toolbox\signal\signal\psd.m

177 Operacje plikowe – cd. Open nazwa.fig (lub inne rozszerzenie, np. .pdf) – uruchamia z aplikacji domyślnej dlmwrite A=['ala']; >> dlmwrite (‘plik.txt', A) -> a,l,a >> dlmwrite ('cos.txt', A, '') -> ala >> save ‘plik.txt' A –ascii e e e+001

178 xlsread – import danych z arkusza ExcelaA = xlsread('testdata2.xls', 1, 'A4:B5') % arkusz 1, komórki A4:B5 xlswrite d = {'Time', 'Temp'; 12 98; 13 99; 14 97}; s = xlswrite('tempdata.xls', d, 'Temperatures', 'E1') textscan – czyta plik tekstowy, zawartość do cell (o typie cell za chwilę) importdata s = importdata('ding.wav') s =data: [11554x1 double] fs: 22050

179 imread – wczytanie grafiki >> a=imread('niebo.jpg'); >> image(a) % lub imagesc(a) – skaluje wykres do pełnej mapy kolorów; Rysunek w oknie wykresów Matlaba imwrite – zapis grafiki BW = imread('text.png'); imwrite(BW,'test.tif'); UWAGA: Przed nazwą pliku zawsze może być adres URL

180

181 Operatory porównania: <, <=, >, >=, ==, ~=Funkcje logiczne  A|B A&B xor(A,B) ~A Operatory porównania: <, <=, >, >=, ==, ~= UWAGA: Podwójny operator logiczny oznacza operacje bitowe: A && B A || B

182 Funkcje różne –cd.

183 Funkcje obsługujące zarządzanie pamięcią.Przestrzeń robocza-obszar pamięci, w której przechowywane są zmienne utworzone w oknie poleceń. who wyświetla listę wszystkich zmiennych znajdujących się aktualnie w pamięci; whos wyświetla listę wszystkich zmiennych wraz z informacją na temat ich rozmiaru i rodzaju; who global whos global wyświetlają informacje o zmiennych globalnych; clear usuwa z pamięci wszystkie zmienne; clear z usuwa z pamięci zmienną z; clear global z usuwa zmienną globalną z; clear all usuwa wszystkie zmienne i funkcje;

184 Funkcje obsługujące polecenia systemowe: dir lub ls cd katalog Katalog bieżący - zapisywane są w nim pliki tworzone podczas pracy z pakietem. Funkcje obsługujące polecenia systemowe: dir lub ls wyświetla pliki w bieżącym lub podanym katalogu (dozwolone jest użycie znaków masek *,?) cd katalog zmienia katalog bieżący na podany delete plik usuwa plik o podanej nazwie pwd wyświetla pełną ścieżkę określającą katalog bieżący !polecenie wykonuje dowolne polecenie systemu operacyjnego (np.. !type matlab.txt wyświetli zawartość pliku matlab.txt)

185 FUNKCJE OBSŁUGUJĄCE POMIAR CZASUclock podaje aktualną datę i czas w postaci sześcioelementowego wektora [rok miesiąc dzień godzina minuta sekunda] date podaje aktualną datę w postaci łańcucha o formacie:’dd-mmm-rrrr’ etime(t2,t2) podaje różnicę czasu, który upłynął między chwilami t1 i t2 (t1, t2-wektory o formacie, jak w poleceniu clock) tic zeruje odmierzanie czasu przed użyciem polecenia toc toc podaje czas (w sek.), który upłynął od momentu użycia polecenia tic

186 INNE TYPY DANYCH

187 TYP CELL Składowymi zmiennej typu CELL są macierze, łańcuchy, inne obiekty typu cell, zmienne obiektowe języka Java, cell(object_Java). Deklaracja: np. cell(1,2) >> A(1,1) = {magic(5)} A = [5x5 double] >> A(1,2) = {{[5 2 8; 7 3 0; 6 7 3] 'Test 1'; [2-4i 5+7i] {17 []}}} [5x5 double] {2x2 cell}

188 Dostęp do elementów w komórceStosujemy nawiasy klamrowe >> A{1,2} ans = [3x3 double] 'Test 1' [1x2 double] {1x2 cell} >> A{1,2}{1,1}

189 Pełne wyświetlenie zawartości komórki, celldisp>> celldisp(A) A{1} = A{2}{1,1} = A{2}{2,1} = i i A{2}{1,2} = Test 1 A{2}{2,2}{1} = 17 A{2}{2,2}{2} = []

190 Graficzne wyświetlenie komórki - cellplotcellplot(A) UWAGA: Jeżeli zmienna A była raz użyta jako cell to nie wystarczy później podstawić jako zwykłą macierz aby usunąć pamięć o typie cell, „clear A” jest konieczne.

191 Typ struct s = struct('field1', values1, 'field2', values2, ...)s = struct('field1', {}, 'field2', {}, ...) struct([]), struct(Java_obj)

192 Przykłady struktur >> s = struct('type', {'big','little'}, 'color', {'red'}, 'x', {3 4}) s = 1x2 struct array with fields: type color x >> s.x ans = 3 4

193 >> s = struct('type', {'big'}, 'color', {'red'}, 'x', {3})>> s.type='small' type: 'small'

194 cell2struct >> c = {'Ania','Kowalska',165; 'Ewa','Kwiatkowska',170} c = 'Ania' 'Kowalska' [165] 'Ewa' 'Kwiatkowska' [170] >> pola = {'imie', 'nazwisko', 'wzrost'} pola = 'imie' 'nazwisko' 'wzrost' >> s = cell2struct(c, pola, 2) s = 2x1 struct array with fields: imie nazwisko wzrost

195 >> s(1) ans = imie: 'Ania' nazwisko: 'Kowalska' wzrost: 165 >> s(1).nazwisko Kowalska

196 Typ sparse (rzadki) S = sparse(A) S = sparse(m,n) Przykład:[i,j,s] = find(S); [m,n] = size(S); S = sparse(i,j,s,m,n); Jeśli ostatni wiersz i kolumna mają niezerowe elementy: S = sparse(i,j,s);

197 Funkcje obsługujące typ rzadkitf = issparse(S) – czy jest? A = full(S) – konwersja na pełny A=spy(S) – jaki wzór OSZCZĘDNOŚĆ PAMIĘCI whos Name Size Bytes Class M_full x double array M_sparse x sparse array Grand total is elements using bytes

198 >> A=round(rand(5,4))Przykład >> S=sparse(A) S = (3,1) (5,1) (1,2) (3,2) (4,2) (5,2) (5,3) (2,4) (3,4) (5,4) >> A=round(rand(5,4)) A = >> spy(S)

199 Łańcuchy Łańcuchy znakowe A='to jest napis'   A = to jest napis

200 Porównywanie łańcuchów k = strcmp('str1','str2') TF = strcmp(S,T)  Przykłady:  strcmp('Yes','No')  ans = strcmp('Yes','Yes') ans = 1

201 Konwersja na łańcuch mat2str - zamiana macierzy w łańcuch. Przydatne w poleceniu „eval” str = mat2str(A) str = mat2str(A,precyzja), precyzja – ilość cyfr po przecinku    num2str – zamiana liczby na string str = num2str(A) str = num2str(A,precision) str = num2str(A,format) int2str  Przykład: title(['Temperatura ' , int2str(T)])

202 Konwersja z łańcucha str2mat(A) – działa tylko łącznie z funkcją evalchar(T1,T2,T3) – zamiana łąńcuchów na macierz, uzupełnienie spacjami do właściwego wymiaru. W przypadku standardowych łańcuchów działa tak samo jak str2mat. Działa tez dla typu cell >> T1='to jest napis'; >> T2='to jest inny napis'; >> T3='Ten tez'; >> char(T1,T2,T3) ans = to jest napis to jest inny napis Ten tez

203 >> cc={'text1','napis_jakis','napis_inny'}>> conv=char(cc) conv = text1 napis_jakis napis_inny >> conv(2,4) ans = i

204 lower - Duże na małe literyupper (odwrotnie) strtok Pierwszy znak w łańcuchu Przykład s = ' To jest przykład.'; [token,rem] = strtok(s) token = To jes rem = t przykład

205 strrep Znalezienie i zamiana łańcucha str = strrep(str1,str2,str3). Przykład: s1 = 'To jest dobry przykład.'; str = strrep(s1,‘dobry',‘zły') str =To jest zły przykład.  findstr Znajdowanie łańcucha w łańcuchu k = findstr(str1,str2)  str1 = 'Find the starting indices of the shorter string.'; str2 = 'the'; findstr(str1,str2)   ans =

206 eval Uruchomienie łańcucha z wyrażeniem Matlaba[a1,a2,a3,...] = eval(‘moja_funkcja(b1,b2,b3,...)’) eval('[a1,a2,a3,...] = moja_funkcja(var)') % NIE POLECANE Przykład >> A='[1 2 3]' A = [1 2 3] >> A=str2mat(A) >> A(1) ans = [

207 feval fplot >> A=eval(str2mat(A)) A = 1 2 3 >> A(1) >> A(1) ans = 1 „A” odzyskuje format matematyczny feval feval('sin',pi/2) – uchwyt do funkcji sinus fplot fplot('tanh',[x0 x_kon])

208 FUNKCJE PRZETWARZAJĄCE ŁAŃCUCHYdeblank(s) usuwa spacje z końca łańcucha; findstr(s1,s2) szuka krótszego z łańcuchów s1 i s2 w dłuższym; zwraca wektor indeksów, od których zaczyna się występowanie krótszego łańcucha; lower(s) / upper(s) zmienia wszystkie litery w łańcuchu na małe / duże; strcat(s1,s2,...) łączy łańcuchy w poziomie z pominięciem spacji na końcu każdego z nich; strcmp(s1,s2) porównuje dwa łańcuchy; jeśli są identyczne, zwraca 1, jeśli nie-0; funkcja rozróżnia wielkość liter;

209 strncmp(s1,s2,n) strvcat(s1,s2,s3) upper(s) strcmpi(s1,s2)porównuje dwa łańcuchy bez rozróżniania wielkości liter; strncmp(s1,s2,n) porównuje n pierwszych znaków w dwu łańcuchach; strvcat(s1,s2,s3) łączy łańcuchy w pionie, dodając na końcu każdego z nich odpowiednią liczbę spacji; zwraca macierz znakową; upper(s) zmienia wszystkie litery w łańcuchu na duże;

210 Interpolacja

211 interp1 interp1 – interpolacja jednowymiarowa yi = interp1(x,Y,xi)yi = interp1(Y,xi) yi = interp1(x,Y,xi,method) yi = interp1(x,Y,xi,method,'extrap') x – wektor argumentów (o znanych wartościach) Y- wektor wartości dla x xi – wektor wartości do interpolacji

212 Metody 'nearest' – metoda najbliższego sąsiada'linear' – liniowa, domyślna 'spline'- krzywa sześcienna składana 'pchip', 'cubic' (to samo) –sześcienna, zachowująca monotoniczność i kształt danych 'v5cubic' -Cubic interpolation używany w MATLAB 5

213 Przykłady Przykład: Interpolacja niedokładnej sinusoidy x=0:10;x=0:10; y=sin(x); xi=0:0.25:10; yi=interp1(x,y,xi); plot(x,y,'o',xi,yi)

214 %interp_metoda x=0:10; y=sin(x); xi=0:0.25:10; yi=interp1(x,y,xi,'linear'); yi2=interp1(x,y,xi,'spline'); yi3=interp1(x,y,xi,'cubic'); yi4=interp1(x,y,xi,'nearest'); plot(x,y,'o',xi,yi,'-', xi,yi2,':',xi,yi3,'-.',xi,yi4,'--') legend('punkty','liniowa','spline','cubic','nearest')

215 Spline    Interpolacja wielomianem wysokiego rzędu – często błędne wyniki (szczególnie poza przedziałem danych). Rózne metody uniknięcia. Metoda cubic spline –każda para punktów ('końcowe') – inny wielomian 3 s-nia. Nieskończona ilość krzywych więc I i II pochodna musi się zgadzać w punktach końcowych. Nachylenie i krzywizna wielomianów musi być ciągła w tych punktach. Dla n punktów rozwiązywanych jest 4(n-1) równań. Jeżeli obliczenia wielokrotnie, ale dla różnych danych podstawowych x– spline bez ostatniego argumentu xi (spline oblicza formę pp), którą można wywoływać wielokrotnie. Wykorzystanie f-cji ppval. Funkcja 'spline' b. dobra, gdy dane źródłowe są z natury dosyć gładkie. W przeciwnym razie czestym błędem jest generowanie nieistniejących maks. i min. Funkcja pchip mniej dokładna dla danych gładkich, ale mniej oscylacyjna

216  Przykład – interp_pchip.mx=[ ]; y=exp(-x/6).*cos(x); cs=spline(x,y); %cubic spline ch=pchip(x,y) %cubic Hermite xi=linspace(0,10); ysi=ppval(cs,xi) %spline yci=ppval(ch,xi) %cubic plot(x,y,'o',xi,ysi,':',xi,yci); legend('dane','spline','cubic')

217 Przykład spline Przykład interp_spline.m: x=0:10; y=sin(x);interp( ...'spline') można zastąpić spline Przykład interp_spline.m: x=0:10; y=sin(x); pp=spline(x,y) xi=0.1:0.5:10; yi=ppval(pp,xi) plot(x,y,'o',xi,yi)

218 Interpolacja dwuwymiarowainterp2 ZI = interp2(X,Y,Z,XI,YI) ZI = interp2(Z,XI,YI) ZI = interp2(Z,ntimes) ZI = interp2(X,Y,Z,XI,YI,method)

219 Przykład [X,Y] = meshgrid(-3:.25:3); Z = X.^2+Y.^2;[XI,YI] = meshgrid(-3:.125:3); ZI = interp2(X,Y,Z,XI,YI); figure mesh(X,Y,Z), hold, mesh(XI,YI,ZI+15) hold off axis([ ]) Uwaga: Jeżeli punkty nie pojawiaja się na siatce geometrycznej w regularny sposób wówczas interpolację przeprowadza się w oparciu o funkcję dalaunay, tsearch, dsearch

220 Metody interpolacji 3D 'nearest' – metoda najbliższego sąsiada‘bilinear' – liniowa, domyślna ‘bicubic‘ – sześcienna, zachowująca monotoniczność i kształt danych

221 [x,y] = meshgrid(-3:1:3); z = peaks(x,y); surf(x,y,z) [xi,yi] = meshgrid(-3:0.25:3); zi1 = interp2(x,y,z,xi,yi,'nearest');

222 Dopasowywanie krzywej wielomianem (ang. fit)Metoda najmniejszych kwadratów. Wybieramy stopień wielomianu n. Jeżeli n=1 – regresja liniowa. polyfit p=polyfit(x,y,n) p-wsp. wielomianu od najwyższej potęgi n-stopień wielomianu

223 Przykład Przykład fitowanie.m x=[0:0.1:1];x=[0:0.1:1]; y=[-.0447, 1.978, 3.28, 6.16, 7.08, , 7.66, 9.56, 9.48, 9.3, 11.2]; n=2; p=polyfit(x,y,n) p = Czyli: y(x)= x x

224 polyval Za pomocą funkcji polyval (oblicza wartość wielomianu o współczynnikach p) xi=x; yi=polyval(p,xi) plot(x,y,'o',xi,yi,':')

225 Narzędzia z menu okna - dopasowanie

226 Dopasowywanie krzywą dowolnego typu Funkcja fitDopasowywanie krzywą dowolnego typu Funkcja fit.m (CurveFitting Toolbox) Przykład: dopasowanie danych krzywą sigmoidalną 5 stopnia

227 opts = fitoptions('Method','Nonlinear','Normalize','Off'); % proba_fit.m E=[1:1:10]'; Y=[ ]'; plot(E,Y,'ko-','LineWidth',2); hold on %fitowanie sigmoidalne opts = fitoptions('Method','Nonlinear','Normalize','Off'); ftype = fittype('a*x.^5./(x.^5+b)','options',opts) f2 = fit(E,Y,ftype) % [fit1,gof1,out1]=fit(E,Y,'poly2') a=f2.a; b=f2.b; Y2=a*E.^5./(E.^5+b); plot(E,Y2,'k--','LineWidth',2); hold off >> proba_fit ftype = General model: ftype(a,b,x) = a*x.^5./(x.^5+b) Warning: Start point not provided, choosing random start point. f2 = General model: f2(x) = a*x.^5./(x.^5+b) Coefficients (with 95% confidence bounds): a = (5.69, 7.295) b = (474.6, 2786)

228

229 Dopasowanie krzywej – optymalizacja (fminsearch-Optimization Toolbox)

230 % fitdemo.m t = (0:.1:2)'; y = [ ]'; plot(t,y,'ro'); hold on; h = plot(t,y,'b'); hold off; title('Input data'); ylim([0 6]) start = [1;0]; options = optimset('TolX',0.1); estimated_lambda = %minimalizacja po pierwszym argumencie f-cji fitfun tzn. lambda

231 function err = fitfun(lambda,t,y,handle)%FITFUN(lambda,t,y,handle) oblicza błąd (err) pomiędzy oryginalnymi wartościami y i obliczonymi (z) na podstawie estymacji funkcją wykładniczą o nieznanych współczynnikach c(1)...c(n). Argument lambda nie jest podawany na wejście, po tym argumencie będzie zachodziła minimalizacja. % Zakładamy, ze nasze dane można przybliżyć funkcją % y = c(1)*exp(-lambda(1)*t) c(n)*exp(-lambda(n)*t) % z n parametrami liniowymi i n nieliniowymi A = zeros(length(t),length(lambda)); for j = 1:length(lambda) A(:,j) = exp(-lambda(j)*t); end c = A\y; z = A*c; err = norm(z-y); set(gcf,'DoubleBuffer','on'); set(handle,'ydata',z) drawnow pause(.04)

232 Dopasowanie krzywej – optymalizacja (fminsearch-Optimization Toolbox)

233 cftool (CurveFitting Toolbox)

234

235

236

237 Rozkład na harmoniczne Analiza w dziedzinie częstotliwości (I)>> f=[ ]; >> omega=2*pi*f; >> sygnal=a(1)*sin(omega(1).*t))+a(2)*sin(omega(2).*t))+... a(3)*sin(omega(3).*t))+a(4)*sin(omega(4).*t)); >> plot(t,sygnal)

238 Transformata FourieraKażdy sygnał sinusoidalny może być przedstawiony jako: fazor Im A   Re Ze wzorów trygonometrycznych:

239 Sygnał można przybliżać sumąDyskretna transformata Fouriera prosta i odwrotna gdzie

240 Transformata Fouriera – fft, ifftTransformata Fouriera: X = fft(x) X = fft(x,N) transformata odwrotna x = ifft(X) x = ifft(X,N) Transformata obliczana dla wektorów o długości N fft2 (dwuwymiarowa), fftn (wielowymiarowa), fftshift (przesuwa zerową częstotliwość w celu wycentrowania widma).

241 function [f,y]=fourier_sampling (A,B, f_sampling,f1,f2,ile_punktow)Przykład fft function [f,y]=fourier_sampling (A,B, f_sampling,f1,f2,ile_punktow) % fourier_sampling (0,100, 200,10,50); t = A:1/f_sampling:B; x = sin(2*pi*f1*t) + sin(2*pi*f2*t); if nargin>5 y = fft(x,ile_punktow); T=length(y)/f_sampling; % t2=[A:1/f_sampling:(T+A-1/f_sampling)]; else y = fft(x); T=B-A; % t2=t; end m = abs(y); f =(0:length(y)-1)/T; subplot(4,1,1), plot(t,x), xlabel('czas'), ylabel('sygnal'), subplot(4,1,2), plot(f,m), xlabel('czestotliwosc'), ylabel('|fft|'), grid on

242 Przykład ifft % ---- filter 50 Hz (ciąg dalszy tej samej funkcji)eps=1e1; gdzie1=find(abs(f-50) gdzie2=find(abs(f-(f_sampling-50)) y([gdzie1,gdzie2])=0; subplot(4,1,3), plot(f, abs(y)) title('PO FILTRACJI','FontWeight','bold') xlabel('czestotliwosc'), ylabel('|fft|'), grid on filter_x=real(ifft(y)); subplot(4,1,4), plot(t2,filter_x,'-r',t2,sin(2*pi*f1*t2),':g') xlabel('czas'), ylabel('sygnal')

243

244 Analiza w dziedzinie czasuKowariancja E(.) – wartość oczekiwana (często przybliżana średnią) C = cov(X) C = cov(x,y) korelacja wzajemna, autokorelacja c = xcorr(x,y), c = xcorr(x)

245 Sygnały losowe Analiza w dziedzinie częstotliwości (II), spectrumSygnały losowe Analiza w dziedzinie częstotliwości (II), spectrum.periodogram – Signal Processing Toolbox Periodogram fft(Cov(Xn+m,Xn)) Fs=1000; t=0 : 1/F : 0.3 ; x=cos(2*pi*t*200)+randn(size(t)); % domyślne wartości Hs=spectrum.periodogram; psd(Hs,x,'Fs',Fs) [Pxx,w] = periodogram(x)

246 Analiza czasowo-częstotliwościowa. FalkiCOEFS = cwt(S,SCALES,'wname') COEFS = cwt(S,SCALES,'wname','plot') Przykład: czest_probk=100; t = 0:1/czest_probk:(10-1/czest_probk); %dlugosc - potega 2 czas=length(t); x = sin(2*pi*5*t) + sin(2*pi*47*t); x(1:floor(czas/2))=x(1:floor(czas/2))/2; plot(t,x) xlabel('czas [s]') ylabel('sin(2*pi*5*t) + sin(2*pi*47*t)') figure coef=cwt(x,1:100,'db4'); surf(coef)

247

248 Dźwięk

249 Nagrywanie dźwięku - wavrecordy = wavrecord(n,Fs) y = wavrecord(...,ch) y = wavrecord(...,'dtype') n – ilość próbek , Fs – częstotliwość próbkowania [próbek/s], (domyślnie Fs = Hz), ch – ilość kanałów we: 1-mono (domyślne), 2-stereo, 'dtype' – typ danych: 'double' (domyślne) - 16 bitów/próbka, 'single‘ bitów/próbka, 'int16‘ - 16 bitów/próbka

250 Odtwarzanie dźwięku - wavplaywavplay(y,Fs) wavplay(...,'mode') ‘mode’: ‘async’ (domyślne) – utrzymany dostęp do linii poleceń ‘sync’ (brak dostępu do linii poleceń w trakcie odtwarzania) Przykład: % nagrywamy 16-bitowy dźwięk przez 5 s. % Częstość próbkowania Fs=11025 Hz., (ilość_próbek=Czas_nagrania * Fs) Fs = 11025; y = wavrecord(5*Fs,Fs,'int16'); % y = wavread('filename') % dźwięk może być też wczytany z pliku % odtwarzamy dźwięk wavplay(y,Fs);

251 wavwrite(y,Fs,'filename') wavwrite(y,Fs,N,'filename')Zapis dźwięku do pliku Nagrany dźwięk można poddać przeróbce (np. filtracji) i tak przetworzony zapisać do pliku wavwrite(y,'filename') wavwrite(y,Fs,'filename') wavwrite(y,Fs,N,'filename')

252 Signal Processing Toolbox –spectrum.music

253 sptool – Signal Processing Toolbox (Filtracja, przetwarzanie, widma)

254 Okno -Filter Design

255 Okno – Signal Browser

256 Problemy algebry liniowej w Matlabie   Metody rozwiązywania układów równań liniowych o właściwej, nadmiernej lub niedostatecznej ilości równań.

257 Identyczna ilość równań i zmiennych.

258 Dzielenie lewostronneAlgebra liniowa stosuje wiele różnych metod do rozwiązania układu równań liniowych np. eliminacji Gaussa, faktoryzacji LU, faktoryzacji Cholesky’ego, Cramera i in. Matlab automatycznie stosuje metodę najbardziej optymalną dla podanej macierzy, wprowadzając operację dzielenia lewostronnego \ Aby być pewnym, że otrzymane rozwiązanie będzie bliskie rzeczywistości można najpierw obliczyć funkcję warunku cond(A). Jeżeli w wyniku nie otrzymamy liczby bardzo dużej (najlepiej jeżeli jest bliska 1), to macierz A-1 ma dobre własności numeryczne i Matlab potrafi dobrze rozwiązać równanie A*XX=b. Inną metodą (mniej polecaną) niż dzielenie lewostronne jest obliczenie macierzy odwrotnej A-1= inv(A).

259 cond >> A=[5 7 4; 8 2 -3; 1 1 -1] A = 5 7 4 8 2 -3 1 1 -1 >> cond(A) ans =

260 Metoda dzielenia lewostronnego>> XX=A\b XX= Sprawdzenie: >> wyn=A*XX wyn=3 2 1

261 Metoda macierzy odwrotnej(metoda wolniejsza) XX=A-1*b >>XX=inv(A)*b XX= 0.6875

262 Metoda macierzy pseudoodwrotnej - pinvMetoda kosztowna, stosowana jedynie wówczas, gdy nie można wyznaczyć macierzy odwrotnej. Na ogół w sytuacjach, o których będzie mowa poniżej. >> XX=pinv(A)*b

263 Metoda Cramera >> A=[5 7 4; 8 2 -3; 1 1 -1] >> b=[3 2 1]'>> AX=[b,A(:,2),A(:,3)] AX = >> AY=[A(:,1),b,A(:,3)] >> AZ=[A(:,1),A(:,2),b]

264 > wyznacznikA=det(A)64 >> wyznacznikAX=det(AX) wyznacznikAX = -4 >> wyznacznikAY=det(AY) wyznacznikAY = 44 >> wyznacznikAZ=det(AZ) wyznacznikAZ = -24

265 >> X=wyznacznikAX / wyznacznikA >> Y=wyznacznikAY / wyznacznikA Y = 0.6875 >> Z=wyznacznikAZ / wyznacznikA Z =

266 Metody rozwiązywania układów równań liniowych o nadmiernej ilości równań.W takiej sytuacji nie ma jednoznacznego rozwiązania (prostą wyznaczają dokładnie 2 punkty). Potrzeba rozwiązania takiego problemu pojawia się wówczas, gdy ilość danych pomiarowych jest większa niż ilość niewiadomych. Wszystkie dane są obarczone błędami (pomiarowymi) i brak podstaw do wskazania, które z nich należy odrzucić. Przykładowo mierzymy zależność y(x). Wiemy, ze jest ona liniowa, ale nie znamy jej równania. Na podstawie wielokrotnych pomiarów chcemy wyznaczyć równanie prostej obarczonej błędem w minimalnym stopniu. Mamy więcej zmiennych niż potrzeba (wystarczyłyby 2 punkty, przy założeniu, że nie są obarczone błędem).

267 Rozwiązanie można uzyskać dwoma metodami:

268 A. Dzielenie lewostronne (metoda najmniejszych kwadratów)Poszukuje się rozwiązania średniokwadratowego, dla którego |A*X-b| jest minimalna. Jest to rozwiązanie metodą najmniejszych kwadratów. (Zakładamy, że A*X=b) >> X=A\b

269 Przykład >> A=[1 2 3; 4 5 6; 7 8 0; 2 5 8];>> b=[366, 804,351,514]’; % Zakładamy, że b=A*X. Mamy 4 równania i 3 niewiadome >> X=A\b x= 114.93 >> e=A*X-b % Obliczamy błąd e= 0.0000 >> norm(e) ans=

270 B. Za pomocą macierzy pseudoodwrotnej (metoda Moore’a – Penrose’a).W metodzie tej wykorzystywana jest funkcja pinv (ang. pseudoinverse). >> X=pinv(A) Obliczona w ten sposób macierz X (pseudoodwrotna do X) spełnia następujące warunki: 1. A*X*A=A 2. X*A*X=X 3. A*X oraz X*A są macierzami hermitowskimi (B jest hermitowska, jeżeli B=B* i B* oznacza macierz sprzężoną)

271 Przykład A=[1 0 0; 0 1 0; 0 0 1; -1 1 0; 0 –1 1; -1 0 1]X=[x1, x2, x3]’ b=[ ]’ A*X=b Rozwiązanie: >> p=pinv(A)

272 p= Columns 1 through 4 Columns 5 through 6 >> X=p*b % odpowiednik A-1*b dla równań dobrze określonych X=1.2500 1.7500 3.0000

273 Porównujemy z dzieleniem lewostronnym:>>X=A\b X= 1.2500 1.7500 3.0000 Otrzymany wynik jest identyczny.

274 Przykładowe zadanie Wielkość y jest funkcją liniową x. Nie znamy równania tej funkcji. Aby je wyznaczyć wykonujemy serię pomiarów (x,y), obarczonych przypadkowym błędem pomiarowym, aby na tej podstawie wyznaczyć równanie prostej y=a*x+c. Uzyskujemy układ n równań: y1=a*x1+c y2=a*x2+c ……. yn=a*xn+c W tym układzie równań nieznane są wartości a, c.

275 ROZWIĄZANIE Układ równań można zapisać macierzowo jako:

276 function nadokr(m1,m2) %Wyznaczanie prostej na podstawie nadokreslonego ukladu równan % (np. wyznaczenie prostej z wielu punktow pomiarowych) x=[1:0.1:7]; y=x*m1+m2; %okreslamy prosta o rownaniu y=m1*x+m2 % obliczamy wspolczynniki dla prostej nie zaburzonej A=[x',ones(size(x,2),1)]; b=y'; disp('Metoda dzielenia lewostronnego - prosta nie zaburzona') Q=A\b; m1=Q(1), m2=Q(2) disp('Metoda macierzy pseudoodwrotnej - prosta nie zaburzona') Q=pinv(A)*b; m1b=Q(1), m2b=Q(2)

277 z=y+randn(1,size(x,2)); % zaburzamy prostą szumem Gaussa (symulacjaz=y+randn(1,size(x,2)); % zaburzamy prostą szumem Gaussa (symulacja % błędu przypadkowego w pomiarach) A=[x',ones(size(x,2),1)]; b=z'; disp('Metoda dzielenia lewostronnego - prosta zaburzona') Q=A\b; m1=Q(1), m2=Q(2) disp('Metoda macierzy pseudoodwrotnej - prosta zaburzona') Q=pinv(A)*b; m1b=Q(1), m2b=Q(2) plot(x,z,'o',x,m1*x+m2,'--r',x,m1b*x+m2b,':g') legend('punkty','dzielenie lewostronne','macierz pseudoodwrotna',0)

278 Wynik >> nadokr(1,2)Metoda dzielenia lewostronnego - prosta niezaburzona m1 =1.0000 m2 =2.0000 Metoda macierzy pseudoodwrotnej - prosta niezaburzona m1b =1.0000 m2b =2.0000 Metoda dzielenia lewostronnego - prosta zaburzona m1 =0.9063 m2 =2.2784 Metoda macierzy pseudoodwrotnej - prosta zaburzona m1b =0.9063 m2b =2.2784

279

280 Metody rozwiązywania układów równań liniowych o niedostatecznej ilości równań.W takiej sytuacji mamy więcej zmiennych niż równań. Taki układ równań nie ma rozwiązania jednoznacznego. Można jednak próbować DOMNIEMYWAĆ możliwe rozwiązanie, przy pewnych założeniach dodatkowych. Rozwiązanie można uzyskać dwoma metodami: a) Dzielenie lewostronne. Znajduje rozwiązanie, które ma maksymalną liczbę zer w wyznaczanym wektorze X b) Za pomocą macierzy pseudoodwrotnej Funkcja pinv(A) znajduje rozwiązanie, w którym długość (norma) wyznaczanego wektora X jest mniejsza niż w innych możliwych rozwiązaniach. Jest to tzw. rozwiązanie o minimalnej normie i ma ono duże znaczenie praktyczne.

281 Przykład >> A=[ ; ; ]; % 4 niewiadome i 3 równania >> b=[366,804,351]'; % metodą dzielenia lewostronnego >> X=A\b  X = 0 >> norm(X) ans =

282 % metodą macierzy pseudoodwrotnej>> Xn=pinv(A)*b Xn = >> norm(Xn) ans =

283 Całkowanie function wynik=oblicz(x) wynik=x.*log(x);>> quad 1,5) ans= >> quadl 1,5) ans=

284 quad (@oblicz, początek_x, koniec_x, tolerancja) funkcje quad, quadl Funkcja quad wykorzystuje regułę Simpson’a (3-punktowa reguła Newton-Cotes’a), dokładna dla wielomianów do 3 stopnia. Funkcja quadl 4-punktową regułę Gauss-Lobatto i 7-punktowe rozszerzeniem Kronroda. Dokładne (odpowiednio) dla wielomianów do 5 i 9 stopnia. Obie metody sa metodami adaptacyjnymi. Dzielą przedział całkowania na podprzedziały o zmiennej długości tak, aby najkrótszy przedział wypadał w miejscu największej zmienności funkcji podcałkowej. quad początek_x, koniec_x, tolerancja) quadl początek_x, koniec_x, tolerancja) tolerancja (jak duży bład) – opcjonalnie

285 tolerancja w całkowaniu>> 1,5,1e-1) ans = >> 1,5,1e-25) Warning: Maximum function count exceeded; singularity likely. (Type "warning off MATLAB:quadl:MaxFcnCount" to suppress this warning.) > In C:\MATLAB6P5\toolbox\matlab\funfun\quadl.m at line 94

286 Przykład zamiany nazw dla quadfunction wynik=calka(nazwa, poczatek_x, koniec_x) % całkuje dowolną funkcję zewnętrzną „nazwa” nazwa_fcji=str2num(sklej); wynik=quad(nazwa_fcji, poczatek_x, koniec_x); >> ile=calka('oblicz',1,5) ile = Inaczej to samo: >> calka >> calka('sin',1,2) ans =0.9564 >> quadl('sin',1,2) ans =

287 Całkowanie za pomocą „trapz”Z = trapz(Y) Z = trapz(X,Y) %Y-wektor wyznaczający przedział dziedziny Metoda NIEADAPTACYJNA >>x=linspace(1,5,5) >>f=x.*log(x); >>trapz(x,f)  ans = >>x=linspace(1,5,15) ans = >>X = 0:pi/100:pi; >>Y = sin(X); >>Z = trapz(X,Y) Z =

288 Całkowanie podwójne q = dblquad(fun,xmin,xmax,ymin,ymax)q = dblquad(fun,xmin,xmax,ymin,ymax,tol) q = dblquad(fun,xmin,xmax,ymin,ymax,tol,method) function out=fxy(x,y) out=y^2*exp(x)+x*cos(y);   >> ans=

289 Różniczkowanie

290 Funkcja diff diff(y) - oblicza różnicę pomiędzy kolejnymi elementami wektora diff(y,n) – różnica obliczana n razy Przykład: >> x=[1, 3, 5, 3, 4, 8] >> diff(x) ans = >> diff(x,2) ans =

291 Różniczkowanie diff(y)./diff(x) – aproksymacja różniczkowania Przykładx=[0:0.1:2*pi]; funkcja=sin(x); pochodna=diff(funkcja)./diff(x); plot(x, funkcja, x(2:end),pochodna,':') legend(' funkcja ', ' pochodna ')

292

293 Różniczkowanie na płaszczyźnie GRADIENTGradient to operator różniczkowy. Gradient pola skalarnego jest polem wektorowym pokazującym kierunek i wielkość zmian pola skalarnego. Zwrot wektora gradientu wskazuje kierunek, w którym pole skalarne rośnie najszybciej. Moduł wektora gradientu jest pochodną dla kierunku największego wzrostu [px,py]=gradient(f) [...] = gradient(F,h1,h2,...) gdzie h1, h2 – przesunięcie w odpowiednim kierunku

294 przykład 1 >> A=[1 3 24 5 78; 3 12 5 23 89; 8 1 3 12 6] A = >> [px,py]=gradient(A) px = py =

295 przykład 2 v = -2:0.2:2; [x,y] = meshgrid(v);z = x .* exp(-x.^2 - y.^2); [px,py] = gradient(z,.2,.2); contour(v,v,z) hold on quiver(v,v,px,py) hold off

296 dywergencja DIV = divergence(A) DIV = divergence(u,v,A)gdzie u, v – wektory przesunięcie w odpowiednim kierunku

297 przykład A = >> [px,py]=gradient(A); >> div=divergence(px,py) div =

298 przykład – obliczenie dywergencji z def.>> [px,py]=gradient(A) >> [pxx,pxy]=gradient(px) >> [pyx,pyy]=gradient(py) >> pxx+pyy ans = >> divergence(px,py) ans =

299 laplasjan Dyskretny laplasjan(U) można obliczyć za pomoca f-cji del2, która na płaszczyźnie oblicza 0.25*laplasjan(U). F-cja del2 wykonuje interpolację liniową wewnątrz przedziałów różniczkowania L = del2(U) L = del2(U,h) L = del2(U,hx,hy) L = del2(U,hx,hy,hz,...)

300 Przykład >> lap=4*del2(A) lap = 81 -1 -23 63 130 NIE ZGADZA SIĘ !! Porównajmy z wynikiem innej metody obliczenia, czyli div(grad(A)) >> [px,py]=gradient(A); div=divergence(px,py) div =

301 inny przykład div(grad(Z)) >> [X,Y]=meshgrid(0:1:6)>> [p1,p2]=gradient(Z) >> div=divergence(p1,p2) div =

302 >> div=divergence(p1,p2)Laplasjan obliczamy z del2 >> lap=4*del2(Z) lap = >> div=divergence(p1,p2) div = Porównanie zawodzi, ale ...

303 >> [X,Y]=meshgrid(0:1:6); >> delta_Z=12*X.^2+6*Y; W tym wypadku jednak znamy funkcję Z(X,Y) opisującą pole skalarne, możemy więc ręcznie obliczyć wynik i porównać która z metod daje lepszy wynik Z=X.^4+Y.^3 , czyli 12*X.^2+6*Y >> [X,Y]=meshgrid(0:1:6); >> delta_Z=12*X.^2+6*Y; >> delta_Z=12*X.^2+6*Y delta_Z =

304 Porównajmy błąd: >> err_lap=sqrt(sum((lap(:)-delta_Z(:)).^2)) div = delta_Z = Porównajmy błąd: >> err_lap=sqrt(sum((lap(:)-delta_Z(:)).^2)) err_lap = >> err_div=sqrt(sum((div(:)-delta_Z(:)).^2)) err_div =

305 Przykład zastosowania operatorów różniczkowych pola w Inżynierii Biomedycznej – modelowanie serca

306 Komórki mięśnia sercowego i ich elektryczna reprezentacjaprzewodnictwo omowe odpowiednik prawa Kirchoffa

307 Równania różniczkowe zwyczajne (ODE)Przykłady równania różniczkowego: Przesunięcie –prędkość, V(t) – jakaś określona funkcja

308 [t,y]=solver(‘odefile’, tspan,y0,options, P1, P2,…)ODE w Matlabie [t,y]=solver(‘odefile’, tspan,y0,options, P1, P2,…) [t,y]- wektor czasu i wartości rozwiazań y(t) solver – nazwa funkcji reprezentującej algorytm rozwiązania równań różniczkowych. Do pierwszych obliczeń zaleca się zastosowanie ode45. Algorytm ten daje poprawne rozwiązania dla dużej grupy równań. tspan – wektor, w którym określa się początkową i końcowa wartość zmiennej niezależnej (zwykle jest to czas). Można też podać wszystkie punkty (w kolejności rosnącej lub malejącej), dla których potrzebne są rozwiązania y(t). Jeżeli znane są warunki końcowe, a nie początkowe to można wykonać całkowanie wstecz (początkowa wartość czasu jest większa od wartości końcowej) y0 – wartości poczatkowe dla układu równań options – struktura, w której polach zapisuje się wartości parametrów całkowania równań różniczkowych. Wartości domyślne można także zmieniać później za pomocą funkcji odeset. P1, p2 – parametry dodatkowe, które można wprowadzić do M-pliku funkcyjnego odefile.

309 Przykład Aby rozwiązać to równanie należy stworzyć funkcję definiującą pochodną pierwszego rzędu function yprim=pochod1(t,y) yprim=-y-5*exp(-t)*sin(5*t);   skrypt rownanie1.m tspan=[0 3]; yzero=1; tspan, yzero); plot(t,y,’*--‘) xlabel t, ylabel y(t)

310 Sztywne równania różniczkoweSą to takie równania, w których przy źle dobranym kroku całkowania lokalne rozwiązanie numeryczne może się bardzo szybko zmieniać, chociaż prawdziwe rozwiązanie zmienia się dużo wolniej. W takich problemach bardzo istotny jest odpowiedni dobór kroku całkowania.

311 Problemy sztywne i wybór solvera Układ Robertsona: model reakcji chemicznej pomiędzy 3 reagentami:  Układ równań różniczkowych zapisujemy macierzowo: function yprim=chem(t,y) yprim=[-0.04*y(1)+1e4*y(2)*y(3); ... 0.04*y(1)+1e4*y(2)*y(3)- 3e7*y(2)^2; ... 3e7*y(2)^2];

312 skrypt robertson.m tspan=[0 3]; yzero=[1;0;0]; subplot(121), plot(ta,ya(:,2),'-*')  xlabel('t'), ylabel('y_2(t)'), title('ode45', 'FontSize',14) [tb,yb]= subplot(122), plot(tb,yb(:,2),'-*') xlabel('t'), ylabel('y_2(t)'), title('ode15', 'FontSize',14)

313

314 solvery ode45 – równania nie-sztywne, metoda Runge-Kutta expilcite, rzędu 4 i 5 ode23 - równania nie-sztywne, metoda Runge-Kutta expilcite, rzędu 3 i 4 ode113 - równania nie-sztywne, metoda liniowa wielokrotnego kroku, expilicite, rząd 1-13 ode15s- równania sztywne, metoda liniowa wielokrotnego kroku, impilicite, rząd 1-5 ode23s - równania sztywne, metoda Rosenbrocka, 1 kroku, expilicite, rząd 2 i 3 ode23t- równania umiarkowanie sztywne, reguła trapezu , implicite, rząd 2 i 3 ode23tb - równania sztywne, metoda Runge-Kutta impilcite, rzędu 2 i 3

315 Równanie różniczkowe drugiego rzędu - przykład:Podstawiamy: y1(t)=(t), y2(t)=d(t)/dt,  Można więc zapisać jako 2 równania pierwszego rzędu

316 [t,y]=ode45(@pend,tspan,yzero); Skrypt ode2rz.m tspan=[0 10]; yzero=[1; 1]; plot(t,y(:,1),'-black',t,y(:,2),':red') legend('y(t)','dy/dt(t)') function yprime= pend(t,y) yprime=[y(2);-sin(y(1))];

317 Parametry rozwiązania Można ustalać bezpośrednio lub za pomocą funkcji odeset options=odeset(‘pole1’,wartość1, ‘pole2’, wartość2,...) >>odeset AbsTol: [ positive scalar or vector {1e-6} ] – max błąd bezwzględny (domyślnie 1e-6) RelTol: [ positive scalar {1e-3} ] - błąd względny NormControl: [ on | {off} ] OutputFcn: [ function ] – nazwa f-cji wywoł. po każdym kroku całk. F-cja własna lub odeplot (rozwiazanie w f-cji czasu), odephas2(wykres faz. 2D lub 3D), odeprint (wydruk obliczonych wartosci) OutputSel: [ vector of integers ] Refine: [ positive integer ] Stats: [ on | {off} ] InitialStep: [ positive scalar ] – początkowawartość kroku całkowania MaxStep: [ positive scalar ] – max długość kroku całkowania (1/10 przedziału tspan) ...

318 Przykład zastosowania ODEPrzykład zastosowania ODE. Model neuronu Hodgkina-Huxleya (Nobel)

319

320

321

322

323 Rozwiązanie problemu metodą EuleraPrzechodzimy do przyrostów skończonych. (Poprawne tylko dla małych przyrostów !!)

324 % M. Turalska IB IV MODEL HH% Rozwiazywanie ukladu rownan rozniczkowych metoda Eulera % maksymalne przewodnictwo w mS/cm^2; 1=K, 2=Na, 3=leak g(1)=36; g(2)=120; g(3)=0.3; % wartosci potencjalow rownowagi dla poszczegolnych jonow Vi(1)=-12; Vi(2)=115; Vi(3)=10.613; % wartosci poczatkowe I_app=0; V=-10; mnh=zeros(1,3); licz=0; % krok dt=0.01; for t=0:dt:200 if t==10; I_app=10; end %wlaczenie zewnetrznego pradu if t==40; I_app=0; end %wylaczenie zewnetrznego pradu if t==90; I_app=8; end

325 % funkcje alfa i beta modelu 1-n, 2-m, 3-halpha(1)=(10-V)/(100*(exp((10-V)/10)-1)); alpha(2)=(25-V)/(10*(exp((25-V)/10)-1)); alpha(3)=0.07*exp(-V/20); beta(1)=0.125*exp(-V/80); beta(2)=4*exp(-V/18); beta(3)=1/(exp((30-V)/10)+1); % stala czasowa tau i wartosc rownowagowa mnh_inf tau=1./(alpha+beta); mnh_inf=alpha.*tau; % calkowanie ! mnh=(1-dt./tau).*mnh+dt./tau.*mnh_inf; % obliczenie przewodnictwa gnmh(1)=g(1)*mnh(1)^4; gnmh(2)=g(2)*mnh(2)^3*mnh(3); gnmh(3)=g(3);

326 %wewnetrzny prad jonowyI=gnmh.*(V-Vi); %potencjal membrany (całkowanie) ! V=V+dt*(I_app-sum(I)); if t > = 0; licz=licz+1; x_plot(licz) = t; y_plot(licz) = V; end plot(x_plot,y_plot); xlabel('Czas [ms]'); ylabel('Potencjal [mV]'); title('POTENCJAL CZYNNOSCIOWY');

327

328 Automatycznie (solver)% Malgorzata Turalska IB IV % MODEL HH - funkcja glowna % w zaleznosci od wartosci zewnetrznego pradu I_app % zalozonej w funkcji 'model' otrzymujemy pojedynczy % sygnal pobudzajacy (I_app=0) lub ciag impulsow (I_app=10) 100],[0.0;0.0;0.0;-10.0]); plot(t, y(:,4)); xlabel('Czas [ms]'); ylabel('Potencjal [mV]'); title('POTENCJAL CZYNNOSCIOWY');

329 function dxdydz=modelHH (x,y,z,k);%wartosci maksymalnego przewodnictwa dla poszczegolnych jonow g_K=36; g_Na=120; g_leak=0.3; % [mS/cm^2] %wartosci potencjalow rownowagi dla poszczegolnych jonow V_K=-12; V_Na=115; V_leak=10.613; %w mV I_app=10; % prad zewnetrzny % n = y(1), m = y(2), h = y(3), V = y(4), dzdydz=[n’, m’, h’, V’] dxdydz(1)=((10-y(4))/(100*(exp((10-y(4))/10)-1)))*(1-y(1))-(0.125*exp(-y(4)/80))*y(1); dxdydz(2)=((25-y(4))/(10*(exp((25-y(4))/10)-1)))*(1-y(2))-(4*exp(-y(4)/18))*y(2); dxdydz(3)=(0.07*exp(-y(4)/20))*(1-y(3))-(1/(exp((30-y(4))/10)+1))*y(3); dxdydz(4)=(-g_K*(y(1)^4)*(y(4)-V_K)-g_Na*(y(2)^3)*y(3)*(y(4)-V_Na) g_leak*(y(4)-V_leak)+I_app); dxdydz=dxdydz';

330

331 Operacje symboliczne Symbolic Math ToolboxW operacjach symbolicznych Matlab korzysta z jądra MAPLE - języka symbolicznego

332 Różniczkowanie symboliczne>> syms x f = sin(5*x) f = sin(5*x) >> diff(f) ans = 5*cos(5*x)

333 Różniczkowanie symboliczne – pochodne cząstkowe>> syms s t >> f = sin(s*t) >> diff(f,t) ans = cos(s*t)*s

334

335 Wybrane operacje symboliczneCałkowanie: R = int(S) R = int(S,v) Obliczanie granicy: limit(F,x,a) limit(F,x,a,'right') F-wyrażenie, x-zmienna, a- granica Suma szeregu: r = symsum(s) r = symsum(s,v) Rozwiązywanie równań g = solve(eq,var)

336 Przykłady >>syms x ; >> limit(sin(x)/x)>>limit(1/x,x,0,'left') >> syms k >> symsum(x^k/sym('k!'), k, 0,inf) %Uwaga: „zwykły” Matlab nie zna „4!” ( factorial(4)) ans =exp(x) >> solve('a*x^2 + b*x + c','b') ans= -(a*x^2+c)/x

337

338

339 Symboliczne rozwiązywanie równań różniczkowychr = dsolve('eq1,eq2,...', 'cond1,cond2,...', 'v') r = dsolve('eq1','eq2',...,'cond1','cond2',...,'v') >> syms y t >> y=dsolve('Dy=-y-5*exp(-t)*sin(5*t)','y(0)=1','t') y =exp(-t)*cos(5*t) >> diff(y,t) ans =-exp(-t)*cos(5*t)-5*exp(-t)*sin(5*t)

340 Równanie różniczkowe 2-rzęduy(0)=1, y’(0)=0 >> y = dsolve('D2y=cos(2*x)-y','y(0)=1','Dy(0)=0', 'x') simplify(y) y = 4/3*cos(x)-1/3*cos(2*x) ans = 4/3*cos(x)-2/3*cos(x)^2+1/3

341 Struktury – powtórzenie

342 Typ struct s = struct(‘pole1', wartosc1, ‘pole2', wartosc2, ...)s = struct(‘pole1', {}, ‘pole2', {}, ...)

343 Przykłady struktur >> s(1,1)=struct('imie','Anka','nazwisko',‘Kowalska','wzrost',165) s = imie: 'Ania‘ nazwisko: 'kowalska‘ wzrost: 165 >> s(2,1)=struct('imie',‘Ola','nazwisko',‘Nowak','wzrost',175) s = 2x1 struct array with fields: imie nazwisko wzrost

344 >> zz=cat(2,s(:).',s2(:).') zz = 1x4 struct array with fields: ans = imie: ‘Ola’ nazwisko: ‘Nowak‘ wzrost: 175 >> s(2,1).imie ans = Ola >> fieldnames(s) ans = 'imie‘ 'nazwisko‘ 'wzrost' >> s.imie ans = Ania ans = Ola >> zz=cat(2,s(:).',s2(:).') zz = 1x4 struct array with fields: imie nazwisko wzrost >> zz=cat(2,s(:)',s2(:)')

345 Programowanie ObiektoweZamiast ograniczać się do typów standardowych dla danego języka (np. w Matlabie: double, char, logical, cell, struct), możemy definiować własne typy zmiennych. Zwykle są to typy złożone (np. w Matlabie - określony rodzaj struktury). 2. Zmienne naszego nowego typu są powoływane do życia za pomocą jakiejś funkcji, napisanej przez nas, która ułatwi (powtarzający się) proces pojawiania się zmiennych takiego samego (nowego) typu. Taka funkcja nie tylko przyspiesza wielokrotne tworzenie zmiennych nowego typu, ale również umożliwia kontrolę poprawności tworzenia nowej zmiennej. W programowaniu obiektowym funkcja taka nazywa się KONSTRUKTOREM

346 3. W programowaniu obiektowym mówimy, że takie zmienne, należące do jednego typu i powołane do życia przez tą samą funkcję konstruktora, należą do jednej KLASY. 4. Na zmiennych należących do jednej klasy zwykle chcemy wykonywać podobne operacje. Tworzymy więc własne funkcje, które będą wykonywały te operacje. Funkcje te również należą do powołanej przez nas KLASY i nazywane są METODAMI. 5. Zakres ważności (widoczność) takiej funkcji ogranicza się do obiektów danej klasy. Można więc stworzyć kilka funkcji o tej samej nazwie, wykonujących podobne działanie, ale dla zmiennych innego typu (co może się wiązać np. ze szczegółami obliczeniowymi) i każdą z tych funkcji przypisać innej klasie. Użytkownik nie będzie sobie zdawał sprawy, że wywołując określoną nazwę funkcji nie zawsze wykonuje tę samą operację. Definiowanie RÓŻNYCH funkcji o tej samej nazwie to PRZEŁADOWANIE.

347 6. Klasy mogą być zorganizowane w strukturę hierarchiczną tak, aby obiekty dziedziczyły swoje metody. Klasy tworzą więc drzewo o zależnościach parent-child, w których dziecko dziedziczy metody klasy rodzica oraz jego pola. Ułatwia to tworzenie dużych projektów, przy których zatrudnione są zespoły programistów, umożliwiając korzystanie z metod wcześniej już zdefiniowanych. 7. W niektórych językach istnieją tez funkcje usuwające obiekty z pamięci. Funkcje te zwane są DESTRUKTORAMI. Matlabie nie istnieją destruktory, obiekty usuwane są standardowo za pomocą polecenia clear Obiekt lub clear classes.

348 klasy w Matlabie Klasę tworzymy w Matlabie poprzezUtworzenie katalogu o nazwie Katalog ten nie może się znajdować na ścieżce domyślnej Matlaba, zdefiniowanej przez zmienna „path” (w przeciwnym razie wszystkie funkcje danej klasy byłyby zawsze widoczne i niemozliwe byłoby przeładowanie funkcji). Zwykle katalog jest umieszczony wewnątrz katalogu danego toolboxa. Tworząc własną klasę w należy uaktualnić ścieżkę: addpath c:\myClasses;

349

350 2. W katalogu „@nazwa_klasy” umieszczamy co najmniej dwie m-funkcje:nazwa_klasy.m display.m Aby klasa miała sens należy też umieścić własne funkcje-metody. Funkcja nazwa_klasy.m jest konstruktorem. W tej funkcji, z otrzymanych argumentów wejściowych, tworzona jest zmienna danej klasy (w sposób określony przez projektanta) i deklarowana jako zmienna danego typu klasa. Obiekt jest jakimś typem strukturalnym. Struktura przechowuje stan obiektu. Funkcja display.m określa sposób wyświetlania zmiennej (np. A) z tej klasy w oknie Command Window, jeżeli zażąda tego użytkownik np. podjąc standardową komendę >> A

351

352 konstruktor function p = polynom(a)Przykład konstuktora klasy polynom, czyli function p = polynom(a) %POLYNOM – konstruktor klasy wielomianów, tworzy wielomian z % wektora a zawierającego współczynniki wielomianu if nargin == 0 p.c = []; p = class(p,'polynom'); else p.c = a(:).'; end

353 funkcja class obiekt = class(obiekt, 'nazwa_klasy')Funkcja ta definiuje zmienną jako obiekt wybranej klasy. Może ona wystąpić tylko w konstruktorze danej klasy Uzyskanie nazwy klasy: str = class(object) >> nazwa=class(wiel) nazwa = polynom

354 Funkcja isa isa – rozpoznaje czy obiekt obj istnieje w klasieK = isa(obj,‘nazwa_klasy’) polynom_obj = polynom([ ]); isa(polynom_obj, 'polynom') ans = 1

355 przykład konstruktora z isafunction p = polynom(a) %POLYNOM – konstruktor klasy wielomianów, tworzy wielomian z % wektora a zawierającego współczynniki wielomianu if nargin == 0 p.c = []; p = class(p,'polynom'); elseif isa(a,'polynom') p = a; else p.c = a(:).'; end

356 inputname przy okazji function y = moja_funkcja (a,b) y=a+b;zmienna1=inputname(1) zmienna2=inputname(2) >> x=7; moja_funkcja(x,5) zmienna1 = x zmienna2 = '' ans = 12

357 metoda display.m function display(p) % POLYNOM/DISPLAYdisp([inputname(1),' = ']) disp(['wielomian(',num2str(p.c(:)'),')'])

358 przykład utworzenia obiektu klasy polynomPowołanie obiektu >> wiel = polynom([ ]) argument f-cji display.m Brak średnika powoduje wykonanie metody display.m w klasie polynom, czyli wyświetlenie jego „wartości” >> wiel = polynom([ ]) wiel = wielomian( )

359 Antyprzykład 1 Co się stanie jeżeli nie zdefiniujemy funkcji display.m ? >> wiel = polynom([ ]) wiel = polynom object: 1-by-1 Matlab wyświetla informację o obiekcie tak jak umie. Tracimy możliwość eleganckiej obsługi naszych zmiennych

360 wiel = p = ... ??? Maximum recursion limit of 500 reached. Use set(0,'RecursionLimit',N) to change the limit. Be aware that exceeding your available stack space can crash MATLAB and/or your computer. Error in ==> polynom.display at 9 p Antyprzykład 2 function display(p) % POLYNOM/DISPLAY disp([inputname(1),' = ']) p >> wiel = polynom([ ])

361 Przykład metody, czyli konwersja typówfunction s = char(p) % POLYNOM/CHAR CHAR(p) łańcuch reprezentujący p.c if all(p.c == 0) s = '0'; else d = length(p.c) - 1; s = []; for a = p.c; if a ~= 0; if ~isempty(s) if a > 0

362 s = [s ' + ']; else s = [s ' - ']; a = -a; end if a ~= 1 | d == 0 s = [s num2str(a)]; if d > 0 s = [s '*']; if d >= 2 s = [s 'x^' int2str(d)]; elseif d == 1 s = [s 'x']; end d = d - 1;

363 Przykład display.m dla polynomfunction display(p) disp([inputname(1),' = ']) disp([' ' char(p) ]) >> p = polynom([ ]) p = x^3 - 2*x - 5

364 Przeładowanie operatorówWszystkie operatory Matlaba można przeładować funkcjami o predefiniowanych nazwach (obok). Wówczas operator (np. „+”) zmienia swoje znaczenie w danej klasie

365 Dziedziczenie Obiekt typu dziecko dziedziczy te same pola, które ma jego rodzic. Dziecko ma także własne pola różne od pól rodzica. Metody z klasy rodzica działają na obiektach z klasy dziecka. Metody z klasy dziecka nie są dostępne (niewidoczne) dla obiektów z klasy rodzica Deklaracja klasy wewnątrz konstruktora obiektu typu dziecko ma 3 argumenty (zamiast 2) childObj = class(childObj, 'childClass', parentObj); W momencie deklaracji klasy dziecka musi być już znany obiekt typu rodzic.

366 Wielo-dziedziczenie Dziedziczenie od więcej niż jednego rodzica:A) obj = class(s,'class_name',parent1,parent2...) B) obj = class(struct([]),'class_name',parent1,parent2... Wywołanie B) z pustą strukturą na początku gwarantuje, że dziecko nie ma żadnych własnych pól, które nie zostały odziedziczone.

367 Pierwszeństwo Klasy zdefiniowane przez użytkownika mają pierwszeństwo przed klasami wbudowanymi w Matlaba. Pierwszeństwem klas można sterować za pomocą funkcji: superiorto(‘inna_klasa1’,...) – nasz obiekt ma pierwszeństwo inferiorto(‘inna_klasa1’,...) – obiekty „innej klasy” mają priorytet przed naszą klasą Funkcje te muszą się pojawić w konstruktorze klasy !! Przykład: inferiorto(‘char’,’double’)

368 Powtórzmy Co to jest Programowanie Obiektowe ?Zamiast ograniczać się do typów standardowych dla danego języka (np. w Matlabie: double, char, logical, cell, struct), możemy definiować własne typy zmiennych. Zwykle są to typy złożone (np. w Matlabie - określony rodzaj struktury). 2. Zmienne naszego nowego typu są powoływane do życia za pomocą jakiejś funkcji, napisanej przez nas, która ułatwi (powtarzający się) proces pojawiania się zmiennych takiego samego (nowego) typu. Taka funkcja nie tylko przyspiesza wielokrotne tworzenie zmiennych nowego typu, ale również umożliwia kontrolę poprawności tworzenia nowej zmiennej. W programowaniu obiektowym funkcja taka nazywa się KONSTRUKTOREM

369 3. W programowaniu obiektowym mówimy, że takie zmienne, należące do jednego typu i powołane do życia przez tą samą funkcję konstruktora, należą do jednej KLASY. 4. Na zmiennych należących do jednej klasy zwykle chcemy wykonywać podobne operacje. Tworzymy więc własne funkcje, które będą wykonywały te operacje. Funkcje te również należą do powołanej przez nas KLASY i nazywane są METODAMI. 5. Zakres ważności (widoczność) takiej funkcji ogranicza się do obiektów danej klasy. Można więc stworzyć kilka funkcji o tej samej nazwie, wykonujących podobne działanie, ale dla zmiennych innego typu (co może się wiązać np. ze szczegółami obliczeniowymi) i każdą z tych funkcji przypisać innej klasie. Użytkownik nie będzie sobie zdawał sprawy, że wywołując określoną nazwę funkcji nie zawsze wykonuje tę samą operację. Definiowanie RÓŻNYCH funkcji o tej samej nazwie to PRZEŁADOWANIE.

370 6. Klasy mogą być zorganizowane w strukturę hierarchiczną tak, aby obiekty dziedziczyły swoje metody. Klasy tworzą więc drzewo o zależnościach parent-child, w których dziecko dziedziczy metody klasy rodzica oraz jego pola. Ułatwia to tworzenie dużych projektów, przy których zatrudnione są zespoły programistów, umożliwiając korzystanie z metod wcześniej już zdefiniowanych. 7. W niektórych językach istnieją tez funkcje usuwające obiekty z pamięci. Funkcje te zwane są DESTRUKTORAMI. Matlabie nie istnieją destruktory, obiekty usuwane są standardowo za pomocą polecenia clear Obiekt lub clear classes.

371 Funkcje Matlaba do pracy na obiektach

372 Współpraca Matlaba z innymi językami

373 Matlab i C++ Funkcja mex tworzy z programu program.c tzw. mex-plik, który można wywoływać z Matlaba. Przykłady na wykorzystanie funkcji mex można znaleźć w dysk:\matlabR12\extern\examples\mex Aby skorzystać z funkcji mex trzeba ją wcześniej zainstalować >> mex -setup

374 Instalacja funkcji mex – wybór kompilatora C>> mex -setup Please choose your compiler for building external interface (MEX) files: Would you like mex to locate installed compilers [y]/n? y Select a compiler: [1] Lcc C version 2.4 in E:\MATLAB701\sys\lcc [0] None Compiler: 1 Please verify your choices: Compiler: Lcc C 2.4 Location: E:\MATLAB701\sys\lcc Are these correct?([y]/n): y Try to update options file: C:\Documents and Settings\gosia\Application Data\MathWorks\MATLAB\R14\mexopts.bat From template: E:\MATLAB701\BIN\WIN32\mexopts\lccopts.bat Done . . .

375 Wywoływanie z Matlaba programu w C>>mex moja_funkcja.c W wyniku tego działania powstaje w bieżącym katalogu plik moja_funkcja.dll, który można wywoływać z Matlaba tak jak każdą inną funkcję Matlaba. Żeby to jednak prawidłowo zadziałało program moja_funkcja.c musi mieć specjalną konstrukcję: zawiera : #include mex.h  ( może też #include matrix.h  i in.) nagłówek program właściwy w C interfejs (mexFunction)

376 Przykład /* moja.c oblicza x^2-x+1/x dla kazdego elementu macierzy (przyklad w oparciu o Hanselmana). Wywolanie: p=moja(n) */ /*naglowek*/ #include mex.h /* to jest wlasciwa funkcja*/ static void moja(double p[], double n[], int r, int c) { int i; for (i=0; i p[i]=n[i]*n[i]-n[i]+1.0/n[i]; } ...

377 /* to jest interfejs do zmiennych w Matlabie */ void MexFunction( int nlhs, mxArray *plhs[], //wyjscie int nrhs, const mxArray *prhs[]) //argumenty-wejscie { double *p, *n, int r,c; //deklaracja wskaźników p,n do arg. we /* przykladowe sprawdzanie bledu wywolania funkcji przez uzytk. if (nlhs>1) mexErrMsgTxt(„Za duzo wyjsciowych wartosci”); /* pobranie rozmiaru macierzy wejsciowej */ r=mxGetM(prhs[0]); c=mxGetN(prhs[0]); /* tworzy macierz zmiennych wyjsciowych */ plhs[0]=mxCreateDoubleMatrix(r,c,mxREAL); /*tworzenie wskaznikow do arg. we */ p=mxGetPr(prhs[0]); n=mxGetPr(prhs[1]); /*wlasciwe wykonanie obliczen*/ moja(p,n,r,c);

378 Uwagi Plik nagłówkowy mex.h (są jeszcze inne) zawierają funkcje (jęz. C), których nazwy zaczynają się od mx... i mex... umożliwiają komunikację pomiedzy Matlabem a C Nazwa funkcji właściwej typu mex to nazwa m-pliku a nie funkcji wewnętrznej w C

379 Java Możliwe jest również bezpośrednie wywoływanie spod Matlaba funkcji Java (np. grafika Matlaba powstała w Javie) Przykład (uzyskanie danych sieciowych): function moj_IP try ja=java.net.InetAddress.getLocalHost; catch error('niemozliwy do ustalenia') end moja_nazwa=ja.getHostName moj_ip=ja.getHostAddress >> moj_IP moja_nazwa = user ea0 moj_ip =

380 I na zakończenie ... Czy program w Matlabie możemy skompilować?

381 Kompilacja programów mcc – Toolbox Compilermcc –m nazwa_pliku.m >> mcc -m RyR2sim

382

383