NUMERICKÁ ANALÝZA PROCESů

1. NUMERICKÁ ANALÝZA PROCESů NAP3 Modely „černých krabiček“ a experimentální identifikace charakteristik lineárních systémů: Přenosové charakteristiky systému. Volterrova rovnice a identifikace přenosu. Konvoluce a dekonvoluce. Fourierova analýza (spojitá a diskrétní Fourierova transformace). Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

2. x(t) y(t) x y t [s] Přenosové charakteristiky NAP3 Do kategorie modelů „černé krabičky“ patří i modely přenosových charakteristik systémů, které vyjadřují vztah mezi vstupy x(t) a výstupy y(t) systému, např. vztah mezi průběhem vlhkosti vzduchu (t) a vlhkostí sušeného materiálu X(t) vztah mezi časovým průběhem teploty potraviny a koncentrací mikroorganizmů vztah mezi časovými průběhy koncentrací na vstupu a výstupu průtočného reaktoru vztah mezi časovým průběhem zavírání ventilu a tlakem v potrubí (hydraulický ráz) vztah mezi průběhy teploty kapaliny a záznamem plášťového termočlánku

3. x(t) y(t) E(t) x y E t [s] Přenosové charakteristiky NAP3 U nelineárních systémů může i nepatrná změna vstupního signálu vyvolat nestability a tím se dostáváme do komplikovaně strukturované říše deterministického chaosu, podivných atraktorů (viz příští přednáška)…. Raději se jí zatím vyhneme: u lineárních systémů je to jednodušší, protože vztah mezi vstupní a výstupní funkcí lze vždy vyjádřit integrálnímvztahem, tzv. konvolucí slovy: odezva y(t) je konvolucí vzruchu x(t) a přenosové funkce E(t), symbolický zápis y=E*x proč to tak je viz následující folie…

4. x=(t) E t [s] Přenosové charakteristiky NAP3 Uvažujme speciální případ vstupní funkce x(t) ve tvaru nekonečně krátkého impulzu o jednotkové ploše. To je tzv. delta funkce (t) rovná nule pro t0 a nekonečnu pro t=0. Její integrál od minus do plus nekonečna je ale 1, takže konvoluční integrál je přímo roven funkci E(t): proto se funkce E(t) nazývá impulzní odezva E(t) (t) Diracovu funkci lze vyjádřit i jako limitu Gaussovy funkce E(t) Experimentálně se impulzní odezva stanoví měřením odezvy systému na kratičký impulz, například mžikový nástřik značkovací látky do vstupního proudu. Impulzní odezva E(t) je pak měřená koncentrace značkovače na výstupu.

5. x(t) y(t) E(t) x()d E(t-)x()d t [s]  Přenosové charakteristiky NAP3 Obecnou vstupní funkci lze nahradit sekvencí nekonečně mnoha nekonečně krátkých pulzů s plochou x()d. Každému takovému impulzu odpovídá impulzní odezva, která je ale posunutá o čas  a násobená plochou vstupního pulzu x()d (hovoříme o lineárních systémech, kde výstup je přímo úměrný vstupu). Součtem těchto elementárních odezev je pak konvoluční integrál E * x to plyne z toho, že impulzní odezva i x(t) jsou pro záporné časy nulové

6. Řešení Volterrovy rovnice NAP3 Jsou dva základní typy problémů Konvoluce, když známe charakteristiku systému E(t) a časový průběh vstupu x(t). Odezvu y(t) stanovíme integrací, nepříjemné je ale to, že pro každý čas t je třeba počítat jiný integrál. Dekonvoluce, když máme například změřené vstupy i výstupy x(t) a y(t) a chceme z nich stanovit impulzní odezvu E(t) (identifikace systému). Pak je třeba řešit integrální rovnici. Nástrojem pro řešení obou těchto problémů je Fourierova transformace Leger

