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

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 reaktoruvztah 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

Přenosové charakteristiky

x(t) y(t)

x

y

t [s]

NAP3 Přenosové charakteristiky

E(t)x(t) y(t)

x

y

t [s]

E

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ím vztahem, tzv. konvolucí

dxtEty )()()(

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…

Diracovu funkci lze vyjádřit i jako limitu Gaussovy funkce

NAP3 Přenosové charakteristiky

E(t)(t)

E(t)

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 t0 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):

( ) ( ) ( ) ( )y t E t d E t

proto se funkce E(t)

nazývá impulzní odezvax=(t)

t [s]

EExperimentá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.

2

lim)( nt

nent

NAP3 Přenosové charakteristiky

E(t)x(t) y(t)

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

0

( ) ( ) ( ) ( ) ( )t

y t E t x d E t x d

x()d

t [s]

E(t-)x()d

to plyne z toho, že impulzní odezva i x(t) jsou pro záporné časy nulové

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

Řešení Volterrovy rovnice

dxtEty )()()(

Leger

Fourierova transformace

Duchamp

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

0

( ) ( )j jj

E t E P t

Ale co to vůbec znamená, když řekneme, že funkce Pm(t) a Pn(t) jsou ortogonální?

Vlastně byste to měli vědět, viz předchozí přednáška ☺.

Například goniometrické funkce (sin, cos) jsou ortogonální v intervalu (-1,1).

Ortogonalita Pm(t)=cos 2mt plyne z1 1

1 1

11

1cos 2 cos 2 (cos 2 ( ) cos 2 ( ) )2

1 sin 2 ( ) sin 2 ( )[ ] 12 2 ( ) 2 ( )

mt ntdt m n t m n t dt

m n t m n tm n m n

pro m=n, jinak =0

Zcela stejně lze dokázat ortogonalitu sin 2mt. Z Eulerovy formule pak plyne i ortogonalita komplexní funkce 2( ) cos 2 sin 2i mt

mP t e mt i mt

Dokaž!!!

i-imaginární jednotka

NAP3 OrtogonalitaFunkce Pm(t) a Pn(t) jsou ortogonální, když integrál jejich součinu (přes určitý interval) je roven nule.

funkce cos jsou dokonce ortonormální, protože integrál jejich

druhé mocniny je roven jedné.

integrálu součinu se říká skalární součin dvou funkcí

Lineární kombinace ortogonálních funkcí Pj(t) je Fourierovým rozvojem libovolné funkce c(t)

0

( ) ( )j jj

c t c P t

NAP3 Fourierův rozvoj

přičemž koeficienty rozvoje jsou funkcí c(t) určeny jednoznačně (a přiřazení se nazývá diskrétní Fourierova transformace)

2

( ) ( )

( )j

jj

c t P t dtc

P t dt

jc~

Dokaž!!!0 1( ) , ,..., nc t c c c

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í 2( ) imt

mP t e (za okamžik dokážeme jejich ortogonalitu i na jiném intervalu než <-1,1>)

Spojitá dopředná Fourierova transformace funkce času na funkci frekvence

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

2 2

0

( ) ( ) ( )ijt iftj

j

c t c e c t c f e df

NAP3 Fourierův rozvoj

čímž se dostáváme k pojmu spojité integrální Fourierovy transformace.

dtetcfc ift2)()(~

2( ) ( ) iftc t c f e df

Zpětná Fourierova transformace z frekvenčního do časového prostoru

Úžasnou výhodou FT je praktická shodnost dopředné a zpětné transformace (použije

se stejný podprogram)!

FT je komplexní funkce i

když je c(t) funkce reálná

Fourierova řada Fourierův integrál

Důkaz (nepřesný, spíše ilustrativní)

NAP3 Fourierova transformace

2 2 2 2

2 ( ) 2 ( )2 ( )

( ) ( ( ) ) ( )( )

( ) ( ) ( ) ( )lim lim2 ( )

sin 2 ( )( ( ) )lim( )

if ift if ift

a ia t ia tif t

a aa

a a

c t c e d e df c e e df d

e ec e df d c di t

a tc dt

sin sin( ( )) ( ) ( )lim

2x x xc t dx c t dx c ta x x

dále dokážeme, že funkce e2if jsou ortonormální a jejich

skalární součin je delta funkce (-t)

Substituce 2a(-t)=x

dx=2 ad

Výkonová spektrální hustota (PSD Power Spectral Density)

22 )(~)(~)( fcfcfPc

Energie nízkých frekvencí

Energie vysokých frekvencí (šum)

NAP3 Fourierova transformace

Graf PSD ukáže dominantní frekvence v signálu a umožní nastavení frekvenčních filtrů pro odstranění šumu

Půvab Fourierovy transformace tkví v tom, že

1. integrál součinu (typu konvoluce nebo korelace) převádí na prostý součin transformovaných funkcí

2. derivaci funkce převádí na násobení transformované funkce frekvencí (viz další přednáška)

3. dá se snadno provést zpětná transformace (snadný přechod z časového do frekvenčního prostoru a vice versa)

NAP3 Fourierova transformace

dytxtCxy )()()( ( ) ( ) ( )xyC f x f y f

dytxtRxy )()()( )(~)(~)(~ * fyfxfRxy

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)

