4. Achiziţia, interpolarea si aproximarea datelor

1. 4. Achiziţia, interpolarea si aproximarea datelor • MATLABul deţine o biblioteca de funcţii dedicate analizei datelor, incluzând si: -reprezentările grafice interactive (GUI); -metode de statistica descriptiva; -coeficienţi si metode de corelare, interpolare sau filtrare; -aproximarea unui set de date cu o functie; -cautarea datelor in tabele; -analiza Fourier si FFT.

2. 4.1. Importul si exportul fisierelor de date • MATLABul furnizează mai multe modalităţi de a importa sau exporta fişiere de date din / în fereastra de lucru (workspace); • Formatul datelor poate fi text, binar sau standard; • Datele text sunt în format ASCII, iar cele binare nu sunt in format ASCII şi nici nu pot fi văzute într-un editor de text; • Formatul binar reprezintă imagini, sunete sau alte informaţii; • Fişierele de date pot fi importate sau exportate şi prin utilizarea unor funcţii MATLAB dedicate. Pentru utilizarea acestor funcţii este necesar să se cunoascăformatul fişierului (txt, dat, csv, mat); • Cea mai simpla metodă, de a importa fişiere de date, utilizează editorulImport Wizard;

3. 4.2. Prelucrarea datelor si calcule statistice • Pentru a crea o matrice C care conţine elementele maxime sau minime ale matricelor de aceeaşi dimensiune A si B, se foloseşte sintaxa: C = max(A, B) sau C = min(A, B). • Daca fişierul de date conţine mai multe coloane, analiza acestora si calculul statistic sunt realizate separat (independent) pentru fiecare coloana. Asta înseamnă ca daca dorim sa calculam, de exemplu, maximul unei matrice, rezultatul va fi un vector linie care conţine valorile maxime ale datelor de pe fiecare coloana.

4. In tabelul următor se prezintă principalele funcţii MATLAB dedicate prelucrării datelor si calculelor statistice.

5. Prelucrarea datelor si calcule statistice. Exemplu • Datele eronate sau cele care nu îndeplinesc o anumita condiţie, impusa de utilizator, pot fi înlocuite, ca in exemplul următor: • Exemplu: Fie matricea A=[9 1000 4 8; 1 6 5 5; 3 2 7 1]; Sa se elimine coloanele matricei A care conţin cel puţin un element ce se abate de la valoarea medie cu mai mult de trei ori abaterea standard (legea celor 3 sigma din statistica). • Solutie: Se scrie un fisier - m cu secventa: m =mean(A); %valoarea medie sigma=std(A); %abaterea standard [n, p]=size(A); e =ones(n, 1); % generează o matrice coloana cu toate elem. egale cu 1 dist=abs(A–e*m); mcond=dist<3*e*sigma; %genereaza o matrice (3 x 4) cu toate elementele unu Y=A(:, find(all(mcond))) celim=p-length(find(all(mcond))) • Daca variabila celim=0, matricea Y va fi identica cu matricea A, deci nu se elimina nici o coloana. Daca Y e diferit de 0, atunci coloanele care nu indeplinesc conditia sunt eliminate din matricea A.

6. Derivarea si integrarea numerica a funcţiilor • Integrarea si derivarea sunt concepte fundamentale pentru rezolvarea unui număr mare de probleme atat in inginerie cat si in ştiinţa. • In unele situaţii nu pot fi obţinute soluţii analitice fiind necesara aplicarea metodelor de integrare si derivare numerica. • Derivata unei funcţii f(x) reprezinta viteza de variaţie a funcţiei in raport cu variabila x, notata cu dx: -Interpretarea geometrica a derivatei intr-un punct este panta tangentei la graficul funcţiei in punctul considerat; • Integrala funcţiei f(x) pe intervalul [a, b] are semnificaţia ariei delimitata de axa ox, curba f(x) si dreptele x = a si x = b

7. Calculul numeric al integralei • Calculul numeric al unei integrale, denumit si cuadratura, se poate face prin: • aproximarea funcţiei de integrat f(x) printr-o alta funcţie g(x) pe intervale finite; • aproximarea funcţiei f(x) cu un set de funcţii liniare pe porţiuni, aria calculându-se ca suma trapezelor care o compun (metoda trapezelor) • Daca aproximarea se face cu funcţii pătratice pe porţiuni, metoda de integrare se numeşte metoda lui Simpson.