7. Fourierovatransformace NAP3 Funkce x(t),y(t),E(t) by bylo možné aproximovat běžnými polynomy (viz regresní modely) a konvoluci vyjádřit vztahy mezi koeficienty těchto polynomů. Jenže: každý polynom „uteče“ k nekonečnu pro t a regresní polynomy vyššího stupně než 7 ani nemají smysl (jejich koeficienty nejsou optimalizací určeny přesně, protože funkce 1,t,t2,t3,…,t7,t8 jsou si dost podobné a nejsou dostatečně lineárně nezávislé. Nakonec vždy převáží vliv zaokrouhlovacích chyb.). Místo těchto funkcí je proto užitečnější použít ortogonální funkce, např. ortogonální polynomy nebo goniometrické funkce a aproximovat funkce x,y,E Fourierovými rozvoji Ale co to vůbec znamená, když řekneme, že funkce Pm(t) a Pn(t) jsou ortogonální? Duchamp Vlastně byste to měli vědět, viz předchozí přednáška☺.

8. Ortogonalita NAP3 Funkce Pm(t) a Pn(t) jsou ortogonální, když integrál jejich součinu (přes určitý interval) je roven nule. integrálu součinu se říká skalární součin dvou funkcí Například goniometrické funkce (sin, cos) jsou ortogonální v intervalu (-1,1). Ortogonalita Pm(t)=cos 2mtplyne z funkce cos jsou dokonce ortonormální, protože integrál jejich druhé mocniny je roven jedné. pro m=n, jinak=0 Dokaž!!! Zcela stejně lze dokázat ortogonalitu sin 2mt. Z Eulerovy formule pak plyne i ortogonalita komplexní funkce i-imaginárníjednotka

9. Fourierův rozvoj NAP3 Lineárníkombinaceortogonálních funkcí Pj(t) je Fourierovým rozvojem libovolné funkce c(t) přičemž koeficienty rozvoje jsou funkcí c(t) určeny jednoznačně (a přiřazení se nazývá diskrétní Fourierova transformace) Dokaž!!! Koeficienty rozvoje lze jak vidno stanovit bez řešení soustavy lineárních algebraických rovnic (porovnej s příkladem regresní analýzy a obyčejných polynomů). Použití ortogonálních funkcí místo P0(t)=1, P1(t)=t, P2(t)=t2,… je často výhodné i v regresní analýze, vezme se prostě jen několik prvních členů Fourierovy řady a jako ortogonální funkce např. Legendreovy, Čebyševovy, Hermiteovy nebo Laguerrovy zobecněné polynomy (odlišnost vyplývá z různé definice skalárního součinu, třeba u Laguerrových polynomů je integrační interval <0,>, zatímco u Legendreových <-1,1>. Volba vhodného typu ortogonální funkce je dána požadovaným definičním intervalem regresní funkce: když bude<0,>, použijeme Laguerrovy, pro interval<-,> Hermiteovy, pro konečný interval Legendreovy nebo Čebyševovy polynomy). V dalším textu se však budeme zabývat jen Fourierovými rozvoji na bázi goniometrických funkcí (za okamžik dokážeme jejich ortogonalitu i na jiném intervalu než <-1,1>)

10. Fourierův rozvoj NAP3 Sumu nekonečného počtu členů Fourierovy řady lze nahradit integrálem. Místo přirozených čísel (indexů j) se použije spojitá veličina, frekvence f. Například Fourierova řada Fourierův integrál čímž se dostáváme k pojmu spojité integrálníFourierovy transformace. FT je komplexní funkce i když je c(t) funkce reálná Spojitá dopředná Fourierova transformace funkce času na funkci frekvence Úžasnou výhodou FT je praktická shodnost dopředné a zpětné transformace (použije se stejný podprogram)! Zpětná Fourierova transformace z frekvenčního do časového prostoru

11. Fourierova transformace NAP3 Důkaz (nepřesný, spíše ilustrativní) dále dokážeme, že funkce e2if jsou ortonormální a jejich skalární součin je delta funkce (-t) Substituce 2a(-t)=x dx=2 ad

12. Fourierova transformace NAP3 Výkonová spektrální hustota (PSD Power Spectral Density) Energie vysokých frekvencí (šum) Energie nízkých frekvencí Graf PSD ukáže dominantní frekvence v signálu a umožní nastavení frekvenčních filtrů pro odstranění šumu