NAP3 Fourierova transformace

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í

NAP3 FT aplikace

)(~)( ),(~)( fEtEfxtx )(~)(~)(~ fEfxfy

( ) ( )y f y t

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í

( ) ( ), y( ) ( )x t x f t y f

( ) ( ) / ( )E f y f x f

( ) ( )E f E t

Vzorkování signáluV reálu (i v numerice) pracujeme jen s bodovou reprezentací funkcí vzorkovaných obvykle s konstantním časovým krokem .Funkci reprezentuje N-bodů dat a N-frekvencí (sudé číslo)

t: t0=0, t1=, t2=2, 3,......, tN-1= (N-1),

f: f-N/2=

21 ,f-N/2+1=

1)121(

N,.. f-1=

N1 , f0=0,

f1=N

1 ,.. fN/2-1=

1)1

21(

N, fN/2=

21

Nyquistova frekvence

Nyquistova frequence je maximální frekvence rozlišitelná v datech vzorkovaných s krokem

2

NAP3

Diskrétní FT (Fourierova Transformace)

DFT má stejné vlastnosti (konvoluce, korelace) jako spojitá Fourierova Transformace.

NAP3

2( ) ( ) iftc f c t e dt

Diskrétní DFT

1

0

/2)()(~N

k

Niknkn etcfc

kde

n=-N/2,-N/2+1,…,0,1,…,N/2

t: t0=0,t1 =,t 2=2 ,3,......,tN-1= (N-1),

f: f-N/2=

21 ,f-N/2+1=

1)121(

N,.. f-1=

N1 , f0=0, f1= N

1 ,.. fN/2-1= 1)1

21(

N, fN/2= 2

1

Nnfn ktk

Diskrétní FT (Fourierova Transformace)

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.

NAP3

1 1 1 12 ( /2)/

/20 0 0 0

( ) (cos( ) sin ) cosN N N N

ik N N ikN k k k k

k k k k

c f c e c e c k i k c k

ale sin(k)=0 pro libovolné k

Fourierovy koeficienty kladných a záporných Nyquistových frekvencí jsou totožné:

1 1 1 12 ( /2)/

/20 0 0 0

( ) (cos sin ) cosN N N N

ik N N ikN k k k k

k k k k

c f c e c e c k i k c k

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.

0 1 2 3 4 5 6 7-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

4sin(2 ( ) )8t

0 1 2 3 4 5 6 7-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

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ů.

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

4cos(2 )8t

3cos(2 )8t )

842sin( t

3sin(2 )8t

Dopředná a zpětná FFT

1 12 2

0 0

/2 /22 2

/2 /2

1 1

kn k nN Ni iN N

n k kk k

kn k nN Ni iN N

k n nn N n N

c c e c e

c c e c eN N

NAP3 Diskrétní FFT (Fast FT)

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,…

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.

( ) ( ) /n n k kc f c c t c

fn

tk

dt

dtetcfc ift2)()(~

2( ) ( ) iftc t c f e df

N

df 1

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 .

Parsevalův teorém

1

0

21

0

2 |~|1||N

kk

N

kk c

Nc

NAP3 Diskrétní FFT (Fast FT)

dffcdttc 22 )(~)(

N

df 1dt kk cfc ~)(~

( )k kc t c

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í).

NAP3 MATLAB příklady aplikací FFT

Braque

FFT funkce v MATLABu

12 ( 1)

1 1

12 ( 1)

1 1

1 1(cos(2 ( 1) ) sin(2 ( 1) ))

1 1 1 1(cos(2 ( 1) ) sin(2 ( 1) ))

nN Ni kN

n k kk k

nN Ni kN

k n nn n

n nc c e c k i kN N

n nc c e c k i kN N N N

Dopředná transformace vektoru c

cf=fft(c,N)

Inverzní transformace

c=ifft(cf,N)

Vektor vzorků dat c(1),c(2),….,c(N)

Vektor vypočtených Fourierových koeficientů cf(1),…,cf(N)

NAP3

FFT funkce v MATLABuNAP3

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.

FFT funkce v MATLABuNAP3

( 1)2

1

k nN iN

n kk

c c e

( 1)( 1)2

1

k nN iN

n kk

c c e

MATLAB (n=1 až N)FFT (n=-N/2 až N/2 index vyjadřuje přímo frekvenci)

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(2i(k-1))=1), takže např.

1 2 3 /2 /2 1 /2 2

0 1 2 /2 1 /2 /2 /2 1

.... .... .... = ....

N N N N

N N N N

c c c c c c cc c c c c c c

1

1 11 2 1 12 20 .... .... 2

cN N

N N N N N

MATLAB

frekvence=

f

index koeficientu v MATLABu

Nyquistova frekvence

)12

(12)1(2)12

(12)12

(12

22/

~N

Nki

k

kiNNki

k

NNki

kMATLAB ecececcN

FFT

takto jsme si definovali FT

a tak je to v MATLABu

FFT MATLAB příklad PSDBěž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ýchfrekvencí komplexně sdružené s koeficienty kladných frekvencí)