8. Functiile MATLAB pentru integrare numerica • Functiile quad si quad8 se apeleaza cu sintaxele: quad(’f’, a, b), quad(’f’, a, b, tol), quad(’f’, a, b, tol, trace) • In care „f” este numele unui fisier functie care descrie functia f(x), a si b sunt limitele de integrare, tol – eroarea relativa admisa dintre 2 pasi consecutivi si trace – controleaza afisarea pe ecran a valorilor intermediare. • Functia trapz poate fi apelata cu una dintre cele doua sintaxe: Z = trapz(X, Y) sau Z = trapz(Y) • Cand functia trapz este apelata cu doua argumente, calculeaza integrala functiei y(x), X fiind un vector coloana care conţine abscisele iar Y este un vector coloana sau o matrice, cu acelasi numar de linii ca X. Z este un scalar daca Y este un vector, sau un vector linie daca Y este o matrice.

9. Integrare numerica. Exemple • Ex.: Sa se calculeze integrala functiei sin(x) pe intervalul [0, π] cu un pas de integrare de π / 100. • Solutie: Cu urmatoarea secventa se obtine rezultatul cerut: X = 0:pi/100:pi;Y = sin(X); Z = trapz(X, Y) • Se obtine rezultatul Z = 1.9998

10. Ex. : Utilizând regula lui Simpson si metoda trapezelor sa se calculeze integrala funcţiei y = sin(x) pe intervalul [0 - π]. • Solutie: Se scrie un fisier – m script cu secvenţa: • %Utilizarea regulii lui Simpson si metoda trapezelor pt. calculul integralei functiei y = f(x) a =0; b=pi; N=4;x = linspace(a, b, 2*N + 1);y = sin(x); for k=1:2*N+1 if k==1 | k==2*N+1 w(k)=1; elseif rem(k, 2)==0 w(k)=4; else w(k)=2; end end • Intsimp=((b-a)/(3*length(x)-1))*sum(y.*w) • Inttrapz=trapz(x,y) • Se obtine rezultatul: • Intsimp = 1.8464 • Inttrapz = 1.9742

11. Functia MATLAB pentru derivare • Functia MATLAB pentru derivare este diff si poate fi apelata cu sintaxele: Y = diff(X) sau Y = diff(X, n) sau Y = diff(X, n, dim) • daca X este un vector linie sau coloana, X = [x(1), x(2),..., x(n)], atunci functia diff(X) va returna vectorul diferentelor elementelor adiacente, adica XD = [x(2)-x(1), x(3)-x(2),..., x(n)-x(n-1)]; • daca X este o matrice atunci diferentele se calculeaza pe fiecare coloana prin scaderea elementului din linia imediat urmatoare (superioara): Unde m reprezinta nr. de linii.

12. Derivarea numerica. Exemple • Ex. 34: Sa se aproximeze derivata numerica a unei functii y(x) pentru X = [0, 1, 2, 3, 4, 5] si Y = [2, 6, 3, 5, 9, 12]. • Solutie: Functia diff permite aproximarea derivatei numerice a unei functii y(x) cu relatia: -Cu urmatoarea secventa MATLAB: X = [0, 1, 2, 3, 4, 5]; Y = [2, 6, 3, 5, 9, 12]; Z = diff(X); W = diff(Y); %dy = W = diff(Y) -Se returneaza rezultatele: Z = 1 1 1 1 1 W = 4 -3 2 4 3

13. Derivarea numerica. Exemple • Ex. 35: Sa se scrie un fisier function, utilizand instructiunea for, care sa realizeze aproximarea derivatei numerice a unei functii y(x) stiind ca X si Y sunt aceleasi ca in exemplul precedent. • Solutie: Se scrie un fisier function cu numele diffn.m cu urmatorul continut: function dY = diffn(X, n) dY=X; for i=1:n dY=diff(dY); end • Dupa rularea fisierului function diffn.m se tasteaza in fereastra de comanda: dy=diff(Y)./diff(X)sau dY • si se obtine acelasi rezultat ca in exemplul 34.