13. Fourierova transformace NAP3 • Půvab Fourierovy transformace tkví v tom, že • integrál součinu (typu konvoluce nebo korelace) převádí na prostý součin transformovaných funkcí • derivaci funkce převádí na násobení transformované funkce frekvencí (viz další přednáška) • dá se snadno provést zpětná transformace (snadný přechod z časového do frekvenčního prostoru a vice versa)

14. Fourierova transformace NAP3 Konvoluce dvou funkcí (určuje jak se změní funkce x(t), když projde filtrem y(t)) Korelace dvou funkcí (vyjadřuje míru shody funkcí x,y při vzájemném posunutí o čas t) Střední doba vstupu x(t)=0.5 Střední doba y(t)=2.1 Střední doba kros-korelace Rxy(t)=1.6 reprezentuje časový posun mezi x(t) a y(t)

15. FT aplikace NAP3 • Výpočet odezvy systému na daný vstupní signál x(t) • Stanoví se FT vstupu a impulzní odezvy (dopředná FT) • FT výstupu je součin FT (součin komplexních čísel) • Výstup se stanoví inverzní FT transformací • Identifikace systému ze změřeného vstupu a výstupu x(t),y(t) • Stanoví se FT vstupu a výstupu (dopředná FT) • FT impulzní odezvy je podíl FT (podíl komplexních čísel) • Impulzní odezva se stanoví inverzní FT transformací

16. 2 Vzorkování signálu NAP3 V reálu (i v numerice) pracujeme jen s bodovou reprezentací funkcí vzorkovaných obvykle s konstantním časovým krokem. Funkci reprezentuje N-bodů data N-frekvencí (sudé číslo) Nyquistova frekvence Nyquistova frequenceje maximální frekvence rozlišitelná v datech vzorkovaných s krokem

17. D D D D t: t =0,t = ,t =2 ,3 ,......,t = (N - 1) , 0 2 N - 1 1 1 1 1 1 1 - - + - f = ,.. f = , f = f: f = ,f = ( ) ,.. f = , f =0, 1 N/2 - 1 N/2 - N/2 - N/2+1 - 1 0 D D D 2 2 N N 1 1 1 1 1 - ( ) D D D N 2 N 2 Diskrétní FT(Fourierova Transformace) NAP3 Diskrétní DFT kde - N 1 å ~ p = D 2 ikn / N c ( f ) c ( t ) e n k = k 0 n=-N/2,-N/2+1,…,0,1,…,N/2 DFTmá stejné vlastnosti (konvoluce, korelace) jako spojitá Fourierova Transformace.

18. Diskrétní FT(Fourierova Transformace) NAP3 N-Fourierových koeficientů reprezentuje funkci c(t) definovanou od - do +, která prochází zadanými body c0, c1,…cN-1 a je periodická s periodou N. Fourierovy koeficienty kladných a záporných Nyquistových frekvencí jsou totožné: ale sin(k)=0 pro libovolné k

19. Záporné frekvence? NAP3 Při práci s DFT působí někdy problém zjistit, jaká frekvence odpovídá indexu vypočteného Fourierova koeficientu. Příklad: Bylo navzorkováno 8 dat s frekvencí 1 Hz, tj. =1s. Nyquistova (maximální) frekvence je 0.5 Hz. Těmto osmi číslům (c0,…,c7) odpovídá vektor devíti koeficientů Fourierovy transformace (c-4,c-3,…c0,…c3,c4) pro frekvence (-4/8,-3/8,-2/8,-1/8,0,1/8,2/8,3/8,4/8). Jenomže Fourierovy koeficienty kladných a záporných Nyquistových frekvencí jsou totožné (c-4=c4) takže těch koeficientů je vlastně jen osm a obvykle se uspořádají do vektoru, který začíná nulovou frekvencí, tj. (c0, c-1, c-2,c-3,c4=c-4, c3, c2 ,,c1), viz MATLAB. Cosinové složky kladných a záporných frekvencí jsou totožné Sinusové složky kladných a záporných frekvencí mají jen obrácené znaménko Když se jedná o FT reálné funkce úplně stačí uvažovat jen kladné (nebo záporné) frekvence, protože Fourierovy koeficienty kladných a záporných frekvencí nejsou nezávislé (jsou komplexně sdružené). FT je transformace 1:1, tzn. 8 reálných čísel se transformuje na jiných 8 reálných čísel, což jsou reálné a imaginární části pouze čtyř Fourierových koeficientů.