0 100 200 300 400 500 6000

10

20

30

40

50

60

70

80

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7-6

-4

-2

0

2

4

6

Normální distribuce náhodných čísel (střední hodnota=0,

variance=1)

DFT z 512ti vzorků (mocnina dvou)

PSD Dominantní frekvence

NAP3

FFT MATLAB příklad PSD

0 100 200 300 400 500 6000

10

20

30

40

50

60

70

80

Dominantní frekvence

NAP3

12 ( 1)

1

12

1

1( ( 1) ) n

nN i kN

n kk

ki tN

n k

c c eN

c t n c eN

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

Nkfk

1

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

tady je kladná a současně záporná Nyquistova

frekvence případ, kdy pouze k-tý Fourierův

koeficient je nenulový

f = 1000*(0:256)/512;

PSD filtr šumuNejjednodušší 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

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7-6

-4

-2

0

2

4

6

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5-4

-3

-2

-1

0

1

2

3

4

0 100 200 300 400 500 6000

10

20

30

40

50

60

70

80

Filtrovaný signál je rekonstruován inverzní FT

y=ifft(Y,512)

Originální signál pro srovnání

Vynulované fourierovy koeficienty

NAP3

n 514-nprosím, vezměte v úvahu, že poslední koeficient Y(512) neodpovídá nulové, ale nejmenší frekvenci.

FFT konvoluce a korelace

0 1 2 3 4 5 60

1

2

3

4

5

6

7

8x 10

-3

0 1 2 3 4 5 60

1

2

3

4

5

6

7

8x 10

-3

for i=1:512c1(i)=t(i)*exp(-2*t(i));c2(i)=t(i)^4*exp(-t(i)*3);endf1=fft(c1,512);f2=fft(c2,512);

for i=1:512c12(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;

Fourierovy koef. pro konvoluci a korelaci

efekt zrcadlení (periodicity FT)

NAP3

Efekt zrcadlení je projev periodicity diskrétní Fourierovy transformace. Dá se potlačit jen volbou hodně dlouhého intervalu vzorkování.

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

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

FFT identifikace E(t)NAP3

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));

(t)y(t)

x(t)

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é)

mtNtNm

NN

etN

tNtE /1

)!1()(

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

0 2 4 6 8 10 12-0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

x(t)

y(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));

anonymní funkce @ se třemi parametry. Linspace je

standardní funkce.

FFT identifikace E(t)NAP3

xf=fft(x);yf=fft(y);ef=yf./xf;e=ifft(ef)/dt;plot(t,e)

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

0 2 4 6 8 10 12-0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

xf=fft(x);yf=fft(y);ef=yf./xf;ef(17:1024-16)=0;e=ifft(ef)/dt;

0 2 4 6 8 10 12-15

-10

-5

0

5

10

15

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

( ) ( ) / ( )E f y f x f

FFT identifikace E(t) WienerNAP3

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.

0 10 20 30 40 50 60 70 80 90 100-15

-10

-5

0

5

10

2

2

|)(~||)(~|)(~

fyfyf c

PSD „zašuměného“

signálu

odhadnutá PSD „čistého“ signálu

touto korekční funkcí se pak násobí Fourierova

transformace změřeného výstupu

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.

FFT posunutí signáluNAP3

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).

0 2 4 6 8 10 12-0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

dt=0.01;t=0:dt:10.23;for i=1:1024 f(i)=(i-1)/(1024*dt);end% vstup tm=1, N=2n=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;

(t)y(t)x(t) xp(t)

)(~)()()(~ 2222 fcedtetcedtetcfc ififtififtp

FT funkce posunuté o

Odezva serie dvou mísičů

Zona pístového toku (zdržení )

x(t)xp(t)

y(t)

Vektor frekvencí

Posunutí o tau

Korelace stimulovaného nebo náhodného signálu ze dvou míst (technicky třeba korelace signálu dvou termočlánků)

dttTtTR )()()( 2112

T1 T2

Ohřívač (zdroj náhodných pulzů)

NAP3 FFT vzájemná korelace

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)

dttTtTR )()()( 2112

Příklad simulovaný v MATLABu

Ohřívač

Náhodný signál posunutý o 100 časových kroků

NAP3 FFT vzájemná korelace

z polohy tohoto vrcholku lze usuzovat na dobu

průchodu a průtok

NAP3 Co je třeba si pamatovat

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

NAP3 Co je třeba si pamatovat

dxtEty )()()(Co je konvoluce, dekonvoluce, přenos a korelace

dtetcfc ift2)()(~ 2( ) ( ) iftc t c f e df

Fourierova transformace

Fourierova transformace konvoluce a korelace ( ) ( ) ( )xyC f x f y f )(~)(~)(~ * fyfxfRxy

Diskrétní Fourierova transformace

1 /22 2

0 /2

1 kn knN Ni iN N

n k k nk n N

c c e c c eN

Nyquistova frekvence 1/2

( ) ( ) ( )xyR t x t y d

Jak se změní signál x(t) po průchodu filtrem E(t) Korelace dvou signálů posunutých o čas t