14. Derivarea si integrarea ecuatiilor diferentiale • O ecuaţie diferenţiala de ordinul 1 are forma: • Unde x este variabila independenta iar y este functia necunoscuta. • Un exemplu de ecuaţie diferenţiala de ordinul 1 este y’=3x2. Cunoscând condiţia iniţiala y(0)=-7.5 va rezulta soluţia y’ = x3-7.5. • Metodele numerice cele mai cunoscute, pentru rezolvarea ecuaţiilor diferenţiale, sunt metoda Euler si metoda Runge-Kutta (ode23, ode45). • Aceste functii se apeleaza cu sintaxa: [x, y] = ode23(’yprim’,x0, xf, y0, tol, trace) • In care x0 si xf sunt valorile iniţiala si finala a variabilei x, y0 este vectorul coloana conţinând condiţiile iniţiale, tol reprezintă toleranta iar trace parametrii care asigura tipărirea rezultatelor intermediare.

15. Integrarea ecuatiilor diferentiale. Exemplu • Ex. Sa se integreze ecuaţia diferenţiala y’ = 3x2 pe intervalul [2, 4] cu condiţia iniţiala y(2) = 0.5 si sa se reprezinte grafic funcţia rezultanta. • Soluţie: Se scrie un fişier function cu următorul conţinut: function dy=g1(x, y) dy=3*x^2; • Care se apelează cu următoarea secvenţa MATLAB, scrisa in fereastra de comanda sau, mai convenabil, intr-un fisier script: [x,yn]=ode23(’g1’, 2, 4, 0.5); %sol. numerica ya = x.^3 - 7.5; %solutia analitica plot(x, yn, ’*’, x, ya); grid

16. 4.3. Interpolarea si aproximarea datelor • Interpolarea unui set discret de date [xi, yi] presupune determinarea unei funcţii f(x) a.i. f(xi)=yi. • Sa presupunem ca dorim sa estimam valoarea y(x) pentru 2 puncte de coordonate (x1, y1) si (x2, y2), unde x1 < x < x2: • Daca punctele sunt unite printr-o dreapta, interpolarea se numeste liniara, • daca sunt unite printr-un polinom de gradul 3, interpolarea se numeste spline cubica. • O alta problema, frecvent întâlnita in inginerie electrica, consta in aproximarea unui set de date cu o funcţie care reprezintă „cea mai buna aproximare”. In acest caz funcţia determinata nu va trece prin toate punctele si va încerca sa elimine posibilele erori de măsurare. • De exemplu, metoda celor mai mici pătrate furnizează cea mai buna aproximare in sensul minimizării pătratului distantelor dintre punctele date si functia de aproximare.

17. Functiile MATLAB pentru interpolarea si aproximarea datelor.

18. Interpolarea functiilor de o singura variabila • Obiectivul interpolării il constituie estimarea valorilor unei funcţii f(x) pentru orice punct x=[x1, x2] al unui set de date. • Curba de interpolare trece prin toate punctele care o definesc. • Legea de interpolare poate fi liniara, cubica sau polinomiala. • Interpolarea spline cubica reprezintă o curba neteda dintre fiecare pereche de puncte, definita de un set de polinoame de gradul trei. De exemplu sase puncte sunt conectate intre ele prin cinci curbe diferite de gradul trei. Se apeleaza cu sintaxa: yi = spline(x, y, xi) unde x si y sunt vectorii care conţin abscisele si ordonatele setului de date (cu pas de eşantionare mare), xi este un vector care conţine noile abscise, de regula cu un pas mai fin, iar yi este vectorul returnat asociat vectorului xi. Obs: Toate legile de interpolare necesita ca elementele lui x sa fie ordonate crescator. In plus, metoda cubic cere ca punctele de pe axa x sa fie situate la distante egale.

19. Pentru a sublinia diferenţa dintre cele doua tipuri de interpolări (liniara si spline) se prezintă exemplul următor cu reprezentare grafica: • Ex: Sa se reprezinte grafic o interpolare liniara si una spline cubica prin 6 puncte. • Soluţie: Cu secvenţa: x =[0, 1, 2, 3, 4, 5]; y =[0, 20, 55, 60, 90, 100]; xi=0:0.1:5; yi=spline(x ,y, xi); plot(x, y, xi, yi, x, y,’o’) title(’Comparatie intre interpolare liniara si spline’); xlabel(’Timpul (s)’), ylabel(’Temperatura (grade)’);grid