20. tk fn všimněte si, že FFT je transformace N-čísel ck na N-koeficientů, která odpovídá vzorkovacímu intervalu =1s. Na to je třeba si dát pozor při výpočtu konvoluce a korelace. Diskrétní FFT (Fast FT) NAP3 Na této foliii se pokouším vysvětlit souvislosti mezi spojitou a diskrétní FT (když si pamatujete spojitou snadno odvodíte diskrétní). Jen si dejte pozor na roli časového kroku . Dopředná a zpětná FFT mezi DFT a FFT není žádný rozdíl-Fast Fourier Transform označuje jen velmi rychlý algoritmus vyčíslení sumy, který se dá použít tehdy, když N je mocnina dvou, tedy např. 512,1024,2048,…

21. Diskrétní FFT (Fast FT) NAP3 Parsevalův teorém Slovy: výkon signálu počítaný z časového průběhu je stejný jako výkon, počítaný z frekvenčního průběhu (proč výkon? u elektrického signálu je výkon na konstantním odporu úměrný kvadrátu napětí).

22. MATLAB příklady aplikací FFT NAP3 Braque

23. FFT funkce v MATLABu NAP3 Dopředná transformace vektoru c cf=fft(c,N) Vektor vypočtených Fourierových koeficientů cf(1),…,cf(N) Vektor vzorků dat c(1),c(2),….,c(N) Inverzní transformace c=ifft(cf,N)

24. FFT funkce v MATLABu NAP3 Na následující folii porovnáme definici Fourierovy transformace uvedenou dříve (s kladnými a zápornými frekvencemi) s definicí, kterou používá MATLAB.

25. MATLAB frekvence= f Nyquistova frekvence index koeficientu v MATLABu FFT funkce v MATLABu NAP3 FFT (n=-N/2 až N/2 index vyjadřuje přímo frekvenci) MATLAB (n=1 až N) takto jsme si definovali FT Tak co je správně? a tak je to v MATLABu FFT Vztahy pro výpočet Fourierových koeficientů dle původní definice FFT a MATLABu nejsou úplně stejné a mohlo by se zdát, že MATLAB operuje s frekvencemi, které jsou i vyšší než Nyquistova frekvence. Oba vztahy jsou ale správně a vypočtené koeficienty jsou totožné (jen jinak uspořádané). Je to důsledek periodicity goniometrických funkcí (exp(2i(k-1))=1), takže např.

26. FFT MATLABpříklad PSD NAP3 Běžnou aplikací FFT je nalezení dominantních frekvencí signálu, utopených v náhodném šumu. Uvažujme data vzorkovaná s frekvencí 1000 Hz. Vytvořte signál obsahující složky s frekvencí 50 Hz a 120 Hz na které je superponován náhodný šum s nulovou střední hodnotou: t = 0:0.001:0.6;x = sin(2*pi*50*t)+sin(2*pi*120*t); y = x + 2*randn(size(t)); Je asi těžké rozeznat základní frekvence v zaznamenaném signálu. Převedeme ho tedy do frekvenční oblasti použitím N=512ti bodů a FFT: Y = fft(y,512); Výkonové spektrum Pyy = Y.* conj(Y) / 512; Smysl má jen první polovina bodů (symetrie je daná tím, že y(t) byla reálná funkce a pro reálnou funkci jsou koeficienty záporných frekvencí komplexně sdružené s koeficienty kladných frekvencí) Normální distribucenáhodných čísel (střední hodnota=0, variance=1) DFT z 512tivzorků (mocnina dvou) Dominantní frekvence PSD

27. Dominantní frekvence FFT MATLABpříklad PSD NAP3 V grafu PSD jsou na horizontální ose jen indexy Fourierových koeficientů, a my vidíme, že ty dominantní frekvence odpovídají Fourierovým koeficientům zhruba v rozsahu indexů 25 až 30 (první vrchol) a 60 až 65 (druhý vrchol). Frekvence odpovídající indexům plynou z definice inverzní transformace f = 1000*(0:256)/512; tady je kladná a současně záporná Nyquistova frekvence případ, kdy pouze k-tý Fourierův koeficient je nenulový k je index vektoru vypočítaných Fourierových koeficientů, začínající od 1 Pro N=512 a =0.001s je tudíž f25=47 Hz, f30=57 Hz, f60=115 Hz, f65=125 Hz

