P
Glushach V.S. UIT-44
Stăpânirea lucrării în MathCad. Dobândirea abilităților în utilizarea transformării Laplace pentru a analiza componentele spectrale ale semnalelor. Studiul scărilor de timp și frecvență ale seriilor temporale și transformării Fourier.
1. Generați o serie temporală de trei sinusoide. Numărul de puncte trebuie să fie 2^n
2. Determinați media și varianța.
3. Facem transformarea directă și inversă F. Suprapunem semnalul de două ori convertit pe graficul seriei temporale originale.
4. Găsiți relația dintre scara seriei de timp de-a lungul axei timpului și transformata Fourier de-a lungul axei frecvenței.
1. Selectați discreția temporală dt și numărul de puncte din seria temporală sub forma nl:= 2 k
Fie k:= 9 nl:= 2 k nl=512 Lungimea probei în timp T:=512
SH 
ag de Or, dat fiind că nl-1
timpul este aproximativ egal cu nl Atunci i:=0..nl-l t. := i*dt
2. Generam semnalul de intrare x ca suma a trei semnale armonice si determinam statisticile de baza ale acestuia.
А1:= 1 f1:= 0,05 xl i:= Al-sin/2*3,14*fl*t i) srl:= mean(xl) srl = 0,012 s1:=stdev(x1) s1=0,706
A2:= 0,5 f2:= 0,1 x2 i:= A2-sin/2*3,14*f2*t i) sr2:= medie(x2) sr2 = 3,792x10 -4 s2:=stdev(x2) s2=0,353
A3:= 0,25 f3:= 0,4 x3 i:= A3-sin/2*3,14*f3*t i) sr3:= medie(x3) sr3 = 3,362x10 -4 s3:=stdev(x3) s3=0,177
x i:= xl i + x2 i + x3 i sry:= medie(x) sry = 0,013 sy:= stdev(x) sy = 0,809
1. Transformată Fourier directă în MathCad F:= fft(x)
Perioada maximă a componentei armonice care poate fi într-o serie de timp este egală cu lungimea eșantionului. Această componentă armonică corespunde frecvenței minime posibile pe scara de frecvență a transformării Fourier frnin și, în consecință, pasului de-a lungul axei frecvenței transformării Fourier df.
Tmax:= Tfrnin:= 
df:= frnin df = 1,953 x 10 -3
Astfel, frecvența minimă și treapta de frecvență a transformării Fourier sunt egale cu frnin = df = 1/T.
Transformata Fourier are un număr de ordonate de frecvență care este jumătate mai mare decât numărul de ordonate ale seriei de timp în timp n2=nl/2 sau, inclusiv punctul zero (la care transformata Fourier nu este definită)
n2:= 1 + 2 k -1 n2 = 257 j:= l..n2
Indicele de frecvență curent se modifică de la j=l la j=n2
În acest caz, frecvența variază de la fmin =df= 1/T Frecvența maximă finax:= n2*df fmax = 0,502
pana la frnax=n2*df Frecventa curenta f i:= i*df
f 1 = 1,953 x 10 -3 f 257 = 0,502
DESPRE 
Vă rugăm să rețineți că transformata Fourier este definită numai pentru frecvențe în intervalul de la f=finin la f=fmax.
În acest caz, vârfurile de pe graficul spectrului Fourier corespund frecvențelor sinusoidelor originale, adică transformata Fourier vă permite să izolați componentele de frecvență ale semnalului. Dar amplitudinile componentelor armonice nu reflectă acum amplitudinile componentelor seriei temporale originale (unde A1=1, A2=0,5, A3=0,25)
Să remarcăm, de asemenea, că la dt = 1 frecvența maximă în spectrul transformării Fourier este frnax = 0,5 oscilații pe scară unității de timp. La dt = 1 sec aceasta corespunde cu fmax = 0,5 Hz. În acest caz, perioada de frecvență maximă este Tfmax=1/0,5=2. Aceasta înseamnă că pentru o perioadă de frecvență maximă există două selecții ale seriei temporale. Aceasta corespunde teoremei lui Kotelnikov, conform căreia pentru a restabili un semnal armonic continuu de la unul discret fără pierderi de informații pentru o perioadă trebuie să existe cel puțin două mostre în timp.
3. Să verificăm coincidența seriilor de timp înainte și după transformarea dublă Fourier. Pentru a face acest lucru, obținem transformata Fourier inversă din transformarea directă rezultată. Ar trebui să coincidă cu seria temporală inițială, care este confirmată de următorul grafic FF:= ifft(F)
Semnificația matematică a transformării Fourier este de a reprezenta semnalul y(x) ca o sumă infinită de sinusoide de forma F(v)sin(vx). Funcția F(v) se numește transformată Fourier sau integrală Fourier sau spectrul Fourier al semnalului. Argumentul său v are semnificația frecvenței componentei semnalului corespunzătoare. Transformarea Fourier inversă transformă spectrul F(V) în semnalul original y(x). Conform definiției,
După cum se poate observa, transformata Fourier este o cantitate esențial complexă, chiar dacă semnalul este real.
Transformată Fourier a datelor reale
Transformarea Fourier este de mare importanță pentru diverse aplicații matematice și pentru aceasta a fost dezvoltat un algoritm foarte eficient numit FFT (Fast Fourier Transform). Acest algoritm este implementat în mai multe funcții Mathcad încorporate care diferă în normalizare.
- fft(y) - vector direct de transformare Fourier;
- FFT(Y) este vectorul transformării directe Fourier într-o normalizare diferită;
- ifft(v) - vector de transformare Fourier invers;
- IFFT(V) este vectorul transformării Fourier inverse într-o normalizare diferită;
- y este un vector de date reale luate la intervale egale de valori ale argumentului;
- v este un vector de date reale din spectrul Fourier luate la intervale egale de valori ale frecvenței.
Argumentul transformării directe Fourier, adică vectorul y, trebuie să aibă exact 2 n elemente (n este un număr întreg). Rezultatul este un vector cu 1+2 n-1 elemente. Invers, argumentul transformării Fourier inversă trebuie să aibă 1+2 n-1 elemente, iar rezultatul său va fi un vector de 2 n elemente. Dacă numărul de date nu coincide cu puterea lui 2, atunci este necesar să completați elementele lipsă cu zerouri.
Orez. 15.24. Date sursă și transformată Fourier inversă (Listarea 15.20)
Un exemplu de calcul al spectrului Fourier pentru suma a trei semnale sinusoidale de amplitudini diferite (prezentat ca o curbă solidă în Fig. 15.24) este dat în Lista 15.20. Calculul se efectuează folosind N=128 de puncte și se presupune că intervalul de eșantionare a datelor хх este egal cu A. În penultima linie a listării, se utilizează funcția încorporată dacă ft, iar în ultima linie valorile corespunzătoare ale frecvențelor Qx sunt corect determinate. Vă rugăm să rețineți că rezultatele calculului sunt prezentate sub forma modulului spectrului Fourier (Fig. 15.25), deoarece spectrul în sine este complex. Este foarte util să comparați amplitudinile obținute și locația vârfurilor spectrului cu definiția sinusoidelor din Lista 15.20.
Listare 15.20. Transformare rapidă Fourier
Orez. 15.25. Transformată Fourier (Listarea 15.20)
Rezultatul transformării Fourier inversă este prezentat sub formă de cercuri în aceeași figură. 15.24, la fel ca datele originale. Se poate observa că, în cazul în cauză, semnalul y(x) este restabilit cu mare precizie, ceea ce este tipic pentru o schimbare lină a semnalului.
Transformarea Fourier a datelor complexe
Algoritmul rapid de transformare Fourier pentru date complexe este încorporat în funcțiile corespunzătoare, al căror nume include litera „c”.
- cfft(y) - vector de transformată Fourier complexă directă;
- CFFT(y) este vectorul transformării Fourier complexe directe într-o normalizare diferită;
- icfft(y) este vectorul transformării Fourier complexe inverse;
- ICFFT(V) este vectorul transformării Fourier complexe inverse într-o normalizare diferită;
- y este un vector de date luate la intervale egale de valori ale argumentului;
- v este un vector de date din spectrul Fourier luate la intervale egale de valori ale frecvenței.
Funcțiile de transformare Fourier reală profită de faptul că, în cazul datelor reale, spectrul este simetric în jurul zero și scoate doar jumătate din el (vezi secțiunea „Transformarea Fourier a datelor reale” de mai sus în acest capitol). Prin urmare, în special, folosind 128 de date reale, au fost obținute doar 65 de puncte din spectrul Fourier. Dacă aplicăm funcția complexă de transformare Fourier la aceleași date (Fig. 15.26), obținem un vector de 128 de elemente. Comparând Fig. 15.25 și 15.26, este posibil să înțelegem corespondența dintre rezultatele transformării Fourier reale și complexe.
Orez. 15.26. Transformată Fourier complexă (continuare din Lista 15.20)
Transformată Fourier bidimensională
Mathcad are capacitatea de a aplica funcții încorporate ale transformării Fourier complexe nu numai matricelor unidimensionale, ci și bidimensionale, adică matrici. Un exemplu în acest sens este prezentat în Lista 15.21 și Fig. 15.27 sub forma unui grafic al liniilor de nivel ale datelor inițiale și spectrul Fourier calculat.
Lista 15.21. Transformată Fourier bidimensională
Orez. 15.27. Date (stânga) și spectrul său Fourier (dreapta) (Listarea 15.21)
Partea 3. Rezolvarea ecuațiilor diferențiale ordinare în Mathcad
Seria Fourier pe un segment arbitrar
Partea 2. Extinderea funcțiilor în seria Fourier
Acțiuni cu numere complexe
Partea 1. Calcule cu numere complexe în Mathcad
Prelegerea nr. 5
Subiect: « Variabile complexe. Extinderea funcțiilor în seria Fourier. Rezolvarea ecuațiilor diferențiale»
În Mathcad, unitatea imaginară i este definită și, prin urmare, sunt definite numere complexe și operații cu acestea.
Z=a+bi– forma algebrică de scriere a unui număr complex.
a – parte reală, b – parte imaginară
Forma exponențială (exponențială) de scriere a unui număr complex,
A – modul, φ – argument (fază)
Forma trigonometrică de scriere a unui număr complex.
Relația dintre mărimi: a=A cos φ b=A sin φ
Z1=a1+j b1, Z2=a2+j b2
a) Adunare (scădere) Z3=Z1±Z2=(a1±a2)+j·(b1±b2)
b) Înmulțirea c·Z1=a·c+j·b·c
Z3=Z1·Z2=(a1·a2-b1·b2)+j·(a1·b2+a2·b1)=A1A2ej(φ1+φ2)
c) Diviziune
d) Ridicarea la puterea n (naturala)
e) Extracția rădăcinii: , unde k =0,1,2…n-1
Aparatul accepta doar radiani!!! radian=grad grad=radian
Exemple:
Funcția f(x) este absolut integrabilă pe intervalul [-p;p] dacă există o integrală. Fiecare funcție absolut integrabilă f(x) pe intervalul [-p;p] poate fi asociată cu seria sa trigonometrică Fourier:
Coeficienții seriei trigonometrice Fourier se numesc coeficienți Fourier și se calculează folosind formulele Euler–Fourier: ,
Să notăm a n-a sumă parțială a seriei Fourier a unei funcții netede pe bucăți f(x) pe intervalul [-p;p]. Abaterea standard este determinată de formula:
Pentru orice funcție mărginită f(x) integrabilă pe [-p;p], suma parțială a seriei sale Fourier este un polinom trigonometric cu cea mai bună aproximare a gradului al n-lea.
Exemplu:
Graficele arată cum converg sumele parțiale ale seriei Fourier. În vecinătatea punctelor de continuitate ale funcției f(x), diferența dintre valoarea funcției în punctul x și valoarea sumei parțiale a seriei în acest punct tinde spre zero ca n®¥, care este pe deplin în concordanţă cu teoria, întrucât în acest caz. Se poate observa, de asemenea, că diferența tinde să se zero cu atât mai rapid cu cât punctul x este situat mai departe de punctele de discontinuitate ale funcției.
Exemplu:
Pentru o funcție netedă pe bucăți pe intervalul [-L;L] al funcției f(x), problema extinderii seriei Fourier pe intervalul [-L;L] prin înlocuire liniară se reduce la problema extinderii funcției. pe intervalul [-p;p]:
Să luăm în considerare simplificările din seria Fourier în diferite condiții de simetrie:
formula (1) formula (2)
Să fie necesar să se găsească o soluție la ecuație
cu condiția inițială. Această sarcină se numește Problema Cauchy . Să extindem funcția dorită într-o serie în apropierea punctului și să ne limităm la primii doi termeni ai expansiunii. Luând în considerare ecuația (1) și notând-o, obținem Această formulă poate fi aplicată de mai multe ori, găsind valorile funcției în tot mai multe puncte noi.
Această metodă de rezolvare a ecuațiilor diferențiale obișnuite se numește metoda lui Euler . Geometric, metoda lui Euler înseamnă că la fiecare pas aproximăm soluția (curba integrală) printr-un segment tangent desenat la graficul soluției la începutul intervalului. Precizia metodei este scăzută și este de ordinul a h. Ei spun că metoda lui Euler este o metodă de ordinul întâi, adică precizia ei crește liniar cu pasul descrescător. h.
Există diverse modificări ale metodei lui Euler pentru a crește acuratețea acesteia. Toate se bazează pe faptul că derivata calculată la începutul intervalului este înlocuită cu valoarea medie a derivatei pe acest interval.
Deoarece rezultatul formulelor de interpolare ale lui Newton și Lagrange este același polinom de ordinul N, erorile lor se comportă în același mod.
Exemplul 3.4. Pentru datele inițiale utilizate în exemplul 3.1, calculăm valoarea polinomului lui Newton. Mai întâi, completați tabelul diferențelor împărțite:
F(xi,xj) |
F(xi ,xj ,xk ) |
F(x0 ,x1 ,x2 ,x3 ) |
z–xi |
||
Folosind formula lui Newton, obținem:
P 3 (1)= –1+0,6 1+(–0,1) 1 (–1)+0,0857 1 (–1) (–2)= –0,129.
3.6 Seria Fourier
Seria Fourier ne permite să studiem atât funcțiile periodice, cât și neperiodice prin descompunerea lor în componente. Curenții și tensiunile alternative, deplasările, viteza și accelerația mecanismelor de manivelă, undele acustice sunt exemple practice tipice de utilizare a funcțiilor periodice în calculele inginerești. În termeni de procesare a semnalului, transformata Fourier ia reprezentarea în serie de timp a unei funcții de semnal și o mapează la un spectru de frecvență. Adică transformă o funcție a timpului într-o funcție a frecvenței; Aceasta este descompunerea unei funcții în componente armonice la frecvențe diferite. Transformarea Fourier poate reprezenta un semnal care variază în timp în funcție de frecvență și amplitudine și oferă, de asemenea, informații despre fază (Fig. 3.4).
Expansiunea seriei Fourier se bazează pe presupunerea că toate funcțiile de importanță practică din intervalul π ≤x≤ π pot fi exprimate sub formă de serii trigonometrice convergente (o serie este considerată convergentă dacă șirul sumelor parțiale compuse din termenii săi converge). ).
Conform ipotezei Fourier, nu există nicio funcție care să nu poată fi extinsă într-o serie trigonometrică. Să extindem funcția f (t) într-o serie pe intervalul [–π, π]
f (t) = a 0 /2 + a 1 cos(t) + a 2 cos(2t) + a 3 cos(3t) + … + b 1 sin(t) + b 2 sin(2t) + b 3 sin (3t)+…,
unde cele n-ele elemente ale seriei sunt exprimate ca
f (t) cos(nt)dt , |
|||||
Orez. 3.4. Ilustrație pentru extinderea seriei Fourier
Se numesc coeficienții a n și b n Coeficienții Fourier, iar reprezentarea funcției f (t) conform formulei (3.1) este Extinderea seriei Fourier. Uneori, o expansiune serie Fourier prezentată în această formă se numește o expansiune serie Fourier reală, iar coeficienții sunt numiți coeficienți Fourier reali (spre deosebire de o expansiune complexă).
Să analizăm expresiile (3.2) și (3.3). Coeficientul a 0 reprezintă valoarea medie a funcției f (t) pe segmentul [–π, π] sau componenta constantă a semnalului f (t). Coeficienții a n și b n (pentru n > 0) sunt amplitudinile componentelor cosinus și sinus ale funcției (semnal) f (t) cu o frecvență unghiulară egală cu n. Cu alte cuvinte, acești coeficienți specifică mărimea componentelor de frecvență ale semnalelor. De exemplu, atunci când vorbim despre un semnal sonor cu frecvențe joase (de exemplu, sunetele unei chitare bas), aceasta înseamnă că coeficienții a n și b n sunt mai mari pentru valori mai mici ale lui n și invers - în frecvență înaltă sunet
vibrațiile (de exemplu, sunetul unei viori) sunt mai mari pentru valori mari de n.
Oscilația celei mai lungi perioade (sau cea mai joasă frecvență), reprezentată de suma a 1 cos(t) și b 1 sin(t) se numește oscilația frecvenței fundamentale sau a primei armonice. O oscilație cu o perioadă egală cu jumătate din perioada frecvenței fundamentale este a doua armonică, o oscilație cu o perioadă egală cu 1/n din frecvența fundamentală este n-armonica. Astfel, prin extinderea funcției f (t) într-o serie Fourier, putem face tranziția din domeniul timpului în domeniul frecvenței. Această tranziție este de obicei necesară pentru a identifica caracteristicile semnalului care sunt „invizibile” în domeniul timpului.
Vă rugăm să rețineți că formulele (3.2) și (3.3) sunt aplicabile pentru un semnal periodic cu o perioadă egală cu 2π. În cazul general, un semnal periodic cu perioada T poate fi extins într-o serie Fourier, apoi segmentul [–T /2, T /2] este utilizat în expansiune. Perioada primei armonice este egală cu T și componentele iau forma cos(2πt /T) și sin(2πt /T), componentele armonicii n sunt cos(2πtn /T) și sin(2πtn /T). ). Dacă notăm frecvența unghiulară a primei armonice ω0 = 2π/T, atunci componentele n-armonicilor iau forma cos(ω0 nt), sin(ω0 nt) și
cos(nt) b sin(nt) , |
||||||||||||||
f(t) |
||||||||||||||
unde coeficienții Fourier sunt calculați folosind formulele |
||||||||||||||
T/2 |
(t ) cos(0 nt )dt , |
T/2 |
f (t ) sin(0 nt )dt . |
|||||||||||
T/2 |
T/2 |
|||||||||||||
Expansiunea seriei Fourier este utilizată pentru analiza armonică sau spectrală a semnalelor periodice. Pentru analiza spectrală a semnalelor neperiodice se utilizează transformata Fourier. Pentru a face acest lucru, reprezentăm seria (3.4) folosind un sistem de funcții de bază sub formă de exponențiale cu exponenți imaginari:
2 jnt |
||||||
f(t) |
C n exp( |
|||||
T/2 |
2 jnt |
|||||
f(t)exp( |
||||||
T/2 |
||||||
Omitând un număr de calcule, scriem expresia (3.6) sub forma
C () f (t ) exp(j t )dt .
Această formulă se numește transformată Fourier directă sau transformata Fourier. De obicei, transformata Fourier este notată cu aceeași literă (numai majuscule) ca și funcția care este aproximată (care este de obicei notă cu o literă mică -
F () f (t ) exp(j t )dt .
Funcția F (ω) se numește funcție densitatea spectrală(sau pur și simplu densitate spectrală, transformată Fourier, imagine Fourier). Gama de valori ale funcției F (ω) în cazul general este mulțimea de numere complexe.
Transformată Fourier inversă , oferind restaurare
Modificarea funcției inițiale f (t) față de funcția de densitate spectrală se calculează după cum urmează
f (t ) F () exp(j t )dt .
Transformată Fourier discretă (DFT, DFT - Discrete Fourier Transform) este una dintre transformările Fourier utilizate pe scară largă în algoritmii de procesare a semnalului digital (modificările sale sunt utilizate în compresia audio în MP3, compresie imagini în JPEG etc.), precum și în alte domenii legate de analiza frecvențelor într-un semnal discret (de exemplu, analog digitizat). Transformata Fourier discretă necesită o funcție discretă ca intrare. Astfel de funcții sunt adesea create prin eșantionare (valori de eșantionare din funcții continue). Dezavantajul acestui algoritm este cantitatea mare de calcule repetitive. Eliminarea acestor operațiuni redundante are ca rezultat așa-numitul algoritm
transformată Fourier rapidă, care este de obicei folosită.
Transformare rapidă Fourier (FFT, FFT) - un algoritm pentru calcularea rapidă a transformatei Fourier discrete (DFT). Adică, algoritmul de calcul ia un număr de acțiuni mai mic decât O(N 2 ), necesar pentru calculul direct (prin formulă) a DFT (N este numărul de valori ale semnalului măsurate pe o perioadă, precum și numărul de componente de descompunere). Uneori, FFT înseamnă unul din algoritmi rapizi, numiți algoritmul decimare frecvență/timp sau algoritmul radix 2.
Pentru a implementa transformarea Fourier în pachetul MathCAD, trebuie să selectați operatorul Fourier în panoul simbolic pentru transformarea directă și invfourier pentru cea inversă. Acest operator trebuie plasat lângă funcția care trebuie convertită, iar singurul parametru trebuie să fie variabila relativ la care va fi convertită această funcție. Exemple de utilizare a afișajului-
noi poza. 3.5 pentru funcția f (t) e 2 t și în Fig. 3.6, unde modulația amplitudine-frecvență este aplicată funcției f (t), iar apoi rezultatul este extins într-o serie.
Orez. 3.5. Exemplu de extindere a seriei Fourier folosind funcția simbolică Fourier
Orez. 3.6. Exemplu de extindere a seriei Fourier folosind funcția simbolică Fourier
MathCAD conține funcții pentru transformarea Fourier rapidă discretă (FFT) și inversarea acesteia. Există două tipuri de funcții pentru transformata Fourier discretă: fft și ifft, cfft și icfft. Aceste funcții sunt discrete: iau vectori și matrice ca argumente și le returnează.
Funcțiile fft și ifft sunt utilizate dacă sunt îndeplinite următoarele condiții: (1) argumentele sunt reale; (2) – vectorul de date are 2m elemente.
În toate celelalte cazuri, sunt utilizate funcțiile cfft și icfft.
Prima condiție trebuie îndeplinită deoarece funcțiile fft și ifft folosesc faptul că pentru date reale a doua jumătate a transformării Fourier este complex conjugată la prima. MathCAD renunță la a doua jumătate a vectorului rezultat, ceea ce economisește timp și memorie în timpul calculelor. Perechea de funcții cfft și icfft nu utilizează simetria în conversie și poate fi folosită pentru numere reale și complexe.
A doua condiție este necesară deoarece perechea de funcții fft și ifft utilizează un algoritm de transformare Fourier rapidă foarte eficient. Pentru acest vector de argument, utilizați-
dat de funcția fft , trebuie să fie format din 2m elemente. Algoritmul pentru funcțiile cfft și icfft acceptă ca argumente vectori și matrice de dimensiuni arbitrare. Pentru transformata Fourier 2D, sunt utilizate numai aceste funcții. Funcțiile fft și ifft, cfft și icfft sunt reciproc inverse, adică adevărate:
și icfft(cfft(v)) v . |
În fig. Figura 3.7 ilustrează utilizarea funcțiilor ff t(v) și ifft(v) pentru un semnal sinusoid care este interferat cu utilizarea funcției rnd(x), care generează numere aleatorii în intervalul de la 0 la x.
Orez. 3.7. Transformată Fourier directă și inversă folosind funcțiile fft și ifft
Aceste grafice arată imaginea Fourier a semnalului c și o comparație a semnalului original x cu cel reconstruit din imaginea Fourier. Mai multe detalii despre analiza Fourier pot fi citite în și.
3.7 Metoda celor mai mici pătrate
În toate metodele de mai sus pentru aproximarea unei funcții, condițiile de interpolare au fost îndeplinite exact. Totuși, în cazurile în care datele inițiale x i, f i, i= 1,...,N, sunt date cu o oarecare eroare, poate fi necesară doar o execuție aproximativă.
Definiția condițiilor de interpolare: |F(x i ) – f i |< . Это условие означает, что интерполирующая функция F(x) проходит не точно через заданные точки, а в некоторой их окрестности, так, например, как это показано на рис. 3.8. Приблизим исходные данные глобальным полиномом. Если решать задачу интерполяции точно, то полином должен иметь степень N . При рассмотрении полинома Лагранжа мы выяснили, что полином N –й степени хорошо приближает исходную функцию только при небольших значениях N .
Orez. 3.8. Îndeplinirea aproximativă a condițiilor de interpolare
Vom căuta un polinom de grad scăzut, de exemplu, P 3 (x)=a 1 +a 2 x+a 3 x 2 +a 4 x 3. Dacă N >4, atunci problema exactă nu are soluții: pentru patru coeficienți necunoscuți (a 1 , a 2 , a 3 , a 4 ), condițiile de interpolare dau N > 4 ecuații. Dar acum nu este necesară îndeplinirea exactă a condițiilor de interpolare; dorim ca polinomul să treacă în apropierea punctelor date. Există multe astfel de polinoame, fiecare dintre acestea fiind definit de propriul său set de coeficienți. Dintre toate polinoamele posibile de acest tip, îl alegem pe cel care are cea mai mică abatere standard la nodurile de interpolare din valorile date, adică. Polinomul trebuie să fie cel mai apropiat de punctele date dintre toate polinoamele posibile de gradul trei în sens metoda celor mai mici pătrate(MNC). La punctul i-lea
linia P 3 (x) se abate de la valoarea f i cu suma (P 3 (x i) – f i). Însumăm abaterile pătrate ale polinomului pe toate punctele i= 1, 2,…, N și obținem funcționalitatea abaterilor pătrate:
G(a1 ,a2 ,a3 ,a4 ) (P3 (xi ) fi )2 |
a2 xi a3 xi 2 a4 xi 3 fi )2 . |
|
Să găsim minimul acestei funcționalități. Pentru a face acest lucru, echivalăm cu zero derivatele sale parțiale în raport cu variabilele a 1, a 2, a 3, a 4. Folosind regulile standard de diferențiere, obținem:
2 (a 1 |
a2 xi a3 xi 2 a4 xi 3 fi ) 0 |
||||||
a3 xi 2 a4 xi 3 fi ) 0 |
|||||||
G 2 xi (a1 a2 xi |
|||||||
2 x i 2 (a 1 |
a2 xi |
a3 xi 2 a4 xi 3 fi ) 0 |
|||||
2 x i 3 (a 1 |
a2 xi |
a3 xi 2 a4 xi 3 fi ) 0 |
|||||
Colectând coeficienții pentru necunoscutele a i, obținem un SLAE față de vectorul necunoscutelor (a 1, a 2, a 3, a 4):
N a1 xi a2 xi 2 |
a3 xi 3 a4 fi |
|||||||
xi 2 a2 |
xi 3 a3 |
xi 4 |
f i x i |
|||||
xi 2 a1 |
xi 3 a2 |
xi 4 a3 |
xi 5 |
f i x i2 |
||||
xi 3 a1 |
xi 4 a2 |
xi 5 a3 |
xi 6 |
f i x i3 |
||||
Sistemul rezultat se numește normal. Pentru a o rezolva, se folosesc metode standard de rezolvare a SLAE-urilor. De regulă, numărul de necunoscute ale sistemului (adică numărul de coeficienți ai funcției de interpolare) este mic, așa că puteți utiliza metode exacte pentru rezolvarea SLAE-urilor, de exemplu, metoda Cramer sau metoda Gauss. Metoda celor mai mici pătrate vă permite să „aproximați” datele originale folosind o combinație liniară a oricăror funcții elementare. Aproximațiile liniarului F (x )=a 1 +a 2 x, , trigonometric F (x )=a 1 sin(x )+a 2 cos(x), exponențial F (x )=a 1 e x sunt adesea folosite
N a1 xi a2
xi a1 xi 2 a2 fi xi
Noi calculăm
xi 2,
f i x i,
înlocuiți-l în normal
Orez. 3.9. Selectarea liniară
5a 1.4a
Dependențe MNC
0,148. Graficul funcției F (x)=-0,04+0,57x este prezentat în Fig. 3,9 cu o linie continuă. Punctele arată datele originale. Se poate observa că funcția liniară găsită aproximează de fapt punctele date.
În MathCAD, metoda celor mai mici pătrate este strâns legată de regresia liniară (y(x) = b + ax), deoarece coeficienții a și b sunt calculați din condiția minimizării sumei erorilor pătrate |b + ax i – y i | . Pentru calcule în MathCAD există două metode care se suprapun:
linia (x,y) returnează un vector cu două elemente de coeficienți de regresie liniară b + ax ;
Mod(x, y) – restul împărțirii lui x la y. Rezultatul are același semn ca x; unghi(x, y) – unghi (în radiani) dintre axa x pozitivă și vectorul (x, y) în planul XY. Argumentele trebuie să fie reale. Returnează o valoare între 0 și 2π. ceil(3.25) = 4 floor(3.25) = 3 mantissa (x) := x − floor(x) mantissa (3.45) = 0.45 Rotunjire tradițională: rotunjire (x) := if(mantissa (x)< 0.5, floor(x) , ceil(x)) roundoff (3.46) = 3 roundoff (3.56) = 4 Рис. 14. Создание функций округления На рис. 14 показано, как из этих функций могут быть сформированы функции округления. 4.4. Дискретное преобразование Фурье В Mathcad входят два типа функций для дискретного прямого и об- ратного преобразования Фурье: fft/ifft и cfft/icfft. Эти функции дискрет- ны: они берут в качестве аргументов и возвращают векторы и матрицы. Они не могут быть использованы с другими функциями. Используйте функции fft и ifft, если выполнены два следующих ус- ловия: аргументы вещественны, и вектор данных имеет 2m элементов. Первое условие необходимо, потому что функции fft/ifft используют тот факт, что для вещественных данных вторая половина преобразова- ния Фурье является комплексно сопряженной с первой. Mathcad отбра- сывает вторую половину вектора-результата. Это сохраняет время и память при вычислениях. Пара функций cfft/icfft не использует симметрию в преобразова- нии. По этой причине необходимо использовать их для комплексных данных. 41 Второе условие требуется, потому что пара функций fft/ifft исполь- зует высоко эффективный алгоритм быстрого преобразования Фурье. Для этого вектор аргумента, используемого с fft, должен иметь 2m эле- ментов. В функциях cfft/icfft использован алгоритм, который допускает в качестве аргументов как матрицы, так и векторы произвольного раз- мера. Когда эта пара функций используется с матрицей в качестве аргу- мента, вычисляется двумерное преобразование Фурье. Следует иметь в виду, что если для прямого преобразования исполь- зована функция fft, то для обратного преобразования необходимо ис- пользовать функцию ifft. Аналогично используются функции cfft/icfft. 4.5. Преобразование Фурье в вещественной области Для вещественных векторов с 2m элементами предпочтительно ис- пользовать функции fft/ifft. Функция fft(v) возвращает дискретное пре- образование Фурье, векторный аргумент которой можно интерпретиро- вать как результат измерений через равные промежутки времени некоторого сигнала. Вектор v должен содержать 2m элементов. Резуль- тат – комплекснозначный вектор размерности 1 + 2m–1. Если v имеет размерность, отличную от 2m, Mathcad выдает сообщение об ошибке "неверный размер вектора". Элементы вектора, возвращаемого fft, вычисляются по формуле n −1 ∑ vk e 2 πi (j n) k . 1 Cj = n k =0 В этой формуле n – число элементов в v, i – мнимая единица. Эле- менты в векторе, возвращенном функцией fft, соответствуют различ- ным частотам. Чтобы восстановить фактическую частоту, необходимо знать частоту измерения исходного сигнала. Если v есть n-мерный век- тор, переданный функции fft, и частота измерения исходного сигнала – fs, то частота, соответствующая Ck k fk = fs. n Обратите внимание, что это делает невозможным обнаружить часто- ты выше частоты измерения исходного сигнала. Это ограничение, нала- гаемое не Mathcad, а самой сутью проблемы. Чтобы правильно восста- новить сигнал по его преобразованию Фурье, необходимо произвести 42 i:= 0 .. 63 xi:= sin π⋅ + rnd (1) − 0.5 i Формирование сигнала: 10 Применяется комплексное преобразование Фурье: c:= fft(x) N:= last (c) N = 32 Обращение преобразования Фурье: z:= ifft(c) N2:= last (z) N2 = 63 j:= 0 .. N k:= 0 .. N2 Графическое представление сигнала zk = xj = 2 –0.499 –0.499 2.34·10 –3 2.34·10–3 0.673 0.673 xi 0 0.659 0.659 1.274 1.274 0.674 0.674 –2 0 20 40 60 80 1.162 1.162 i 0.613 0.613 Фурье-образ 0.179 0.179 4 –0.044 –0.044 0.489 0.489 –0.69 –0.69 cj 2 –1.079 –1.079 –0.777 –0.777 –0.849 –0.849 –1.334 –1.334 0 0 10 20 30 40 j Рис. 15. Быстрые пр6еобразования Фурье в Mathcad 43 измерения исходного сигнала с частотой, по крайней мере, вдвое боль- шей, чем ширина полосы частот. Подробное обсуждение этой пробле- мы содержится в специальных курсах. Функция ifft(v) возвращает обратное дискретное преобразование Фурье. Вектор v должен иметь 1 + 2m элементов, где m – целое. Резуль- тат есть вектор размерности 2m+1. Аргумент v – вектор, подобный созданному функцией fft. Чтобы вы- числить результат, Mathcad сначала создает новый вектор w, комплекс- но сопряженный v, и присоединяет его к вектору v. Затем Mathcad вы- числяет вектор d, элементы которого вычисляются по формуле n −1 ∑ wk e−2πi(j n)k . 1 dj = n k =0 Это та же самая формула, что и для fft, кроме знака минус в функции экспоненты. Функции fft и ifft – точные обращения. Для всех веще- ственных v справедливо ifft(fft(v)) = v. Пример использования прямого и обратного преобразований Фурье приведен на рис. 15. 4.6. Альтернативные формы преобразования Фурье Определения преобразования Фурье, рассмотренные выше, не явля- ются единственно возможными. Например, часто используются следу- ющие определения прямого и обратного преобразований Фурье: n n ∑ f (τ)e−2πi(ν n)τ ; f (τ) = ∑ F (ν) e () . 1 2 πi τ / n ν F (ν) = n τ=1 v =1 Эти определения реализованы во встроенных функциях FFT/IFFT и ICFFT. Они отличаются от быстрого преобразования Фурье следующим: вместо коэффициента 1 n перед обеими формулами стоит коэф- фициент 1/n и коэффициент 1 в обратном преобразовании; знак минус появляется в показателе экспоненты прямого преобразо- вания и исчезает в формуле обратного. 4.7. Кусочно-непрерывные функции Кусочно-непрерывные функции полезны для управления ветвлени- ями и остановками вычислительных процессов. Имеются пять функций 44 Использование условных операторов 2 x:= −2 , − 1.8 .. 2 f (x) := x − 1 g (x) := if(f (x) >0 , f (x) , 0) g(x) este egal cu f(x) când f(x) > 0, în caz contrar 0 4 4 2 f (x) g(x) 2 0 2 0 2 0 2 2 0 2 x x h (x) := if(x ≥ 1 , f (x) , − f (x)) altfel –f(x) 5 h(x) 0 Continuați calculele până când condiția este îndeplinită 5 2 0 2 2 întrebare − A< err x −2 N:= 100 i:= 0 .. N a:= 1000 quess 0:= 10 err:= 10 quess i + a quess i quess i+ 1:= until (quess i) − a − err , 2 2 N2:= last (quess) − 1 j:= 0 .. N2 j= quess j = (quess j)2 = 0 10 100 Число итераций N2 = 5 1 55 3.025·10 3 answer:= quess N2 2 36.591 1.339·10 3 3 31.96 1.021·10 3 answer = 31.623 4 31.625 1·10 3 5 31.623 1·10 3 Рис. 16. Условные выражения в Mathcad 45 Mathcad, относящихся к этому классу. Функция if полезна для выбора одного из двух значений, определяемого условием. Ступенчатая функ- ция Хевисайда Ф(х) и символ Кронекера δ(m, n) во многом аналогичны функции if. Функция until используется, чтобы управлять процессом итераций. Функция if(cond, tval, fval) возвращает значение tval, если cond отли- чен от 0 (истина) и возвращает fval, если cond равен 0 (ложь). Обычно в качестве аргумента cond выбирается булево выражение вида w = z, x >y, x< y, x ≥ y, x ≤ y, w ≠ z. Можно объединять булевы операторы, чтобы записать более сложные условия. Например, условие (x < 1) ⋅ (x >0) acționează ca „și” logic, care returnează 1 numai dacă x este între 0 și 1. În mod similar, expresia (x > 1) + (x< 0) действует подобно логическому "или", возвращающему 1, если x >1 sau x< 0, и 0, если x заключено между 0 и 1. Функция until (x, z) возвращает z, пока выражение x не становится отрицательным; должно содержать дискретный аргумент. Функция until позволяет останавливать вычисления для последовательных значений дискретного аргумента. Функция until полезна в итеративных процес- сах с определенным условием сходимости. На рис. 16 приведены примеры использования функций if и until. Функция Хевисайда эквивалентна следующей функции: Ф (x) := if (x < 0,0,1) Символ Кронекера δ(m, n) возвращает 1, если m = n; иначе 0. Оба аргумента должны быть целочисленными. Символ Кронекера эквива- лентен функции δ (m, n) := if (m = n,1,0) Ступенчатая функция Хевисайда может быть использована для со- здания импульса шириной w: pulse (x, w) := Ф (x) − Ф (x − w) Можно определить также две полезные функции lowpass и highpass. Они обе являются фильтрами – умножение на них какого-либо сигнала 46 вырезает из этого сигнала кусок вокруг точки x, имеющий ширину 2w. Разница состоит в том, что lowpass оставляет только вырезанный ку- сок, highpass – все, кроме вырезанного куска. lowpass (x, w) := pulse (x+w, 2 ⋅ w) highpass (x, w) := 1 − pulse (x+w, 2 ⋅ w) 4.8. Статистические функции Для вычисления статистических оценок случайных совокупностей чисел в Mathcad могут использоваться следующие функции: mean(A) – возвращает среднее значение элементов массива А раз- мерности m × n по формуле m −1 n −1 ∑ ∑ Aij ; 1 mean(A) = mn i =0 j =0 var(A) – возвращает дисперсию элементов массива А размерности m × n согласно формуле m −1 n −1 ∑ ∑ Aij − mean(A) 1 2 var(A) = ; mn i =0 j =0 stdev(A) - возвращает среднеквадратичное отклонение (квадратный корень из дисперсии) элементов m × n массива А stdev(A) = var(A). 4.9. Плотности распределения вероятности Эти функции показывают отношение вероятности того, что случай- ная величина попадает в малый диапазон значений с центром в задан- ной точке, к величине этого диапазона. В Mathcad имеются функции семнадцати плотностей вероятностей. Отметим только некоторые из них: dnorm(x, µ, σ) – возвращает плотность вероятности нормального рас- пределения 1 (x − µ) 2 dnorm(x, µ, σ) = exp − , 2πσ 2σ 2 47 в котором µ и σ есть среднее значение и среднеквадратичное отклоне- ние, σ >0; dunif(x, a, b) – calculează densitatea de probabilitate a unei distribuții uniforme 1, x ∈ , dunif(x, a, b) = b − a 0, x ∉ în care a și b sunt granițe interval de puncte, a< b. 4.10. Функции распределения Эти функции возвращают вероятность того, что случайная величи- на меньше или равна определенному значению. Функция распределе- ния вероятности – это функция плотности вероятности, проинтегриро- ванная от минус бесконечности до определенного значения. Приведем две из них: pnorm(x, µ, σ) – возвращает функцию нормального распределения со средним µ и среднеквадратическим отклонением σ (σ >0); punif(x, a, b) – returnează funcția de distribuție uniformă. a și b sunt valorile la limită ale intervalului (a< b). Mathcad имеет ряд функций для генерирования случайных чисел, имеющих разнообразные распределения вероятностей. Приведем две из них: rnorm(m, µ, σ) – возвращает вектор m случайных чисел, имеющих нормальное распределение (σ >0); runif(m, a, b) – returnează un vector de m numere aleatoare având o distribuție uniformă, în care a și b sunt punctele limită ale intervalului (a< b). Остальные встроенные статистические функции и их описания мож- но посмотреть, выбрав команду Функция из меню Вставка. 4.11. Интерполяция и функции предсказания Интерполяция заключается в использовании значений некоторой функции, заданных в ряде точек, чтобы предсказать значения между ними. В Mathcad можно или соединять точки данных прямыми линия- ми (линейная интерполяция) или соединять их отрезками кубического полинома (кубическая сплайн-интерполяция). 48 В отличие от функций регрессии, обсуждаемых в следующем разде- ле, функции интерполяции определяют кривую, точно проходящую че- рез заданные точки. Из-за этого результат очень чувствителен к ошиб- кам данных. Если данные зашумлены, следует рассмотреть возможность использования регрессии вместо интерполяции. Для линейной интерполяции используется функция linterp(vx, vy, x), которая по векторным данным vx и vy возвращает линейно интерполи- руемое значение y, соответствующее третьему аргументу x. Аргументы vx и vy должны быть векторами одинаковой длины. Вектор vx должен содержать вещественные значения, расположенные в порядке возраста- ния. Эта функция соединяет точки данных отрезками прямых, созда- вая, таким образом, ломаную линию. Интерполируемое значение для конкретного x есть ордината y соответствующей точки ломаной. Пример линейной интерполяции показан на рис. 17. Кубическая сплайн-интерполяция позволяет провести кривую через набор точек таким образом, что первые и вторые производные кривой непрерывны в каждой точке. Эта кривая образуется путем создания ряда кубических полиномов, проходящих через наборы из трех смежных то- чек. Кубические полиномы состыковываются друг с другом, чтобы об- разовать одну кривую. Чтобы провести кубический сплайн через набор точек: создайте векторы vx и vy, содержащие координаты x и y, через кото- рые нужно провести кубичный сплайн. Элементы vx должны быть рас- положены в порядке возрастания; вычислите вектор vs:=cspline(vx, vy). Вектор vs содержит вторые про- изводные интерполяционной кривой в рассматриваемых точках. Чтобы найти интерполируемое значение в произвольной точке, ска- жем х0, вычислите interp(vs, vx, vy, x0), где vs, vx и vy – векторы, опи- санные ранее. Обратите внимание, что можно сделать то же самое, вычисляя interp(cspline(vx, vy),vx,vy, x0). Пример использования кубической сплайн-интерполяции приведен на рис. 17 внизу. 49 Линейная интерполяция i:= 0 .. 5 VXi:=i VYi:=vd(1) VXi = VYi = –3 linterp(VX, VY, 1.5) = 0.389 0 1.268·10 1 0.193 linterp(VX, VY, 3.75) = 0.705 2 0.585 linterp(VX, VY, 4.1) = 0.758 3 0.35 4 0.823 x:= 0 , 0.1.. 5 5 0.174 1 linterp(VX , VY , x) 0.5 VYi 0 0 2 4 6 x , VX i Кубическая сплайн-интерполяция i:= 0 .. 5 VXi:= i VYi:= rnd (1) VS:= lspline (VX, VY) interp (VS, VX, VY, 1.5) = 0.188 interp (VS, VX, VY, 3.75) = 0.868 interp (VS, VX, VY, 4.1) = 0.989 1 VYi = 0.71 interp(VS , VX , VY , x) 0.304 0.5 VYi 0.091 0.147 0.989 0 0 2 4 6 0.119 x , VX i Рис. 17. Примеры интерполяции 50