20. Aproximarea datelor prin metoda celor mai mici pătrate • Pentru ca aproximarea sa fie considerata „cea mai buna” suma pătratelor distantelor de la fiecare punct la curba (dreapta) aproximata (linie sau polinom) trebuie sa fie minima. Este posibil ca nici un punct sa nu se găsească pe curba aproximata. Acest lucru separa foarte clar aproximarea de interpolare. • Aproximarea datelor printr-o linie dreapta se numeşte regresie liniara iar aproximarea unor date printr-un polinom se numeşte regresie polinomiala. • Regresia liniara minimizează suma pătratelor dintre dreapta de aproximare si punctele date. Măsura calităţii unei aproximări liniare este data de suma pătratelor, după cum reiese si din ecuaţia: • Determinarea parametrilor m si n ai dreptei de aproximare y=mx+n se face utilizând funcţia polyfit. • Determinarea celei mai bune aproximări a unui set de date (x, y) cu un polinom de ordinul n, poate fi apelata cu sintaxa:p = polyfit(x, y, n)

21. Aproximarea datelor prin metoda celor mai mici pătrate. Exemple. • Exemplu: Sa se aproximeze un set de date (x, y), utilizând metoda celor mai mici pătrate. • Soluţie:Se scrie un fişier MATLAB cu secvenţa următoare, obţinându-se graficul din figura: x = [0, 1, 2, 3, 4, 5]; y = [0, 20, 55, 60, 90, 100]; coef = polyfit(x, y, 1); %determinarea parametrilor m = coef(1); n = coef(2); %parametrii dreptei de aproximare y1 = m*x + n; sum_p = sum((y-y1).^2; %minimizarea sumei patratelor distantelor plot(x, y1, x, y, ’*’)

22. Aproximarea datelor prin metoda celor mai mici pătrate. Exemplu. • Funcţiapolyfitreturnează coeficienţii (ai) ai polinomului p(x), care in punctele precizate de vectorul x are, in sensul celor mai mici pătrate, valorile date de vectorul y. • Ex: Fie polinomul p(x) = x3-6x2+11x-6, peste care este suprapus un zgomot cu distribuţie normala. Aproximaţi datele rezultate cu un polinom de gradul 3 utilizând metoda celor mai mici pătrate. Reprezentaţi grafic datele cu zgomot si polinomul aproximat. • Solutie: Cu secvenţa următoare se obţine graficul din figura: p = [1, -6, 11, -6];x = 0 : .25 : 4; y = polyval(p, x)+randn(size(x)); %generarea polinomului p(x) peste care s-a %suprapus un zgomot cudistributie normala c = polyfit(x, y, 3); % aproximeaza un polinom de ord. 3 poli3 = polyval(c, x); % evaluarea polinomului c pt. valorile precizate de x plot(x, poli3, x, y,’*’);grid

23. Analiza Fourier si Transformata Fourier Rapida (FFT) • Analiza Fourier este extrem de utila pentru procesarea datelor, deoarece descompune un semnal intr-un sir de componente sinusoidale de frecvente diferite, făcând trecerea din domeniul timp in domeniul frecventa, realizând calculul amplitudinii si fazei variabilelor (datelor, semnalelor) transformate. • Pentru eşantionarea datelor vectoriale, analiza Fourier utilizează transformata Fourier discreta (discrete Fourier transform-DFT). • Transformata Fourier rapida - FFT (Fast Fourier Transform) este un algoritm foarte eficient pentru calcularea transformatei Fourier, sau a transformatei Fourier discrete (DFT). • De asemenea, FFT este o unealta utila pentru filtrarea datelor, procesarea semnalelor sau pentru procesarea imaginilor in domeniul frecventei si pentru estimarea spectrului puterii.

24. Analiza seriei Fourier • Daca g(t) este o functie periodica, cu perioada Tp: • Si daca in orice interval finit g(t) are cel puţin un număr de discontinuitati si un numar finit de minimuri si maximuri (condiţiile lui Dirichlet), atunci g(t) poate fi exprimata cu o serie de functii sinusoidale: • Unde w0 = 2*pi / T, a0 / 2 este componenta continua a seriei si reprezintă valoarea medie a funcţiei g(t) pe o perioada.

25. 4.4.3. Transformata Fourier • Unealta matematica pentru analiza unui semnal in domeniul frecventei este Transformata Fourier care poate lua diferite forme in funcţie de semnalul analizat. • Ceea ce au in comun aceste semnale este faptul ca sunt alcătuite dintr-un număr de componente sinusoidale de frecvente diferite, fiecare având o anumita amplitudine si faza initiala. • Transformata Fourier face conversia unui semnal din domeniul timp intr-un semnal discret in domeniul frecventa. • Daca g(t) este un semnal neperiodic exprimat ca funcţie de timp, transformata Fourier a funcţiei g(t) este data de expresia integrala: …..(18) • Valoarea semnalului g(t) poate fi obţinuta din expresia transformatei Fourier (ec. 18) utilizând relaţia transformatei Fourier inverse:

26. Transformata Fourier • Daca g(t) este un semnal continuu si periodic: • Unde Tp este perioada semnalului, va rezulta transformata Fourier a functiei g(t) sub forma: • iar cn se poate exprima prin relatia:

27. Transformata Fourier rapida (FFT) • Transformata Fourier rapida (FFT) este o metoda eficienta de calcul a transformatei Fourier discrete. • FFT reduce numărul de calcule matematice necesare pentru calculul transformatei Fourier discreta (DFT). • FFT poate fi utilizata pentru calculul spectrului puterii unui semnal, pentru filtrarea digitala a semnalelor sau pentru obţinerea corelaţiei dintre 2 semnale. • Functia MATLAB pentru calcularea transformatei Fourie discrete (DFT) a unui vector x, utilizand algoritmul de calcul al transformatei Fourier rapide (FFT) este:fft(x). • Functia fft(x, N) se utilizează pentru a obţine FFT pentru un vector de lungime N (cu N puncte). • Daca X este o matrice, functia fftreturnează transformata Fourier a fiecărei coloane. • Daca lungimea lui X este mai mica decat N, atunci matricea X este completata cu valoarea zero pana la dimensiunea N. Daca lungimea lui X este mai mare decât N, atunci secvenţa matricei X este trunchiata (redusa, micsorata).

28. MATLAB furnizează o colecţie de funcţii pentru calcularea si determinarea transformatei Fourier:

29. FFT - Exemple • Ex. : Sa se calculeze si sa se reprezinte grafic amplitudinea unui semnal x = sin(2π*15*t) utilizand functia fft. b) Sa se reprezinte pe acelasi grafic faza si amplitudinea semnalului. • Solutie: Se scrie un fisier script cu urmatorul continut, care va genera graficul din figura urmatoare: %Acest fisier calculeaza amplitudinea transformatei Fourier rapide a unei secvente t =(0:99) / 100; % vectorul timp x =sin(2*pi*15*t)+sin(2*pi*40*t); %semnalul achizitionat y =fft(x); % calculeaza DFT a semnalului m =abs(y); %amplitudinea semnalului f=(0:length(y)-1)'/length(y)*100; % vectorul frecventei subplot(211),plot(t,x),xlabel('timpul (s)'),ylabel('x(t)'), title(’Semnalul achizitionat in domeniul timp’) subplot(212),plot(f,m),xlabel('frecventa(Hz)'),ylabel('Amplitudinea semnalului'),title('Amplitudinea unui semnal calculata cu functia fft') p =unwrap(angle(y)); %functia angle calculeaza faza semnalului si functia unwrap elimina discontinuitatile fazei figure, subplot(211), plot(f,m), ylabel('Amplitudinea semnalului'), grid; subplot(212), plot(f,p*180 / pi), xlabel('Frecventa (Hz)'), ylabel('Faza (grade)'),grid