28. PSD filtr šumu NAP3 Nejjednodušší způsob filtrace šumu je potlačení Fourierových koeficientů vyšších frekvencí, např. komponent Y(65),Y(66),…Y(449) (přiřaď těmto prvkům nuly). Výsledná PSD je n 514-n prosím, vezměte v úvahu, že poslední koeficient Y(512) neodpovídá nulové, ale nejmenší frekvenci. Vynulované fourierovykoeficienty Originální signál pro srovnání Filtrovaný signál je rekonstruován inverzní FT y=ifft(Y,512)

29. FFT konvoluce a korelace NAP3 for i=1:512 c1(i)=t(i)*exp(-2*t(i)); c2(i)=t(i)^4*exp(-t(i)*3); end f1=fft(c1,512); f2=fft(c2,512); Fourierovy koef. pro konvoluci a korelaci for i=1:512 c12(i)=f1(i)*f2(i); r12(i)=f1(i)*conj(f2(i)); end cc=ifft(c12,512)*0.001; rr=ifft(r12,512)*0.001; efekt zrcadlení (periodicity FT) při konvoluci a korelaci se počítá součin FT a výsledek je třeba násobit , při dekonvoluci se počítá podíl FT a výsledek se musí dělit  Inverzní Fourierova transformace Efekt zrcadlení je projev periodicity diskrétní Fourierovy transformace. Dá se potlačit jen volbou hodně dlouhého intervalu vzorkování.

30. FFT identifikace E(t) NAP3 Následující folie ilustrují základní techniku zpracování experimentů s radioizotopovými indikátory, používanou pro identifikaci impulsní odezvy průtočných systémů. Desítky konkrétních aplikací v průmyslu (fluidní spalování, pece, reaktory, absorbéry, kolony, filtry, aerační nádrže, mlýny, zásobníky) jsou uvedeny v monografii Thýn J.,Žitný R.:Analysis and diagnostics of industrial processes by radiotracers. ČVUT Praha 2000

31. x(t) y(t) (t) x(t) y(t) FFT identifikace E(t) NAP3 Uvažujme modelový systém, tvořený kaskádou 4 ideálně promíchávaných nádob. Impulzní odezva E(t) takové kaskády se dá vyjádřit v analytickém tvaru (třeba Laplaceovou transformací, ale to teď není důležité) N-počet míchaných nádob, tm je střední doba prodlení (celkový objem/průtok) Na vstupu budeme měřit signál x(t) (třeba koncentraci značkovací látky) a na výstupu signál y(t). Tyto signály vytvoříme uměle: x(t) jako odezvu na mžikový nástřik do systému tvořenému dvěma mísiči (N=2). Výstup y(t) pak musí být impulzní odezva kaskády šesti identických misičů (N=2+4) se střední dobou odpovídající součtu středních dob tm vstupního a identifikovaného systému. Takto vygenerované signály x(t) a y(t) pak zatížíme náhodným šumem anonymní funkce @ se třemi parametry. Linspace je standardní funkce. dt=0.01; t=0:dt:10.23; n=2; tm=1; x=n^n.*t.^(n-1)/(factorial(n-1)*tm^n).*exp(-n.*t/tm); n=6; tm=3; y=n^n.*t.^(n-1)/(factorial(n-1)*tm^n).*exp(-n.*t/tm); x=x+0.01*randn(size(t)); y=y+0.005*randn(size(t)); úplně totéž se dá napsat v MATLABu na 4 řádky e=@(t,tm,n) n^n.*t.^(n-1)./(factorial(n-1)*tm^n).*exp(-n.*t/tm) t=linspace(0,10.23,1024); x=e(t,1,2)+0.01*randn(size(t)); y=e(t,3,6)+0.005*randn(size(t));

32. FFT identifikace E(t) NAP3 Impulzní odezvu identifikovaného systému získáme dekonvolucí, tj. aplikací FFT na signály vstupu i výstupu (vektory Fourierových koeficientů xf(1:1024), yf(1:1024)), Fourierovy koeficienty funkce E(t) jsou prostě jen podíl koeficientů yf(i)/xf(i) a impulzní odezva se získá zpětnou FFT xf=fft(x); yf=fft(y); ef=yf./xf; e=ifft(ef)/dt; plot(t,e) Výsledek je divoce rozkmitaná funkce mající jen pramálo společného s teoretickou impulzní odezvou (modrá křivka n=4;tm=2;et=n^n.*t.^(n-1)/(factorial(n-1)*tm^n).*exp(-n.*t/tm); ). Dekonvoluce je totiž tzv. špatně podmíněný problém, kdy i malá odchylka výstupního signály y(t) zapříčiní velkou odchylku řešení integrální rovnice. Šum výsledku je možné potlačit vynulováním vyšších frekvencí, např. teoretická impulzní odezva xf=fft(x); yf=fft(y); ef=yf./xf; ef(17:1024-16)=0; e=ifft(ef)/dt;

33. FFT identifikace E(t) Wiener NAP3 V předchozím případu jsme filtrovali šum z výsledné impulzní odezvy. Možná trochu správnější je filtrovat šum jen z výstupní funkce y(t) a to tak, že potlačíme vyšší frekvence v koeficientech Fourierovy transformace yf tzv.Wienerovou filtrací (její teoretický popis najdete např. v Teukolski et al.: Numerical recipes). Opět se vychází z grafu PSD zašuměné odezvy, z něhož se odhadnou frekvence „užitečného“ signálu a frekvence, které je třeba potlačit. odhadnutá PSD „čistého“ signálu log PSD to je jen ukázka prvních 100 frekvencí PSD. Červená čára je extrapolací šumu a umožňuje odseparovat šum a užitečný signál. PSD „zašuměného“ signálu touto korekční funkcí se pak násobí Fourierova transformace změřeného výstupu

34. x(t) xp(t) y(t) (t) FFT posunutí signálu NAP3 Transformace „posunutého“ signálu v čase se realizuje jen násobením FT exponenciálou. To usnadňuje modelování systémů v nichž se například objevuje dopravní zpoždění (a které působí problémy při numerickém řešení soustavy diferenciálních rovnic). FT funkce posunuté o  Odezva serie dvou mísičů Zona pístového toku (zdržení ) x(t) xp(t) dt=0.01; t=0:dt:10.23; for i=1:1024 f(i)=(i-1)/(1024*dt); end % vstup tm=1, N=2 n=2; tm=1; x=n^n.*t.^(n-1)/(factorial(n-1)*tm^n).*exp(-n.*t/tm); xf=fft(x); tau=-2; xfp=xf.*exp(2*pi*tau*complex(0,1).*f); yf=xfp.*xf; xp=ifft(xfp); y=ifft(yf)*dt; Vektor frekvencí y(t) Posunutí o tau

35. Ohřívač (zdroj náhodných pulzů) T1 T2 FFT vzájemná korelace NAP3 Korelace stimulovaného nebo náhodného signálu ze dvou míst(technicky třeba korelace signálu dvou termočlánků) Z korelační funkce lze vyhodnotit časový posun signálů a z něho pak třeba rychlost, průtok,… Ale třeba též deformace povrchu zatížené součásti snímané digitálními kamerami (image processing)

36. Ohřívač FFT vzájemná korelace NAP3 Příklad simulovaný v MATLABu Náhodný signál posunutý o 100 časových kroků z polohy tohoto vrcholku lze usuzovat na dobu průchodu a průtok

37. Co je třeba si pamatovat NAP3 Přednáška byla věnována Fourierově transformaci a jejím aplikacím. Zapamatujte si alespoň to, co je na následující folii

38. Co je třeba si pamatovat NAP3 Jak se změní signál x(t) po průchodu filtrem E(t) Korelace dvou signálů posunutých o čas t Co je konvoluce, dekonvoluce, přenos a korelace Fourierova transformace Fourierova transformace konvoluce a korelace Diskrétní Fourierova transformace Nyquistova frekvence 1/2