30. FFT – Reprez. graf. a ampl si fazei unui semnal

31. FFT - Exemple • Ex.: Sa se scrie un fisier function care sa calculeze transformata Fourier rapida a unui semnal de forma: x=2*sin(2*pi*50*t)+1.5*sin(2*pi*100*t)+sin(2*pi*200*t)+0.5*sin(2*pi*350*t);si care sa returneze vectorii frecventa si amplitudine. Sa se reprezinte pe acelaşi grafic amplitudinea si faza funcţie de frecventa. • Solutie: Se scrie următorul fişier function cu numele testfft: function [f,X]=testfft t =0:.001:0.39; %crearea vectorului timp fe=1 / (t(2)-t(1)); %frecventa de esantionare x=2*sin(2*pi*50*t)+1.5*sin(2*pi*100*t)+sin(2*pi*200*t)+0.5*sin(2*pi*350*t); Xt=fft(x); %calcularea DFT a secventei Xm=abs(Xt); %calcularea amplitudinii FFT N=length(x);X=Xm(1:N / 2+1);Xp=Xt(1:N / 2+1);P=unwrap(angle(Xp)); f =[0:N / 2]*fe / N; %crearea vectorului frecventa subplot(311), plot(t,x), subplot(312), plot(f,X), subplot(313), plot(f,P*180 / pi)

32. FFT - Exemple