ZÁKLADY TEORIE SIGNÁLŮ A SOUSTAV

VŠB TU Ostrava, Fakulta elektrotechniky a informatiky, Katedra měřící a řídící techniky

ZÁKLADY TEORIE SIGNÁLŮ A SOUSTAV

Pavel Nevřiva

2007

1

PŘEDMLUVA Toto skriptum je věnováno základním metodám, používaným při analýze signálů a soustav. Snahou autora bylo, aby zde výklad aplikačních stránek analýzy signálů a soustav vycházel z názorného teoretického základu.

OBSAH SEZNAM SYMBOLŮ ÚVOD 1 ANALÝZA SIGNÁLŮ SE SPOJITÝM ČASEM 1.1 Základní definice 2 ANALÝZA SIGNÁLU V ČASOVÉ OBLASTI 2.1 Korelační analýza signálu 3 ANALÝZA SIGNÁLU V KMITOČTOVÉ OBLASTI 3.1 Spojité spektrum signálu 3.2 Diskrétní spektrum signálu 4 ANALÝZA SOUSTAV SE SPOJITÝM ČASEM 4.1 Lineární časově invariantní soustava se soustředěnými parametry 5 ANALÝZA LTIL SOUSTAVY V ČASOVÉ OBLASTI 5.1 Časové odezvy LTIL soustavy 5.2 Popis LTIL soustavy diferenciální rovnicí 6 ANALÝZA LTIL SOUSTAVY V KMITOČTOVÉ OBLASTI 6.1 Popis LTIL soustavy kmitočtovým přenosem 6.2 Transformace signálu LTIL soustavou LITERATURA

3 5 9 9 16 16 20 21 40 61 61 74 74 82 99 99 106 106

2

SEZNAM SYMBOLŮ DFT{w[n ]}

{Wk }, diskrétní Fourierova transformace signálu w[n]

DFT-1{Wk }

inversní diskrétní Fourierova transformace k diskrétní Fourierově transformaci signálu w[n] energie signálu w(t ) spojité spektrum energie signálu w(t )

E E (ω) fs F F{w(t )}

frekvence vzorkování soustava, dynamická soustava W (ω) , Fourierova transformace signálu w(t )

F -1{W (ω)} FFT FFT{w[n]}

inversní Fourierova transformace k Fourierově transformaci signálu w(t ) Fast Fourier Transform DFT{w[n ]}, výpočet provedený funkcí fft MATLAB

FFT -1{Wk } fft wFS (t ) g (t ) h (t ) H H (s ) H (ω) i IFFT ifft ± j k LTI LTIL L{w(t )}

DFT-1{Wk }, výpočet provedený funkcí ifft MATLAB funkce fft MATLAB vyjádření signálu w(t ) ve tvaru Fourierovy řady přechodová charakteristika, odezva soustavy na Heavisideův skok impulsová charakteristika, odezva soustavy na Diracův impuls LTIL, lineární časově invariantní soustava se soustředěnými parametry přenos soustavy, L{h(t )} kmitočtový přenos soustavy, F{h(t )} index Inverse Fast Fourier Transform funkce ifft MATLAB −1 index lineární časově invariantní soustava H , lineární časově invariantní soustava se soustředěnými parametry W (s ) , Laplaceova transformace signálu w(t )

L−1{W (s )} n O P P(a, b )

inversní Laplaceova transformace k Laplaceově transformaci signálu w(t ) diskrétní časová proměnná, t = nTs , tzv. diskrétní čas, index výstupní soustava střední výkon signálu w(t ) střední výkon signálu w(t ) za časový interval (a, b )

p (t ) P (ω) PAB PAB (t )

okamžitý výkon signálu w(t ) spojité spektrum výkonu signálu w(t ) střední vzájemný výkon signálů w A (t ) a wB (t ) okamžitý vzájemný výkon signálů w A (t ) a wB (t )

3

{Pm } R (τ) R AB (τ ) s S t T Ts u (t ) w(t ) wrms

spektrum výkonu periodických signálů w(t ) korelační funkce signálu w(t ) vzájemná korelační funkce signálů w A (t ), wB (t ) komplexní proměnná Laplaceovy transformace stav soustavy čas perioda perioda vzorkování vstupní signál, vstup reálný signál se spojitým časem, signál efektivní, RMS, hodnota signálu w(t )

wrms (a, b )

efektivní, RMS, hodnota signálu w(t ) za časový interval (a, b )

ws (t ) w

ideálně vzorkovaný signál w(t ) střední hodnota signálu w(t )

w(a, b )

střední hodnota signálu w(t ) za časový interval (a, b )

W (s ) W (t ) W (ω) {w m }

Laplaceova transformace signálu w(t ) náhodný proces Fourierova transformace signálu w(t ) spektrum periodického signálu w(t )

{w m }

amplitudové spektrum periodického signálu w(t )

y aa (t ) yb (t ) ybb (t ) yC (t ) y N (t ) y P (t ) δ(t ) η(t ) θ(ω) {θ m }

řešení stavové rovnice, stav soustavy výstup nebo celková odezva soustavy, řešení dif. rovnice počáteční podmínky řešení diferenciální rovnice počáteční podmínky řešení diferenční rovnice odezva soustavy při nulovém vstupním signálu přechodná odezva soustavy odezva soustavy při nulovém počátečním stavu ustálená odezva soustavy úplné řešení homogenní rovnice úplné řešení nehomogenní rovnice libovolné partikulární řešení rovnice Diracův impuls (jednotkový impuls) Heavisideova funkce (jednotkový skok) spojité spektrum fáze signálu w(t ) diskrétní spektrum fáze signálu w(t )

(•)∗

komplexně sdružená k (•)

x(t ) y (t ) y (0 − ),y& (0 − ),K y [− 1], y [− 2 ],K ya (t )

4

ÚVOD Změna fyzikální veličiny v čase je studována v mnoha aplikačních úlohách. Často předpokládáme, že v proměnlivosti sledované veličiny je uložena určitá informace. V těchto skriptech nazýváme tuto veličinu signálem. Elektrické napětí, elektrický proud, tlakové impulsy, teplotní změny, jsou nejběžnější signály technické praxe. V těchto skriptech studujeme signály, které jsou popsané svým matematickým modelem. Záleží na konkrétní aplikační oblasti, jaká je fyzikální interpretace určitého matematického modelu signálu. Signály jsou generovány, transformovány a absorbovány fyzikálními elementy, agregáty, strukturami. V těchto skriptech tyto jednotky nazýváme soustavami. Typickým reprezentantem soustav může být skupina RC obvod, tranzistor, zesilovač, mikrofón, vysílač, satelit, přijímač, reproduktor. Jiným reprezentantem může být skupina přístrojů teploměr, tlakoměr, vlhkoměr, anemometr. V těchto skriptech popíšeme soustavu jejím matematickým modelem. Jaká je fyzikální interpretace určitého matematického modelu soustavy záleží opět na konkrétní aplikační oblasti. S pojmy signál a soustava se setkáváme prakticky ve všech oblastech aplikací techniky, počínaje přístroji používanými v domácnostech a konče u nejsložitějších technických zařízení. Vývoj moderní technologie je podložen pokrokem v rozvoji teorie a využití výsledků teorie signálů a soustav. Tato skripta nejsou věnována obecné teorii soustav jako systémů. Je nesporné, že každý systém lze rozložit na subsystémy a že každý systém lze považovat za subsystém určitého většího systému a že složitost této dekomposice a komposice má své meze. Touto filosofickou otázkou se zde nezabýváme. Soustava studovaná v těchto skriptech má jednoduchou strukturu, má jeden vstupní signál u (t ) a jeden výstupní signál y (t ) a je ukázána na obrázku níže. Takovéto schéma kurzu přímo vede k otázkám: Je možné navrhnout a realizovat jednotnou metodu analýzy signálů a soustav za situace, kdy daný signál nebo soustava mohou mít svoji vlastní fyzikální podstatu, princip činnosti, způsob využití? V situaci kdy, například, tisíce 5

fyziků na celém světě pracuje na vývoji různých metod pro vyhodnocení záznamů pozitronové rezonance, magnetické rezonance, CT, EEG? V situaci, kdy počet článků věnovaných automatické regulaci krokových motorů dosahuje několika set položek za rok? V situaci, kdy lze uvést tisíce takovýchto příkladů? Je rozumné studovat obecné přístupy k řešení úlohy analýzy signálů a soustav za situace, kdy je v praxi vývoj nových přístrojů a zařízení v rukou úzce specializovaných odborníků? Kdy je pro dosažení pokroku v aplikaci teorie signálu a soustav potřeba široká znalost a výzkum v dané specializované aplikační oblasti? Čtenář může být ujištěn, že ve skriptech najde základní časem prověřený kurz, věnovaný úvodu do problematiky. Skripta čtenáře seznámí s terminologií a s principy, které jsou společné snad všem otázkám analýzy signálů a soustav. Usnadní mu studium více specializované literatury. Čtenář se seznámí s moderními metodami analýzy signálů a soustav pomocí počítače.

Jednoduchá struktura s LTIL soustavou diskutovanou v těchto skriptech

Základní postupy analýzy signálů a soustav jsou následující. Signál budeme popisovat zejména jeho vhodnými časovými průměry definovanými z časového průběhu signálu. Soustavu budeme popisovat vhodnými vztahy mezi jejím výstupním signálem a počátečním stavem, zejména ale mezi jejím výstupním signálem a vstupním signálem.

6

Analýza v časové oblasti, analýza v kmitočtové oblasti jsou pojmy budoucí diskuse. Teorie analýzy signálů a soustav má dva základní cíle. Základní úlohou analýzy signálu je popsat vlastnosti signálu jeho souhrnnými charakteristikami. Základní úlohou analýzy soustav je určit výstup soustavy buzené jejím počátečním stavem a vstupním signálem. Celkové nebo integrální charakteristiky signálu mohou být časové průměry, integrály a jiné hodnoty vypočítávané nebo měřené z časových průběhů signálů během určitého časového intervalu. Čím více informace tyto charakteristiky obsahují, tím bývají cennější. Spektrum signálu obsahuje více informace o signálu než jeho spektrum výkonu. Spektrum výkonu obsahuje více informace o signálu než jeho střední výkon. Střední výkon signálu může obsahovat více informace o signálu než jeho střední hodnota. Samozřejmě, že nemusí být vždy snadné, dokonce nemusí být ani vždy možné, určit spektrum nebo jiné konkrétní charakteristiky daného signálu. Pokrok v moderní technologii nám umožňuje záznam každého měřeného signálu v prakticky neomezeném počtu bodů. Vzhledem k tomu, že je text skript orientován na úvod do problematiky analýzy signálů a soustav, omezíme diskusi na soustavu s jedním vstupním signálem a jedním výstupním signálem (SISO, single input-single output soustava). Budeme studovat lineární, časově invariantní soustavu, pro kterou, jak uvidíme, platí princip superpozice, studium omezíme na soustavu se soustředěnými parametry (LTIL, Linear Time invariant Lumped System). Takovouto soustavu lze popsat jednoduchým matematickým modelem. Ve skriptech ji popíšeme pomocí diferenciální rovnice. Kurz podává přehled základní metodiky analýzy signálů a soustav. Cílem kurzu je, aby absolvent poznal, co je fyzikální podstatou matematického popisu dynamického chování signálů a soustav v přírodě. Analýza signálů a soustav a její studium jsou podporovány řadou softwarových produktů. Pro tento kurz jsme zvolili MATLAB. Samotný MATLAB má řadu knihoven, nástrojů, 7

toolbox(ů), specializovaných na práci se signály a soustavami. S ohledem na úvodní charakter tohoto kurzu se ve skriptech použijí pouze tři z nástrojů, které MATLAB nabízí. Jsou to jádro MATLAB, na které odkazujeme jako na MATLAB Basic Library, dále MATLAB Simulink Toolbox a MATLAB Symbolic Toolbox. I z těchto tří softwarových produktů jsou v textu použity pouze základní příkazy a funkce. Při použití příkazů a funkcí vyšších úrovní by zápis řešení většiny v textu uvedených příkladů vyšel ještě stručnější. Každý kurz by měl sestávat z teoretického výkladu a z praktických cvičení. Skripta jsou sestavena tak, aby podporovala jak teoretický výklad, tak praktická procvičení látky v simulačním prostředí MATLAB. Simulační úlohy jsou zvoleny tak, aby bylo možno jak úlohy, připravené pro řešení v MATLAB Simulink Toolbox, tak úlohy připravené pro řešení v MATLAB Symbolic Toolbox snadno převést do MATLAB Basic Library.

8

1

ANALÝZA SIGNÁLŮ SE SPOJITÝM ČASEM

Definice 1.0.1

Základní úloha analýzy signálu. Je dán signál w(t ) . Základní úlohou analýzy signálu je popsat vlastnosti signálu jeho souhrnnými charakteristikami.

1.1

Základní definice

Definice 1.1.1

Signál w(t ) se nazývá signál se spojitým časem, pokud časová proměnná t nabývá hodnot z množiny reálných čísel.

Definice 1.1.2

Okamžitá hodnota signálu w(t ) Okamžitá hodnota signálu w(t ) v čase t = t i je w(t i )

Definice 1.1.3

Střední hodnota w signálu w(t ) w

Definice 1.1.4

=

1 T →∞ T lim

+∞ 2 ∫− ∞ w (t ) dt

(1.1.2)

Okamžitý výkon P(t ) signálu w(t )

P(t ) Definice 1.1.6

(1.1.1)

Energie E signálu w(t )

E = Definice 1.1.5

T 2

∫−T 2 w(t ) dt

= w 2 (t )

(1.1.3)

Střední výkon P signálu w(t ) P

=

1 +T/ 2 2 ∫ w (t ) dt T → ∞ T − T/ 2 lim

(1.1.4)

9

Definice 1.1.7

Okamžitý vzájemný výkon PAB (t ) signálů w A (t ) a wB (t ) PAB (t )

Definice 1.1.8

= w A (t )wB (t )

(1.1.5)

Střední vzájemný výkon PAB signálů w A (t ) a wB (t ) PAB

=

1 +T/ 2 ∫ w A (t )wB (t ) dt T → ∞ T − T/ 2 lim

(1.1.6)

Příklad 1.1.1

Výpočet výkonu Uvažujme velmi jednoduchou soustavu tvořenou jediným elektrickým odporem. Soustava je popsána pouze jedním parametrem, kterým je velikost R elektrického odporu. Vstupním signálem u (t ) soustavy nechť je napětí na odporu v(t ) = V cos(ω0 t ) . Výstupním signálem y (t ) soustavy nechť je proud odporem i (t ) = I cos(ω0 t ) .

i (t ) =

u (t ) R

Poznamenejme, že v tomto případě je vstupní a výstupní signál soustavy ve fázi, mezi vstupním a výstupním signálem není fázový posun.

Obrázek 1.1.1

Příklad 1.1.1. Elektrický odpor. Fyzikální model a schéma soustavy.

Střední hodnota časového průběhu napěťového vstupního signálu je

10

v

1 T 2 = lim ∫ v(t ) dt T → ∞ T −T 2 1 T 2 = lim ∫ V cos (ω0t )dt T → ∞ T −T 2

=

1 T0

(1.1.7)

T0 2 V cos (ω0t )dt = 0 2

∫−T0

kde T0 = 2π / ω0 . Obdobně, y ≡ i = 0 . Tento výsledek ovšem známe předem i bez výpočtu, každý ví, že střední hodnota kosinusovky je rovna nule. Okamžitý výkon vstupního signálu v(t ) = V cos(ω0 t ) je

1 Pv (t ) = (V cos(ω0 t ) )2 = V 2 (1 + cos(2ω0 t ) ) 2

(1.1.8)

Poznamenejme, že okamžitý výkon napětí na odporu měřený ve fyzikálních jednotkách (watt) je P (t ) p (t ) = v R

(1.1.9)

Obdobně, Pi (t ) = (I cos(ω0 t ) )2 =

1 2 I (1 + cos(2ω0 t ) ) 2

(1.1.10)

Poznamenejme, že okamžitý výkon proudu protékajícího odporem měřený ve fyzikálních jednotkách je p(t ) = RPi (t )

(1.1.11)

Střední výkon vstupního signálu v(t ) = V cos(ω0 t ) je Pv

1 +T/ 2 2 ∫ v (t ) dt T →∞ T −T/ 2 1 +T0 / 2 (V cos( ω0t ))2 dt = ∫ 2 − T / T0 0 =

lim

V 2 1 +T0 / 2 (1 + cos( 2ω0t )) dt = 2 T0 ∫−T0 / 2 =

(1.2.12)

V2 2 11

Poznamenejme, že střední výkon napětí na odporu měřený ve fyzikálních jednotkách je P p= v R

Obdobně,

Pi =

I2 2

(1.1.13)

Poznamenejme, že střední výkon proudu protékajícího odporem měřený ve fyzikálních jednotkách je p = RPi Okamžitý vzájemný výkon mezi vstupním signálem u (t ) ≡ v(t ) = V cos(ω0 t ) a výstupním signálem y (t ) ≡ i (t ) = I cos(ω0 t ) odporu je 1 Pvi (t ) = (V cos(ω0 t ) )(I cos(ω0 t ) ) = VI ( 1 + cos(2ω0 t ) ) 2

(1.1.14)

Poznamenejme, že okamžitý vzájemný výkon mezi vstupním signálem u (t ) ≡ v(t ) = V cos(ω0 t ) a výstupním signálem y (t ) ≡ i (t ) = I cos(ω0 t ) na odporu se

1 měřením ve fyzikálních jednotkách nemění, pvi (t ) = VI (1 + cos(2ω0t ) ) . 2 Střední vzájemný výkon mezi vstupním signálem u (t ) ≡ v(t ) = V cos(ω0 t ) a výstupním signálem y (t ) ≡ i (t ) = I cos(ω0 t ) odporu je Pvi

1 +T/ 2 ∫ v(t )i(t ) dt T → ∞ T − T/ 2 1 +T0 / 2 (V cos(ω0t ) )(I cos(ω0t ) ) dt = T0 ∫− T0 / 2 VI +T0 / 2 (1 + cos(2ω0t ) ) dt = 2T0 ∫−T0 / 2 VI = 2 =

lim

(1.1.15)

Poznamenejme, že okamžitý vzájemný výkon mezi vstupním signálem u (t ) ≡ v(t ) = V cos(ω0 t ) a výstupním signálem y (t ) ≡ i (t ) = I cos(ω0 t ) odporu se měřením ve fyzikálních jednotkách nemění, p =

VI . 2

Obrázek 1.1.2 ukazuje časové průběhy signálů v(t ) = V cos(ω0 t ) a i (t ) = I cos(ω0 t ) , konkrétně signálů v(t ) = 4 cos(1t ) a i (t ) = 2 cos(1t ) a časový průběh jejich okamžitého

12

vzájemného výkonu. Obrázek 1.1.2 byl generovaný sekvencí příkazů MATLAB

t = -10:0.001:10; v = 4*cos(1*t); i = 2*cos(1*t); Pvi = v.*i; subplot(311) plot(t,v) ylabel('v(t)') grid subplot(312) plot(t,i) ylabel('i(t)') axis([-10 10 -3 3]) grid subplot(313) plot(t,Pvi) ylabel('Pvi(t)') grid xlabel('t')

Obrázek 1.1.2

Příklad 1.1.1. Výpočet okamžitého vzájemného výkonu.

13

Definice 1.1.9

Efektivní hodnota wef signálu w(t ) wef

=

1 +T/ 2 2 ∫ w (t ) dt T → ∞ T − T/ 2 lim

(1.1.16)

= P

Definice 1.1.10 Energetický signál w(t )

Energetický signál w(t ) splňuje podmínky E ≠ 0 a E < ∞ , kde E je energie signálu E =

∞ 2 ∫− ∞ w (t ) dt

(1.1.17)

Definice 1.1.11 Výkonový signál w(t )

Výkonový signál w(t ) splňuje podmínky P ≠ 0 and P < ∞ , kde P je střední výkon signálu

1 T 2 w 2 (t ) dt ∫ − T 2 T →∞ T

P = lim

(1.1.18)

Definice 1.1.12 Diracův impuls δ( x ) je distribuce, která splňuje následující podmínky

δ( x ) = 0 x ≤ 0 −

δ( x ) = 0 x ≥ 0 +

(1.1.19)

∞

∫− ∞ δ(x − x1 ) dx = 1 Platí ∞ ∫− ∞ f (x )δ(x − x1 )dx = f (x1 )

(1.1.20)

∞ ∫− ∞ f (x + x1 )δ(x )dx = f (x1 )

(1.1.21)

též

kde f ( x ) je jednoznačná ohraničená funkce, definovaná v x1

14

Decibel

Pojem decibel byl původně definovaný jako jednotka zesílení výkonu v elektrických obvodech. Nejdříve, tato definice již zastarala, byla jako jednotka výkonového zesílení přijata jednotka bell, nazvaná podle Alexandera Grahama Bella, 1847-1922, amerického (narozeného ve Skotsku) vynálezce telefonu. Logaritmická stupnice jednotky decibel umožňuje pracovat s hodnotami fyzikálních veličin nacházejícími se ve velikých rozsazích klasických lineárních stupnic a kompenzuje exponenciální charakteristiky mnoha fyzikálních a fyziologických zákonů a vztahů. Py ⎛ Py ⎞ ⎜ ⎟ = log10 ⎜P ⎟ Pu ⎝ u ⎠ bell kde

(1.1.22)

u (t ) je vstupní signál y (t ) je výstupní signál

Postupem doby bylo zjištěno, že tato jednotka je příliš velká pro praktické používání a proto byla definována jednotka desetkrát menší, jeden decibel (označená jako dB). Konkrétně, výkonové zesílení obvodu v decibelech je

⎛ Py ⎜ ⎜P ⎝ u

P ⎞ ⎟ = 10 log10 y ⎟ Pu ⎠ dB

(1.1.23)

Střední výkon signálu je druhá mocnina efektivní hodnoty signálu. V decibelech tedy můžeme vyjádřit poměr výkonu P signálu w(t ) k výkonu N šumu n(t ) , měřených v některém bodu soustavy, jako wef P ⎛P⎞ = 20 log10 ⎜ ⎟ = 10 log10 N nef ⎝ N ⎠ dB

(1.1.24)

Poznamenejme, že míra decibel má univerzální použití. Míra decibel se často používá pro vyjádření poměru výkonu, respektive hladiny výkonu, vzhledem k určité referenční hladině.

15

Příklad 1.1.2

Decibely výkonu vzhledem k hladině výkonu 1 mW Decibely výkonu vzhledem k hladině výkonu 1 mW jsou definovány vztahem

⎛ skutecna hladina vykonu (watt) ⎞ dBm = 10 log10 ⎜ ⎟ 10 -3 ⎝ ⎠ = 30 + 10 log10 [ skutecna hladina vykonu (watt)]

(1.1.25)

kde "m" v dBm označuje miliwattovou referenční hladinu. Laboratorní rádiové frekvenční signálové generátory jsou obvykle kalibrovány v dBm.

2

ANALÝZA SIGNÁLU V ČASOVÉ OBLASTI

2.1 Korelační analýza signálu Definice 2.1.1

Korelační funkce R (τ) energetického signálu w(t ) R(τ ) = ∫

∞

−∞

Definice 2.1.2

w(t ) w(t + τ ) dt

(2.1.1)

Korelační funkce R (τ) výkonového signálu w(t ) 1 T 2 ∫ w(t ) w(t + τ)dt T →∞ T −T 2

R (τ) = lim

(2.1.2)

Příklad 2.1.1

Korelační funkce R (τ) výkonového signálu w(t ) Uvažujme signál w(t ) = w1 (t ) + w2 (t ) . w1 (t ) je náhodný šum, generovaný příkazem MATLAB w1=4*randn(1,100001);

16

w2 (t ) je harmonický signál, generovaný sekvencí příkazů MATLAB t=0:1:100000; w2=1.414*cos(.2*t);

Signál w(t ) je generovaný příkazem MATLAB w = w1+w2;

Signály w1(t ) , w2 (t ) , w(t ) jsou vykresleny sekvencí příkazů MATLAB subplot(511) plot(t,w1) ylabel('w1(t)') %xlabel ('t') axis([0 1000 min(w1+w2) max(w1+w2)]) subplot(512) plot(t,w2) ylabel('w2(t)') %xlabel ('t') axis([0 1000 min(w1+w2) max(w1+w2)]) w=w1+w2; subplot(513) plot(t,w) ylabel('w(t)=w1(t)+w2(t)') %xlabel ('t') axis([0 1000 min(w1+w2) max(w1+w2)])

Signály w1(t ) , w2 (t ) , w(t ) ukazuje Obrázek 1. 3 A,B,C. Sekvence příkazů MATLAB [R,tau]=xcorr(w,50000,'unbiased'); subplot(514) plot(tau,R) ylabel('R(tau)') %xlabel('tau') axis([-500 500 min(R)-5 max(R)+5])

vypočítává a zobrazuje korelační funkci R (τ) výkonového signálu w(t ) , zobrazenou na Obrázku 2.1.1D. Z korelační funkce R (τ) lze detekovat harmonickou složku w2 (t ) signálu w(t ) , která je v šumu w1 (t ) skrytá a není z časového průběhu signálu w(t ) ihned patrná.

17

Obrázek 2.1.1

Příklad 2.1.1. Korelační funkce R (τ) výkonového signálu w(t ) = w1 (t ) + w2 (t )

Sekvence příkazů MATLAB subplot(515) plot(tau,R) ylabel('R(tau)') %xlabel('tau') axis([-500 500 -2 2])

vykresluje korelační funkci v mezích, které nejsou závislé na velikosti náhodného šumu.

Definice 2.1.3

Vzájemná korelační funkce R AB (τ ) energetických signálů wA (t ), wB (t )

R AB (τ ) = ∫

∞

w (t ) wB (t + τ ) dt −∞ A

(2.1.3)

18

Příklad 2.1.2

Vzájemná korelační funkce R AB (τ ) energetických signálů wA (t ), wB (t ) Uvažujme energetické signály w A (t ) a wB (t ) generované sekvencí příkazů MATLAB t=0:1:1000; wA=sign(t-100)-sign(t-200); wB= cos(0.1*t).*(sign(t-300)-sign(t-400));

Signály w A (t ) a wB (t ) jsou zobrazeny na Obrázku 2.1.2 A, B. Vzájemná korelační funkce R AB (τ ) energetických signálů w A (t ) a wB (t ) vypočtená příkazem MATLAB [RAB, tau]=xcorr (wB,wA,'none'); je zobrazená na Obrázku 2.1.2 C. Z průběhu vzájemné korelační funkce R AB (τ ) lze detekovat časovou závislost signálů w A (t ) a wB (t ) . Obrázek 2.1.2 byl generován následující sekvencí MATLAB

subplot(311) plot(t,wA) ylabel('wA(s)') xlabel('t(s)') axis([0 1000 min(wA)-1 max(wA)+1]) grid subplot(312) plot(t,wB) ylabel('wB(s)') xlabel('t(s)') axis([0 1000 min(wB)-1 max(wB)+1]) grid subplot(313) plot(tau,RAB) ylabel('R(tau)') xlabel('tau(s)') axis([-500 500 min(RAB)-10 max(RAB)+10]) grid

19

wA(t)

3 2 1 0 -1

0

100

200

300

400

500 t(s)

600

700

800

900

1000

0

100

200

300

400

500 t(s)

600

700

800

900

1000

-400

-300

-200

-100

0 tau(s)

100

200

300

400

500

wB(t)

2 0

R(tau)

-2

Obrázek 2.1.2

Definice 2.1.4

80 60 40 20 0 -20 -500

Příklad 2.1.2. Vzájemná korelační funkce R AB (τ ) energetických signálů wA (t ), wB (t )

Vzájemná korelační funkce R AB (τ ) výkonových signálů w A (t ), wB (t ) 1 T 2 ∫ wA (t ) wB (t + τ)dt T →∞ T −T 2

R AB (τ) = lim

3.

(2.1.4)

ANALÝZA SIGNÁLU V KMITOČTOVÉ OBLASTI

V tomto kurzu se nezabýváme náhodnými signály, šumy. Zbývající signály jsme rozdělili na signály energetické a na signály výkonové periodické. U obou typů signálů budeme předpokládat, že splňují Dirichletovy podmínky, viz dále, pro práci metodami Fourierovy analýzy. Energetický signál vyjádříme pomocí Fourierovy transformace, popíšeme ho pomocí jeho spojitého spektra.

20

Výkonový periodický signál vyjádříme Fourierovou řadou, popíšeme ho pomocí jeho diskrétního spektra. Poznámka

Harmonický signál w(t ) je daný rovnicí

w(t ) = wmax cos(ωt + θ)

(3.0.0)

kde wmax je konstantní parametr, amplituda harmonického signálu ωt + θ je fáze harmonického signálu ω je konstantní parameter, úhlový kmitočet harmonického signálu θ je konstantní parametr, počáteční fáze harmonického signálu t je nezávisle proměnná harmonického signálu, čas

Poznámka:

Úhlový kmitočet ω má rozměr t −1 . Jeho fyzikální jednotka je

sekunda na mínus prvou, s −1 . Název jeho fyzikální jednotky je radián za sekundu,

rad.s−1 .

Kmitočet f má rozměr t −1 . Jeho fyzikální jednotka je sekunda

na mínus prvou, s −1 . Název jeho fyzikální jednotky je hertz.

3.1

Spojité spektrum signálu

Postačující, ne však nutnou, podmínkou, aby byl signál analyzovatelný pomocí Fourierovy transformace, tj popsán svým spojitým spektrem je, aby signál splňoval v nekonečném časovém intervalu tzv. Dirichletovy podmínky dané definicí Definice 3.1.1. Uvidíme dále, že tyto podmínky jsou blízké Dirichletovým podmínkám požadovaným pro konvergenci Fourierových řad, popisujících diskrétní spektrum výkonových periodických signálů.

Definice 3.1.1

Postačující Dirichletovy podmínky, aby funkce y (t ) měla Fourierovu transformaci jsou následující 1.

y (t ) musí být absolutně integrovatelná na intervalu (− ∞, ∞ ) , tj.

21

∞

∫−∞ y(t ) dt ≤ K < ∞

(3.1.1)

kde K je konečná konstanta y (t ) může mít na libovolném konečném časovém intervalu

2.

pouze konečný počet nespojitostí prvého řádu a konečný počet maxim a minim. Připomeňme, že má-li funkce v určitém bodě nespojitost prvého řádu, znamená to, že zde má jak konečnou limitu zleva, tak konečnou limitu zprava, a tyto limity si nejsou rovny.

Věta 3.1.1

Časový průběh libovolného fyzikálního, tj. fyzikálně realizovatelného signálu (tj. signálu s konečnou energií) může být popsán spojitou reálnou funkcí, která splňuje Dirichletovy podmínky.

Z uvedeného vyplývá, že fyzikání signály mají teoreticky vždy své Fourierovy transformace. Nevyplývá z toho ale závěr, že je dokážeme vždy určit.

V současné literatuře se bez výjimky kosinový průběh signálů vyjadřuje pomocí součtu komplexních exponenciál. Připomeňme si proto úvodem, jak se vyjádří harmonický signál pomocí komplexní proměnné. Pochopení tohoto zápisu je nezbytné pro pochopení celé další látky

Jde přitom pouze o geometrickou a fyzikální interpretaci Eulerova vztahu

cos x =

e jx + e − jx 2

(3.1.2)

22

Im j

j*sin(x)

e jx 1

x 0

Re 2*cos(x)

cos(x)

Obrázek 3.1.1 A cos x =

e jx + e − jx 2

Im jA

Ae

θ

0 ω

j (ωt + θ )

A

t=0

jA sin (ωt + θ) t = 0

Re 2 A cos(ωt + θ) t = 0

A cos (ωt + θ ) t = 0

Obrázek 3.1.1 B A cos (ωt + θ ) =

Ae j (ωt + θ ) + Ae − j (ωt + θ ) 2

23

Im

A j (ωt + θ ) e t=0 2 ω

θ

0

Re

A cos(ωt + θ ) t = 0

A − j (ωt + θ ) e t=0 2

Obrázek 3.1.1 C A cos (ωt + θ) =

A j (ωt + θ ) A − j (ωt + θ ) e + e 2 2

Im w m e jθ m e jω m t t = 0

ωm

0

θm

Re

wm, max cos (ω m t + θ m ) t =

w m e − jθ m e − jω m t t = 0

Obrázek 3.1.1 D wm (t ) = wm, max cos (ω m t + θ m )

24

wm (t ) = wm ,max cos(ω m t + θ m ) wm ,max j ω m t + θ m wm ,max = e + 2 2 wm ,max jθ w m ,max = e m e jω m t + 2 2

(

=

w m e jθ m e jω m t

= w m e jω m t

)

e − j (ω m t + θ m )

(3.1.3)

e − jθ m e − jω m t

+ w m e − jθ m e − jω m t + w − m e − jω m t

Poznámka: Uveďme ještě některé užitečné vztahy pro manipulaci s harmonickým signálem v reálné oblasti. V reálné oblasti harmonický signál zapisujeme obvykle některým z následujících zápisů:

wm (t ) = =

wm, max cos(ωmt + θ m )

(3.1.4)

wm, max (cos ωmt cos θ m − sin ωm t sin θ m )

= am cos ωmt + bm sin ωm t ze vztahu vidíme, že lze harmonický signál s fázovým posunem složit z neposunutých kosinusového a sinusového signálu. Pro jejich amplitudy snadno odvodíme am = wm,max cos θm bm = − wm,max sin θm

(3.1.5)

a naopak 2 + b2 wm, max = am m b θ m = − tan −1 m am

Definice 3.1.2

(3.1.7)

Harmonický element dwi (t ) dwi (t ) =

Věta 3.1.2

(3.1.6)

1 W (ωi ) cos (ωi t + θ i ) dω π

(3.1.8)

Každý fyzikální signál w(t ) (tj. signál s konečnou energií) může být vyjádřen inversní Fourierovou transformací w(t ) =

1 ∞ W (ω)e jωt dω = F −1 {W (ω )} 2 π ∫− ∞

(3.1.9)

t ∈ (− ∞ , ∞ )

25

kde komplexní Fourierova transformace W (ω) signálu w(t ) je dána spojitým spektrem W (ω ) =

∞

∫− ∞ w(t )e

− jωt

dt = F {w(t ) }

(3.1.10)

ω ∈ (− ∞ , ∞ )

W (ω) je komplexní funkce reálné proměnné ω .

Poznámka: Uveďme souvislost (3.1.9) a (3.1.10) se zápisem Fourierovy transformace v reálném oboru.

Ze základního kurzu matematiky je známo, že každý signál w(t ) pro který platí předchozí podmínky, může být vyjádřen inversní Fourierovou transformací 1 ∞ [ a(ω)cos ωt +b(ω)sin ωt ] dω 2π ∫0 1 ∞ = A(ω)cos(ωt + θ(ω)) dt 2π ∫0

w(t ) =

(3.1.11)

t ∈ (− ∞, ∞ )

kde

a(ω) = A(ω)cos θ(ω) b(ω) = − A(ω)sin θ(ω)

A(ω) = a 2 (ω) + b 2 (ω) ; A ≥ 0 b(ω) θ(ω) = − tan −1 ; θ(ω) ∈ (− π;+ π] a(ω) ω ∈ [0, ∞ )

(3.1.12)

(3.1.13) (3.1.14)

Pro reálné koeficienty a (ω) a b(ω) reálné proměnné ω platí vztahy Fourierovy transformace ∞ w(t )cos ωt dt −∞ ∞ b(ω) = 2 w(t )sin ωt dt −∞

a (ω) = 2 ∫

(3.1.15)

∫

(3.1.16)

A(ω) = wmax (ω) je jednostranné spojité spektrum amplitudy signálu w(t ) . θ(ω) je jednostranné spojité spektrum fáze signálu w(t ) . Jednostranné spojité spektrum signálu w(t ) ukazuje, že na kmitočtu ωi ≥ 0 signál přispívá ke svému celkovému průběhu přírůstkem

26

1 1 A(ωi )cos(ωit + θ(ωi )) dω = wmax (ωi )cos(ωit + θ(ωi )) dω 2π 2π

(3.1.11) (3.1.17)

Pomocí Eulerova vztahu vzorce (3.1.2)a vztahů (3.1.15), (3.1.16) snadno rovnici (3.1.11) upravíme do dnes již výhradně používaných tvarů (3.1.9) a poté i odvodíme vztah pro (3.1.10): e j ω t + e − jω t a (ω)cos ωt + b(ω)sin ωt = a (ω) 2 e jωt + e − jωt 2 a (ω) − jb(ω) jωt = e 2 = a (ω)

Označíme W (ω) =

e j ωt − e − j ωt + b(ω) 2j − +

e j ω t − e − jω t 2 a (ω) + jb(ω) − jωt e 2

jb(ω)

a(ω) − jb(ω) 2

(3.1.18)

(3.1.19)

Podle (3.1.12) je W (− ω) =

a(ω) + jb(ω) 2

a W (− ω) = W * (ω)

(3.1.20)

a(ω)cos ωt + b(ω)sin ωt = W (ω)e jωt + W (− ω)e − jωt

(3.1.21)

Získáme

Integrál (3.1.11) nyní jen upravíme: w(t ) = = = = = =

1 ∞ A(ω) cos(ωt + θ(ω)) dω 2π ∫0 1 ∞ [a(ω) cos ωt + b(ω) sin ωt ] dω 2π ∫0 1 ∞ W (ω)e jωt + W (− ω)e − jωt dω ∫ 0 2π 1 ∞ 1 ∞ W (ω)e jωt dω + ∫ W (− ω)e − jωt dω ∫ π 0 2π 0 1 ∞ 1 0 W (ω)e jωt dω + ∫ W (ω)e jωt dω ∫ π −∞ 2π 0 1 ∞ W (ω)e jωt dω 2π ∫− ∞

[

]

(3.1.22)

t ∈ (− ∞ , ∞ )

27

Využili jsme skutečnosti, že pro W (ω) platí ∞

∫

0 ω j t W (ω)e dω = W (− ω)e − jωt dω

∫

(3.1.23)

−∞

0

Konečně odvodíme vztah pro W (ω) : W (ω) =

a(ω) − jb(ω) 2 ∞

∞

−∞ ∞

−∞

∫ w(t )cos ωt dt − j ∫ w(t )sin ωt dt

=

(3.1.24)

∫ w(t )(cos ωt − j sin ωt )dt

=

−∞ ∞

∫ w(t )e

=

− jωt dt

−∞

ω ∈ (− ∞, ∞ )

Získáváme klasický tvar Fourierovy transformace (3.1.10) a odpovídající klasický tvar zpětné Fourierovy transformace (3.1.9):

W (ω) =

∞

∫− ∞ w(t )e

− jωt dt

= F {w(t ) }

(3.1.10)

(3.1.27)

(3.1.25)

( 3.1.09)

(3.1.26)

ω ∈ (− ∞ , ∞ )

1 ∞ W (ω)e jωt dω = F −1 {W (ω )} 2 π ∫− ∞ t ∈ (− ∞ , ∞ )

w(t ) =

Definice 3.1.3

Fourierova transformace W (ω ) signálu w(t ) se nazývá (dvoustranné spojité) spektrum W (ω ) signálu w(t )

W (ω) = ∫

∞

−∞

w(t ) e − jωt dt

(3.1.27)

ω ∈ (− ∞ , ∞ )

28

Definice 3.1.4

Dvoustranné spojité spektrum amplitudy W (ω) signálu w(t ) W (ω) = Re 2 {W (ω)}+ Im 2 {W (ω)}

Definice 3.1.5

(3.1.28)

Dvoustranné spojité spektrum fáze θ(ω) signálu w(t ) θ(ω) = arg(W (ω))

(3.1.29)

Platí: W (ω) = W (− ω) θ(ω) = −θ(− ω)

Dvoustranné spojité spektrum ukazuje, že na fyzikálním kmitočtu ωi , fyzikální kmitočet ωi je vždy nezáporný, signál přispívá ke svému celkovému průběhu přírůstkem 1 1 W (− ωi )e − jωi t dω + W (ωi )e jωi t dω = 2π 2π 1 wmax (ωi ) − jθ(ωi ) − jωi t 1 wmax (ωi ) jθ(ωi ) jωit e e dω + e e dω = 2π 2 2π 2 ⎡ e − jθ(ωi )e − jωi t + e jθ(ωi )e jωit ⎤ 1 ⎥ dω = wmax (ωi )⎢ 2π 2 ⎢⎣ ⎥⎦ 1 wmax (ωi )cos(ωit + θ(ωi )) dω 2π tedy stejným přírůstkem jako v případě jednostranného spojitého spektra. To je ovšem očekávaný výsledek. Signál se nemění, zavedl se pouze formálně jiný popis signálu. Definice 3.1.6

Dvoustranné spojité spektrum energie E (ω) signálu w(t )

E (ω) = W (ω)

Definice 3.1.7

2

(3.1.30)

Dvoustranné spojité spektrum výkonu P(ω) signálu w(t ) , tzv. spektrální hustota výkonu P(ω) signálu w(t )

P(ω) = lim

WT (ω)

T E (ω) = lim T T T →∞ T →∞

2

(3.1.31)

kde WT (ω) je spojité spektrum energie signálu w(t ) za časový 2

interval T . 29

Dále již pracujeme, až na ilustrační Poznámku u výkladu diskrétních spekter, jen s dvoustrannými spektry, a tuto skutečnost proto nezdůrazňujeme. Tam kde je to zřejmé ze symboliky také neuvádíme, že se jedná o spektra spojitá.

Příklad 3.1.1

Fourierova transformace pravoúhlého impulsu ⎛t⎞ Je dána kladná konstanta τ , symbolem Π⎜ ⎟ označujeme jednotkový pravoúhlý ⎝τ⎠ impuls o délce trvání τ sekund, definovaný τ τ ⎛t⎞ Π⎜ ⎟ = 1 , − ≤ t < 2 2 ⎝τ⎠ ⎛t⎞ Π⎜ ⎟ = 0 jinde ⎝τ⎠ ⎛t⎞ Pravoúhlý impuls w(t ) = Π ⎜ ⎟ je nakreslený na Obrázku 3.1.2 . ⎝τ⎠

(3.1.32)

⎛t⎞ Příklad 3.1.1 Pravoúhlý impuls w(t ) = Π ⎜ ⎟ ⎝τ⎠

Obrázek 3.1.2

Impuls má Fourierovu transformaci W (ω ) = =

∞

∫− ∞ w(t )e τ2

∫− τ 2 1e

− jωt

− jωt

dt

dt

e − jω τ 2 e + jω τ 2 = − − jω − jω 2 ⎛ ωτ ⎞ = sin ⎜ ⎟ ω ⎝ 2 ⎠

(3.1.33)

30

Vyjádříme W (ω) pomocí funkce sinc (λ ) definované vztahem (3.1.34) sinc(λ ) =

sin (πλ ) πλ

(3.1.34)

Toto vyjádření dává ⎛ τω ⎞ W (ω) = τ sinc ⎜ ⎟ ⎝ 2π ⎠

(3.1.35)

⎛t⎞ Pro w(t ) = Π ⎜ ⎟ získáme spojité spektrum W (ω) uvedené na Obr. 1.6. ⎝5⎠ Na Obr. 3.1.3 je naposledy, pro srovnání, uvedeno i jednostranné spojité spektrum tohoto signálu.

Obrázek 3.1.3

⎛t⎞ Dvoustranné a jednostranné a spojité spektrum pravoúhlého impulsu Π⎜ ⎟ ⎝5⎠

31

Příklad 3.1.2

Fourierova transformace exponenciálního signálu Nyní uvažujme signál w(t ) = 1e

−2 t

(3.1.36)

Analytický tvar Fourierovy transformace lze opět najít mnoha způsoby. Uveďme zde kompletní řešení, dané programem v prostředí MATLAB Symbolic Toolbox.

syms t omega; w = exp(-2*abs(t)); subplot (211) ezplot(w,[-3,3]); axis([-3 3 0 1]); grid W = Fourier(w,t,omega); subplot(212) ezplot(W,[-3,3]); grid

Obrázek 3.1.4

% variables % signal

% Fourier transform

Příklad 3.1.2. Graf exponenciálního signálu (A) a jeho Fourierovy transformace určené pomocí MATLAB Symbolic Toolbox (B)

32

Příklad 3.1.3

Numerický výpočet Fourierovy transformace

Grafický průběh Fourierovy transformace konkrétního signálu snadno určíme numericky. Prostudujte si dále uvedený jeden takovýto jednoduchý program pro MATLAB. Je v něm zapsáno řešení příkladu 3.1.1. Navíc je v něm uvedena i zpětná Fourierova transformace, tj. rekonstrukce signálu z jeho spojitého spektra. Ověřte si tato řešení a zkuste si též vypočítat spojitá spektra a rekonstrukce i jiných signálů. %Na intervalu promenne t=(-50,50) se studuje %signal w(t)o trvani 100s. Obsahuje jednotkovy sudy impuls %o delce 5s a vysce 1. %Hleda se FT W(omega) signalu w(t). %W(omega) se urci prostym vypoctem tak, %ze se pro integraci signal w(t) vyjadri ve 1001 bodech hodnotou %vektoru w a vysledek integrace se ulozi do 1001 bodu vektoru W. %Zvolime, ze sojite spektrum budeme hledat na intervalu pro %omega vetsi nez -10s-1 a mensi nez 10s-1. %poloha clenu ve vektoru w: %w(1)=w(t=-50), w(501)=w(t=0), w(1001)=w(t=50) %casova vzdalenost jednotlivych clenu ve vektoru w: %delta t = dt=0.1s %vztah mezi infexem nn clenu w(nn) ve vektoru w a casem t %t=(nn-501)/10.s %poloha clenu ve vektoru W: %W(1)=W(omega=-10s-1), % W(501)=W(omega=0s-1), % W(1001)=W(omega=10s-1) %kmitoctova vzdalenost jednotlivych clenu ve vektoru W: %delta omega = domega=0.02s-1 %vztah mezi indexem kk clenu W(kk) ve vektoru W a kmitoctem omega: %omega=(kk-501)/50.s-1

33

%Program pro vypocet FT: w=[zeros(1,475) ones(1,51) zeros(1,475)]; %generuje se signal w(t) t=-50:0.1:50; subplot(311); plot(t,w) %vykresluje se signal w(t) grid; axis([-50 50 -1 2]); xlabel('t'); ylabel('w(t)'); for kk=1:1:1001; W(kk)=0; %nuluje se vektor W, kam se ulozi FT end; for kk=1:1:1001; for nn=1:1:1001; W(kk)=W(kk)+w(nn)*exp(-j*(kk-501)/50*(nn-501)/10)*0.1; %FT end; end; omega=-10:0.02:10; subplot(312); plot(omega,abs(W)) %vykresluje se W(omega) grid xlabel('omega'); ylabel('abs(W(omega))')

%Rekonstrukce signalu z jeho FT %Pro rekonstrukci signalu slouzi jeho vyse spoctena FT %Program pro vypocet zpetne FT for nn=1:1:1001; wrek(nn)=0; %nuluje se vektor wrek, kam se ulozi rekonstruovany signal end; for nn=1:1:1001; for kk=1:1:1001; wrek(nn)=wrek(nn)+1/(2*pi)*W(kk)*exp(j*(nn-501)/10*(kk501)/50)*0.02;%zpetna FT end; end; subplot(313); plot(t,abs(wrek)); %Vykresluje se rekonstruovany signal. %Vypocital se jen z nizkych frekvenci. %Hrany rekonstruovaneho signalu wrek(t) jsou proto oproti w(t) zkreslene. grid axis([-50 50 -1 2]); xlabel('t'); ylabel('wrek(t)');

34

Obrázek 3.1.5

Příklad 3.1.3. Numerický výpočet přímé a zpětné Fourierovy transformace. ⎛t⎞ Jednotkový impuls Π⎜ ⎟ (A), jeho amplitudové spojite spektrum (B) ⎝5⎠ a zpětná rekonstrukce (C).

Vlastnosti Fourierovy transformace

Uvažujme signály w(t ) a u (t ) které mají Fourierovy transformace W (ω ) = F{w(t )}

(3.1.37)

U (ω ) = F{u (t )}

(3.1.38)

Ze základních vlastností Fourierovy transformace uveďme tyto:

35

Linearita

Fourierova transformace je lineární, protože F{a u (t ) + b w (t )} =

∞

∫ [a u (t ) + b w (t )]e

− jωt

dt

−∞ ∞

∞

= a ∫ u (t ) e − jωt dt + b ∫ w (t ) e − jωt dt −∞

(3.1.39)

−∞

= aF{ u (t )} + bF{ w (t )} Posun v časové oblasti F{w(t − t0 )} = W (ω )e − jωt0

(3.1.40)

Posun v kmitočtové oblasti

{

} (

F w(t )e jω0t = W ω − ω 0

)

(3.1.41)

Dualita

f (− ω ) =

1 F{W (t )} 2π

(3.1.41)

Změna časového měřítka

F{w(a t )} =

1 ⎛ω ⎞ W⎜ ⎟ , a ≠ 0 a ⎝a⎠

( 3.1.42)

Derivování v časové oblasti

⎧⎪ d n w( t ) ⎫⎪ F⎨ = ( jω )n W (ω ) , n = 1, 2, L n ⎬ ⎪⎩ dt ⎪⎭

(3.1.43)

Integrování v časové oblasti t 1 F⎧⎨∫ w(λ )dλ ⎫⎬ = W (ω ) + π W (0)δ (ω ) ⎩ −∞ ⎭ jω

(3.1.44)

36

Příklad 3.1.4

Fourierova transformace trojúhelníkového impulsu

Uvažujme trojúhelníkový impuls w(t ) uvedený na Obrázku 3.1.6 A. Nejdříve najdeme derivaci w& (t ) signálu w(t ) . Derivace w& (t ) je zobrazena na Obrázku 3.1.6 B. Vidíme, že w& (t ) =

Obrázek 3.1.6

2 ⎛t +T 4⎞ 2 ⎛t −T 4⎞ ⎟ ⎟ − Π⎜ Π⎜ T ⎜⎝ T 2 ⎟⎠ T ⎜⎝ T 2 ⎟⎠

(3.1.45)

Příklad 3.1.4. Signál w(t ) a jeho derivace w& (t )

Fourierova transformace W& (ω ) signálu w& (t ) může být určena s použitím transformačních dvojic z Tabulky 1. 1 v kombinaci s větou o posunu v časové oblasti. Získáme

37

Tω ⎞ ⎡ ⎛ jTω ⎞ ⎛ ⎛ jTω ⎞⎤ W& (ω ) = ⎜ sinc ⎟ ⎢exp⎜ ⎟ − exp⎜ − ⎟ 4π ⎠ ⎣ ⎝ 4 ⎠ 4 ⎠⎥⎦ ⎝ ⎝ Tω ⎞ ⎛ Tω ⎞ ⎛ = ⎜ sinc ⎟ ⎜ j2sin ⎟ 4π ⎠ ⎝ 4 ⎠ ⎝

(3.1.46)

Odtud, protože w(t ) je integrál w& (t ) , uplatněním věty o integraci v časové oblasti získáváme W (ω ) =

1 ⎛ Tω ⎞⎛ Tω ⎞ ⎜ sinc ⎟⎜ j2sin ⎟ + π W (0 )δ (ω ) jω ⎝ 4π ⎠⎝ 4 ⎠

=

2 ⎡ sin (Tω 4 ) ⎤ Tω ⎢ ⎥ sin ω ⎣ Tω 4 ⎦ 4

=

T sin 2 (Tω 4 ) 2 (Tω 4 )2

(3.1.47)

Nalezli jsme tedy transformační dvojici

F{w(t )} =

T ⎛ Tω ⎞ sinc 2 ⎜ ⎟ 2 ⎝ 4π ⎠

(3.1.48)

Konvoluce

F {u (t ) ∗ w(t )} = U (ω )W (ω )

(3.1.49)

Fyzikální signál se spojitým časem je vždy spojitý. V inženýrské praxi se setkáváme s modely signálů se spojitým časem, které jsou nespojité. Je-li w(t ) v t 0 nespojitý, pak inverzní Fourierova transformace generuje střední hodnotu z w(t 0 + ) a w(t 0 − ) tj w(t 0 + ) + w(t 0 − ) 1 ∞ = W (ω ) e jω t dω ∫ − ∞ 2 2π

(3.1.50)

bez ohledu na to, jaká je hodnota w(t 0 ) . V inženýrské praxi mají dva signály se spojitým časem, jejichž funkční popis se od sebe liší jen v izolovaných bodech, na systém se spojitým 38

časem stejný vliv a mohou být považovány za stejný signál. Můžeme tedy předpokládat, že mezi signálem w(t ) a jeho Fourierovou transformací W (ω ) je jednoznačný vztah. Stručná tabulka transformačních dvojic Fourierovy transformace je uvedena níže. Mnoho dalších potřebných transformačních dvojic čtenář najde pomocí MATLAB Symbolic Toolbox. Dalším zdrojem pro vyhledání transformačních dvojic Fourierovy transformace je INTERNET, zde pod heslem Fourierova transformace pairs.

signál

w(t )

W (ω)

Konstanta

1

δ (ω )

Diracův impuls

δ (t )

1

Diracův impuls v t=t0

δ (t − t0 )

Heavisideův skok

η(t )

Signum t

sgn (t)

Pravoúhlý impuls

Π⎛⎜ t ⎞⎟

T sinc

Exponenciální signál

e − atη (t )

( jω + a )−1

Fázor

e j (ω0t +θ )

2π e jθ δ (ω − ω0 )

Kosinus

cos(ω0 t + θ )

π e jθδ(ω − ω0 ) + e − jθδ(ω + ω0 )

n=∞

2π k =∞ 2π δ (ω − nω0 ) , ω0 = ∑ T k = −∞ T

⎝T ⎠

Časový sled Diracových impulsů

∑ δ (t − nT )

n = −∞

Tabulka 3.1.1

e − jω t 0

π δ (ω ) +

1 jω

2 jω

Tω 2π

[

]

Transformační dvojice Fourierovy transformace

39

3.2

Diskrétní spektrum signálu

Definice 3.2.1

Periodický signál w(t ) ≡ w A (t ) vyhovuje rovnici w(t + kT A ) = w(t )

(3.2.1)

k = L, − 2, − 1, 0, + 1, + 2,L kde je základní perioda signálu w(t )

TA

fA =

1 TA

je základní frekvence signálu w(t )

ω A = 2π f A

je základní úhlový kmitočet signálu w(t )

V těchto skriptech aproximujeme pomocí diskrétního spektra periodický signál o periodě T A . Poznamenejme, že pomocí diskrétního spektra můžeme aproximovat v konečném časovém intervalu (a, a + T A ) i každý jiný fyzikální signál se spojitým časem.

Poznámka

Je-li signál w(t ) periodický s periodou T A , je spojité spektrum W (ω) signálu w(t ) W (ω) =

m =∞

∑ w mδ(ω − mω A )

(3.2.2)

m = −∞

kde 1 TA je základní úhlový kmitočet signálu w(t ) 1 a +TA wm = w(t )e− jωmt dt ∫ TA a jsou koeficienty Fourierovy řady funkce w(t ) ω A = 2π

(3.2.3)

Spojité spektrum periodického signálu obsahuje Diracovy impulsy. Pro praktické výpočty to není pohodlné. Bylo proto definováno a zavedeno tzv. diskrétní spektrum signálu. Informace o periodickém signálu, obsažená ve spektru je plně obsažena i v jeho diskrétním spektru.

40

Definice 3.2.2

Dirichletovy podmínky. Dirichletovy podmínky jsou postačující podmínky pro to, aby y (t ) byla rozvinutelná do konvergentní Fourierovy řady. y (t ) je na daném konečném časovém intervalu, respektive na periodě T A , absolutně integrovatelná, tj. TA

∫0

y (t ) dt ≤ K < ∞

(3.2.4)

kde K je konečná konstanta y (t ) má na daném konečném časovém intervalu, respektive na

periodě T A , konečný počet nespojitostí prvého řádu a konečný počet maxim a minim.

Věta 3.2.1

Každý signál w(t ) , který je na daném konečném časovém intervalu

(a, a + T A ) fyzikální, (tj. signál, který na daném konečném časovém intervalu nese konečnou energii) lze na tomto časovém intervalu (a, a + T A ) vyjádřit pomocí komplexní Fourierovy řady wFS (t ) =

m =∞

∑

m = −∞

w me jωmt

(3.2.5)

kde komplexní koeficienty w m (fázory) jsou dány vztahem wm =

1 a +TA w(t )e− jωmt dt TA ∫a

(3.2.6)

kde ω m = mω A ωA =

2π TA

ω A je základní úhlový kmitočet wFS (t ) T A je perioda wFS (t )

41

Poznámka: Uvedeme souvislost (3.2.5) a (3.2.6) se zápisem Fourierovy řady v reálném oboru. Ze vztahů (3.1.4), (3.2.7)

(3.1.4) wm (t ) = wm, max cos(ωmt + θ m ) = wm, max (cos ωmt cos θ m − sin ωm t sin θ m )

(3.2.7)

= am cos ωmt + bm sin ωm t vidíme, že harmonický signál s fázovým posunem lze složit z neposunutých kosinusového a sinusového signálu. Pro jejich amplitudy snadno odvodíme a m = wm, max cos θ m bm = − wm, max sin θ m

(3.1.5)

(3.2.8)

(3.1.6)

(3.2.9)

(3.1.7)

(3.2.10)

a naopak 2 + b2 wm, max = am m b θ m = − tan −1 m am

Ze základního kurzu matematiky je známo, že každý periodický signál w(t ) s periodou T A pro který platí předchozí podmínky, může být v reálném oboru vyjádřen inversní Fourierovou řadou ve tvaru w(t ) = a0 +

∞

∑ am cos mω At + bm sin mω At

m =1 ∞

(3.2.11)

∑ am cos mω At + bm sin mω At

=

m =0

tj. w(t ) =

∞

∑

m=0

Am cos(mω At + θ m )

(3.2.12)

kde 2 + b2 Am = a m m b θ m = − tan −1 m am

Am ≥ 0

; ;

θ m ∈ (− π;+ π]

(3.2.13) (3.2.14)

Pro reálné koeficienty a m a bm platí vztahy 42

2 a +T A w(t ) cos mω At dt T A ∫a 2 a +T A bm = w(t )sin mω At dt T A ∫a am =

(3.2.15) (3.2.16)

Všimněme si, že a0 =

1 a +T A w(t ) dt T A ∫a

(3.2.17)

Množina amplitud harmonických složek signálu Am = wmax, m tvoří jednostranné diskrétní amplitudové spektrum, respektive jednostranné amplitudové spektrum signálu w(t ) . Množina počátečních fází θ m tvoří jednostranné diskrétní fázové spektrum, respektive jednostranné fázové spektrum signálu w(t ) . Jednostranné diskrétní spektrum signálu ukazuje, že na kmitočtu ω m ≥ 0 signál přispívá ke svému celkovému průběhu přírůstkem Am cos(ω m t + θ m ) = wmax, m cos(ω m t + θ m ) (3.2.18) Pomocí Eulerova vzorce vztahy (3.2.11), (3.2.12), (3.2.15-16) snadno upravíme do dnes již výhradně používaných tvarů (3.2.5) a (3.2.6). Vyjdeme ze vztahu (3.2.11) pro w(t ) , který opíšeme do vztahu (3.2.19) w(t ) =

∞

∑ am cos mω At + bm sin mω At

(3.2.11)

(3.2.19)

m =0

Upravíme e jmω A t + e − jmω A t a m cos mω At + bm sin mω At = a m 2

e jmω A t − e − jmω A t + bm 2j

e jmω A t + e − jmω A t e jmω A t − e − jmω A t = am − jbm 2 2 a − jbm jmω A t a m + jbm − jmω A t e e = m + 2 2 (3.2.20) Označíme

a − jbm wm = m 2 am + jbm w −m = 2

(3.2.21)

Je tedy w m = w*− m

(3.2.22)

43

Zopakujeme, že označujeme mω A = ω m

(3.2.23)

Po tomto označení je e jmω A t = e jω m t e

− jmω A t

=e

(3.2.24)

− jω m t

Je tedy e jω m t = ⎛⎜ e − jω m t ⎞⎟ ⎝ ⎠

*

(3.2.25)

Získáme am cos ωmt + bm sin ωmt = w me jωmt + w − me − jωmt

(3.2.26)

Sumu (3.2.19) nyní již jen upravíme: w(t ) = = =

∞

∑ am cos mω At + bm sin mω At

m =0 ∞

∑

m =0 ∞

∑

w m e j ω m t + w − m e − jω m t

(3.2.27)

w me jωmt

m = −∞

Což je vztah (3.2.6)

w FS (t ) =

m =∞ jω t ∑ wm e m m = −∞

(3.2.6)

(3.2.27)

Využili jsme skutečnosti, že pro w me jωmt platí *

w me jωmt = ⎛⎜ w − me − jωmt ⎞⎟ ⎝ ⎠

(3.2.28)

Konečně odvodíme vztah (3.2.5) pro w m . Podle (3.2.21) je a − jbm wm = m 2

(3.2.21)

(3.2.29)

Dosadíme za a m a bm ze vztahů (3.2.15) a (3.2.16):

44

a − jbm wm = m 2 2 a +TA [w(t )cos mω At − jw(t )sin mω At ] = dt TA ∫a 2 2 a +TA ⎛ cos ωmt − j sin ωmt ⎞ = w(t )⎜ ⎟ dt ∫ 2 TA a ⎝ ⎠ =

(3.2.6)

(3.2.30)

1 a +TA w(t )e − jωmt dt TA ∫a

Což je hledaný vztah (3.2.6). Poznamenejme, že je-li w(t ) ohraničená periodická funkce která má v každé periodě nejvíce konečný počet lokálních maxim a minim a konečný počet bodů nespojitosti prvého řádu (Dirichletovy podmínky), pak Fourierova řada wFS (t ) funkce w(t ) konverguje k w(t ) ve všech bodech kde je w(t ) spojitá, a konverguje ke střední hodnotě limity w(t ) zleva a limity w(t ) zprava v bodě nespojitosti w(t ) . V místě nespojitosti se objevují překmity, dosahující

pro n → ∞ minimální hodnoty 18 % z hodnoty nespojitosti. Tento jev se nazývá Gibbsův jev. Je-li w(t ) fyzikální signál periodický na konečném časovém intervalu, pak Fourierova řada signálu w(t ) konverguje k w(t ) . Je-li w(t ) libovolná periodická funkce, která splňuje Dirichletovy podmínky, lze její integrál najít integrací členů její Fourierovy řady. Je-li w(t ) spojitá periodická funkce, která splňuje Dirichletovy podmínky a splňuje-li Dirichletovy podmínky také její derivace, pak, existuje-li, lze její derivaci najít derivováním členů její Fourierovy řady. Definice 3.2.3

Diskrétní spektrum periodického signálu w(t ) je uspořádaná množina komplexních koeficientů

{wm }

(3.2.31)

45

Definice 3.2.4

Diskrétní spektrum amplitudy periodického signálu w(t ) je uspořádaná množina hodnot

{ wm }

( 3.2.32)

Diskrétní spektrum signálu ukazuje, že na fyzikálním kmitočtu ωm ≥ 0 signál přispívá ke svému celkovému průběhu přírůstkem

wm (t ) = w m e jωmt + w − m e − jωmt wmax,m j (ω t + θ ) wmax,m − j (ω t + θ ) m m e m m + e 2 2 = wmax,m cos(ωmt + θm ) = Am cos(mω A t + θm )

Příklad 3.2.1

Diskrétní Fourierovo spektrum časového sledu pravoúhlých impulsů Uvažujme časový sled pravoúhlých impulsů w(t ) o periodě T0 = 20 , amplitudě 1, šířce 10. Jeho základní úhlový kmitočet je ω0 = 2π

1 π = = 0.31416 . Abychom se T0 10

při demonstraci principu analytického postupu výpočtu vyhli práci s komplexními koeficienty volíme signál jako sudou funkci času. Sled impulsů je nakreslen na Obrázku 3.2.1 A. Vyjádříme w(t ) Fourierovou řadou. Vypočítáme w0 =

1 5 1 w(t) dt = = 0.5 ∫ 20 -5 2

π π ⎡ − jm π 5 jm 5 ⎤ − jm t 5 1 1 1 5 ⎢e 10 dt = 10 − e 10 ⎥ = 1 sin ⎛⎜ m π ⎞⎟ 1e− jmω0t w(t) dt = wm = e ∫ ∫ ⎢ ⎥ mπ ⎝ 2 ⎠ π 20 −5 20 −5 − 20 jm ⎢⎣ ⎥⎦ 10

po dosazení

46

1 π = w −2 = 0 −1 = w −3 = 3π = w −4 = 0

w1 = w −1 = w2 w3 w4

w 5 = w −5 =

1 5π L

Dosadíme do vztahu (3.2.5)

w(t ) =

m = +∞

∑

w mexp( j ω mt )

m = −∞

π ⎞ π π ⎞ π π ⎞ ⎛ π ⎛ ⎛ − j t ⎟ − 1 ⎜ j3 t − j3 t ⎟ 1 ⎜ j5 t − j5 t ⎟ 1 ⎜ j 10 t 10 ⎟ + L 10 ⎟ + w(t ) = 0.5 + ⎜ e + e 10 ⎟ + ⎜ e 10 + e ⎜ e 10 + e π⎜ ⎟ 3π ⎜ ⎟ 5π ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠

Koeficienty jsou vykresleny na Obrázku 3.2.1 D jako funkce indexu m , odpovídajícího úhlovému kmitočtu ω = m ω0 = m

π ω a frekvenci f = 2π 10

Výše uvedenou rovnici lze ihned přepsat na kosinovou Fourierovu řadu, π π ⎞ π π ⎞ π ⎞ ⎛ ⎛ ⎛ π − j5 t ⎟ − j3 t ⎟ 1 ⎜ j5 t − j t ⎟ − 1 ⎜ j3 t 1 ⎜ j 10 t 10 ⎟ + 10 ⎟ + L w(t ) = 0.5 + ⎜ e + e 10 ⎟ + ⎜ e 10 + e ⎜ e 10 + e π⎜ ⎟ ⎟ 5π ⎜ ⎟ 3π ⎜ ⎠ ⎝ ⎠ ⎝ ⎠ ⎝ π −1 π π 1 1 = 0 .5 + 2 cos t + 2 cos 3 t + 2 cos 5 t + L 10 3π 10 5π 10 π

(3.2.33) Analytickou integraci v komplexním oboru obvykle nahradíme, ve složitějších úlohách, numerickými výpočty na počítači. Pro výpočet Fourierových koeficientů w k pak můžeme použít MATLAB a nějakou jednoduchou numerickou metodu. Můžeme ovšem také použít

47

analytický MATLAB Symbolic Toolbox, atd. V těchto skriptech demonstrujme pro nalezení diskrétního spektra signálu použití Diskrétní Fourierovy transformace, konkrétně jejího algoritmu fft tzv. Fast Fourier Transform. Fast Fourier Transform, Rychlá Fourierova Transformace, je algoritmus a software, který počítá koeficienty Wk tzv. Discrétní Fourierovy Transformace, DFT, signálu s diskrétním časem w[n] . Signál w[n] může vzniknout ekvidistantním vzorkováním, odečtem hodnot, signálu w(t ) , potom w[n] = w(nTs ) , kde Ts je perioda vzorkování signálu, n = L ,−2,−1,0,1,2,L Nejdříve vzorkujeme jednu periodu signálu se spojitým časem w(t ) , N −1 ⎤ ⎡ T T ⎤ ⎡ N −1 t ∈ ⎢− , ⎥ ≡ ⎢− Ts , Ts ⎥ v N bodech. Vzorkovací interval Ts volíme konstantní. 2 2 ⎦ ⎣ 2 2⎦ ⎣ Počet bodů N volíme lichý. Vzorkováním signálu w(t ) získáváme signál s diskrétním časem w[n] , w[n] = w(nTs ) . Hodnoty w[n] uložíme do N rozměrného vektoru MATLAB, zde do vektoru s názvem w,následovně: část pro t ≥ 0 , poté, opět vzestupně, část pro t < 0 , tj. w[0], w[1],L, w[( N − 1) / 2], w[− ( N − 1) / 2],L, w[− 1]

(3.2.34)

Nyní je naším postupným cílem najít N koeficientů Wk výše zmíněné Diskrétní Fourierovy Transformace, DFT, signálu s diskrétním časem w[n] Vztah pro DFT{w[n ]}je vztah ( 1. 43 )

Wk =

N −1 2

∑ w[n]e − j 2πkn N

−

, k =−

N −1 2

N −1 N −1 L − 2,−10,1,2,L , 2 2

(3.2.35)

Okamžitě vidíme, že vztah (3.2.35) aproximuje numerický výpočet Fourierova integrálu (3.2.6). Odtud již snadno najdeme hledané koeficienty w k diskrétního Fourierova spektra spojitého původního signálu w(t ) : T 1 w k ≅ s Wk = Wk T N

(3.2.36)

48

V MATLABu nemáme k disposici funkci, která by počítala přímo sumu (3.2.35) od − do +

N −1 2

N −1 , máme zde funkci fft, Fast Fourier transform, která počítá tuto sumu od nuly 2

do N − 1 . Lze však ukázat, že Wk je periodické s periodou N , Wk = Wk + N . Hledaný MATLAB vektor Wk s DFT členy Wk získáme tedy například tak, že nejdříve spočteme členy Wk v jiné části

periody a to DFT koeficienty Wk ≡ Wfft k signálu w[n] pro k = 1,2, L, N ,

Wfft k ≡ Wk =

N −1

∑ w[n]e − j 2πkn N

n=0

, k = 1,2, L, N

(3.2.37)

uložíme výsledky do vektoru Wfftk a poté vektor Wfftk přerovnáme do vektoru Wk. Pro celý výpočet použijeme pouze dvě instrukce. Velmi výkonnou funkci MATLAB fft, kde ovšem fft znamená Fast Fourier Transform a funkci fftshift pro přerovnání, on je to častý požadavek, koeficientů. Algoritmus je velmi jednoduchý: Wfftk=fft(w); Wk=fftshift(Wfftk);

(3.2.38)

kde Wfftk=fft(w) počítá DFT koeficienty Wfft k ≡ W k pro k = 0,1,L, N Wk=fftshift(Wfftk) přerovnává koeficienty Wfft k ≡ W k generované fft(w)do

MATLAB vektoru W s DFT členy Wk pro k = −

N −1 N −1 L − 2,−10,1,2, L, 2 2

Zapišme jednu periodu našeho signálu w[n] do vektoru w v MATLABu následovně: w=[ones(1,5) 0.3 zeros(1,9) 0.3 ones(1,4)];

Zvolili jsme pro numerickou integraci tedy tyto parametry: počet vzorků signálu w(t )

N = 21

49

interval vzorkování signalu w(t )

Ts = 1

Odpovídající program v MATLAB může mít tvar: T=20; wt=[zeros(1,50) ones(1,100) zeros(1,100) ones(1,100) zeros(1,100) ones(1,100) zeros(1,50) ]; t=-30:.1:29.9; %vykreslení signálu w(t) v subplotu (411) subplot(411) plot(t,wt) axis([ -30 30 -1 2]); xlabel('t') ylabel('w(t)') grid N=21; %ulozeni jedne periody signalu w[n] Ts=T/(N-1); w=[ones(1,5) 0.3 zeros(1,9) 0.3 ones(1,5)]; Wfftk=fft(w); subplot(412) stem(abs(Wfftk)) axis([ 1 21 0 20]); xlabel('k') ylabel('abs(Wfftk)') grid

%vypocet FFT signalu w[n] %vykresleni FFT signalu w[n]

Wk=fftshift(Wfftk); subplot(413) k=-10:1:10; stem(k,abs(Wk)); axis([-10 10 0 20]); xlabel('k') ylabel('abs(Wk)') grid

%vypocet DFT signalu w[n] %vykresleni DFT signalu w[n]

wk=Wk*Ts/T; subplot(414)

%vypocet disk. Four. spektra signalu w(t) %vykresleni disk. Four. spektra signalu w(t)

k=-10:1:10; stem(k, abs(wk)) axis([-10 10 0 1]); xlabel('k') ylabel('abs(wk)') grid

50

w(t) abs(Wfftk)

2 1 0 -1 -30

-20

-10

0 t

10

20

30

20 10 0

2

4

6

8

10

12

14

16

18

20

abs(wk)

abs(Wk)

k 20 10 0 -10

-8

-6

-4

-2

0 k

2

4

6

8

10

-8

-6

-4

-2

0 k

2

4

6

8

10

1 0.5 0 -10

Obrázek 3.2.1

Příklad 3.2.1 FFT, DFT a diskrétní Fourierovo spektrum časové posloupnosti pravoúhlých impulsů

Obrázek 3.2.1A ukazuje časový průběh časové posloupnosti pravoúhlých impulsů w(t ) Obrázek 3.2.1B ukazuje FFT amplitudové spektrum Wfft k vzorkovaných pravoúhlých

impulsů w(t ) Obrázek 3.2.1C ukazuje DFT amplitudové spektrum Wk vzorkovaných pravoúhlých impulsů w(t )

Obrázek 3.2.1D ukazuje diskrétní Fourierovo amplitudové spektrum w k pravoúhlých impulsů w(t )

Poznamenejme, že signál wFS (t ) lze generovat z Wk přímo použitím funkce MATLAB ifft, Inverzní Rychlá Fourierova transformace. Ukáže se, že aproximace zadaného sledu

pravoúhlých impulsů w(t ) konstantní složkou a deseti harmonickými je příliš hrubá.

51

Definice 3.2.5

Diskrétní spektrum fáze periodického signálu w(t ) je uspořádaná množina reálných koeficientů

{θ m }

(3.2.39)

kde θm = arg[w m ] Definice 3.2.6

(3.2.40)

Diskrétní spektrum výkonu periodického signálu w(t ) je uspořádaná množina reálných koeficientů

{Pm } = Definice 3.2.7

m

2

(3.2.41)

Konstantní složka w0 (t ) = w0 periodického signálu w(t ) se základní periodou T A je rovna w0 =

Definice 3.2.8

{w }

1 TA 2 w(t ) dt TA ∫−TA 2

(3.2.42)

Harmonické složky w m (t ), m = 1, 2,L periodického signálu w A (t ) jsou wm (t ) = w − m exp(− jωmt ) + w m exp(+ jωmt ) = wm, max cos (ωmt + θm )

(3.2.43)

kde w m = w ∗− m =

1 TA 2 w(t )exp(− jωmt ) dt TA ∫−TA 2

(3.2.44)

jsou koeficienty komplexní Fourierovy řady (3.2.5)

Příklad 3.2.2

Fourierova analýza harmonických složek signálu Uvažujme signál w(t ) generovaný programem MATLAB t=0:0.01:2*pi w1=-cos(t); w2=0.7*cos(10*t); w3=1.2*sin(12*t); w=w1+w2+w3;

52

Signál w(t ) je zobrazen na Obrázku 3.2.2 A . Z Obrázku 3.2.2 A stěží rozeznáme kmitočty harmonických složek, ze kterých se signál skládá. K odhalení kmitočtů harmonických složek použijeme Fourierovu analýzu. Příkaz fft programového prostředí MATLAB W=fft(w);

vypočítává DFT koeficienty

k

Wk

k

0e+002 *

Wk 0e+002 *

1

-0.0022

8

0.0095 - 0.0070i

2

-3.1420 - 0.0114i

9

0.0134 - 0.0094i

3

0.0071 - 0.0014i

10

0.0252 - 0.0132i

4

0.0060 - 0.0021i

11

2.2017 + 0.0528i

5

0.0060 - 0.0030i

12

-0.0213 - 0.0472i

6

0.0066 - 0.0040i

13

0.1430 - 3.7682i

7

0.0076 - 0.0053i

14

0.0093 + 0,0514i

...

...............

...

..............

ukázané na Obrázku 3.2.2 B1, B2 Významné jsou koeficienty W2 , W11 a W13 a jejich komplexně sdružené koeficienty W629 , W620 a W618 odpovídající úhlovým kmitočtům ω = 1, ω = 10 a ω = 12 .

53

Obrázek 3.2.2

Příklad 3.2.2 Analýza harmonických složek signálu

Zde T = 2π , krok jsme zvolili Ts = 0.01 , odtud perioda FFT N = 2π / 0.01 . Odpovídající komplexní koeficienty Fourierova spektra signálu se spojitým časem w(t ) w1 , w10 , w12 , w −1 , w −10 , w −12 , jsou = W2 N = − 0.5001 − 0.0018i w10 = W11 N = 0.3504 + 0.0084i w12 = W13 N = 0.0228 − 0.5997i w1

w −1 = W629 N = − 0.5001 + 0.0018i w −10 = W620 N = 0.3504 − 0.0084i w −12 = W618 N =

0.0228 + 0.5997i

54

Tyto koeficienty popisují hledané harmonické složky w1 (t ) = −1 cos(1t )

w2 (t ) = 0.7 cos(10t ) w3 (t ) = 1.2 sin (12t )

(3.2.45)

signálu w(t ) , ukázané na Obrázku 3.2.2 C, D, E. Čtenář si může postup ověřit na variantách obdobných úloh. Výsledný signál wFS (t ) = w1 (t ) + w2 (t ) + w3 (t ) je zakreslený na Obrázku 3.2.2 F. Úplný program v MATLABu může být například následující

Ts=0.01; T=2*pi; N=T/Ts; t=0:Ts:2*pi; w1=-1*cos(t); w2=0.7*cos(10*t); w3=1.2*sin(12*t); w=w1+w2+w3;

% time step Ts % signal period

subplot(611) plot(t,w) axis([0 2*pi -3 3]) ylabel('w(t)') grid W=fft(w); subplot(623) stem(abs(W)) axis([1 30 0 1000]) ylabel('abs(Wk)') grid

% fft MATLAB function

subplot(624) stem(abs(W)) axis([600 630 0 1000]) ylabel('abs(Wk)') grid w1=W(2)/N.*exp(j*t) + W(629)/N.*exp(-j*t); w2=W(11)/N.*exp(j*10*t) + W(620)/N.*exp(-j*10*t); w3=W(13)/N.*exp(j*12*t) + W(618)/N.*exp(-j*12*t); subplot(613) plot(t,w1) axis([0 2*pi -3 3])

55

ylabel('w1(t)') grid subplot(614) plot(t,w2) axis([0 2*pi -3 3]) ylabel('w2(t)') grid subplot(615) plot(t,w3) axis([0 2*pi -3 3]) ylabel('w3(t)') grid wFS=w1+w2+w3; subplot(616) plot(t,wFS) axis([0 2*pi -3 3]) ylabel('wFS(t)') grid

Příklad 3.2.3

Numerický výpočet diskrétního Fourierova spektra %Urcuje se diskretni Fourierovo spektrum w(k) periodickeho signalu w(t). % %************************************************************************** % Aby se predeslo problemum s citelnosti v oznaceni promennych, oznacime v % tomto prikladu w(k)=c(k) %************************************************************************** % %Signal w(t)je tvoren casovou posloupnosti pravouhlych impulsu f(t), %kde f(t)je sudy jednotkovy pravoúhlý impuls o delce 10s a stride 5s. %Pocitaji se koeficienty c(k)diskretniho Fourierova signalu w(t) % %Signal w(t) budeme studovat pro periodu t=[-5,5]. %Ulozime 1001 time-ekvidistantnich hodnot signalu w(t) do MATLAB %vektoru w (1,1001). % %K disposici mame 1001 vzorku signalu w(t). Muzeme tedy urcit az %1001 koeficientu c(k), pro k=-500,-499,...,-2,-1,0,1,2,...,499,500 % %Urceme si maximální hodnotu omega pro maximální k v c(k): %fmax = fsampling/2.Hz = 100/2.Hz = 50.Hz %omegamax = 2.pi.fmax.rad/s = 2.pi.100/2.rad/s = 100pi.rad/s % %omega0 = 2.pi.f0.rad/s = 2.pi/T.rad/s = 2.pi/10.rad/s %dt = 10/1000.s = 1/100.s % %Ulozime 1001 kmitoctove-equidistantnich hodnot koeficientu spektra %c(k) do MATLAB vektoru c(1,1001). %

56

%Posice clenu vektoru w : %w(t=-5)=w(nn=1), w(t=0)=w(nn=501), w(t=5)=w(nn=1001) %Vztah mezi casem t v w(t) a indexem nn v MATLAB vektoru w(nn): %t=(nn-1)/100-5 % %Posice clenu vektoru c : %c(k=-500)=c(kk=1), c(k=0)=c(kk=501), c(500)=c(kk=1001) %Vztah mezi indexem k v c(k) a indexem kk v MATLAB vektoru c(kk): %k=kk-501 % %MATLAB program: % omega0=2*pi/10; dt=1/100; w=[zeros(1,250) ones(1,501) zeros(1,250)]; %definice signalu w(t) n=0:1:1000; subplot(411); plot(n/100-5,w) %vykresleni signalu w(t) grid; axis([-5 5 -1 2]); xlabel('t'); ylabel('w(t)'); for kk=1:1:1001; c(kk)=0; %nulovani vektoru c end; for kk=501:1:1001; for nn=1:1:1001; c(kk)=c(kk)+w(nn)*exp(-j*(kk-501)*omega0*((nn-1)/1000*10-5))*dt; end; c(kk)=c(kk)/10; end; for kk=1:1:500; c(kk)=conj(c(1002-kk)); end; k=-15:1:15; %rozsak k pro vykresleni c(k) subplot(412) stem(k,abs(c(k+501))) %vykresleni spektra c(k) grid xlabel('k'); ylabel('abs(c(k))') % %Rekonstrukce signalu w(t) z jeho spektra c(k)vypocitaného vyse uvedenym %programem. %wrec(t)je rovno FS-1{c(k)}. %MATLAB program: % for nn=1:1:1001; wrec(nn)=0; %nulovani vektoru wrec end; for nn=1:1:1001; for kk=1:1:1001; wrec(nn)=wrec(nn)+c(kk)*exp(j*kk*omega0*((nn-1)/1000*10-5)); end; end; subplot(413) t=-5:0.01:5; plot(t,abs(wrec)) %vykresleni signalu wrec(t) grid axis([-5 5 -1 2]); xlabel('t');

57

wrec10(t)

wrec(t)

abs(c(k))

w(t)

ylabel('wrec(t)'); %Rekonstrukce signalu w(t) z jeho deseti spektralnich slozek c(k), %vypocitanych vyse uvedenym programem. %MATLAB program: for nn=1:1:1001; wrec10(nn)=0; %nulovani vektoru wrec10 end; for nn=1:1:1001; for kk=491:1:511; wrec10(nn)=wrec10(nn)+c(kk)*exp(j*kk*omega0*((nn-1)/1000*10-5)); end; end; subplot(414) t=-5:0.01:5; plot(t,abs(wrec10)) %vykresleni signalu wrec10(t) grid axis([-5 5 -1 2]); xlabel('t'); ylabel('wrec10(t)');

2 1 0 -1 -5

-4

-3

-2

-1

0 t

1

2

3

4

5

1 0.5 0 -15

-10

-5

0 k

5

10

15

2 1 0 -1 -5

-4

-3

-2

-1

0 t

1

2

3

4

5

2 1 0 -1 -5

-4

-3

-2

-1

0 t

1

2

3

4

5

Obrázek 3.2.3 Numerický výpočet diskrétního Fourierova spektra a rekonstrukce signálu

58

Definice 3.2.9

Náhodný signál je deterministický signál, který lze popsat náhodnou funkcí W (t ) , nazývanou náhodný proces.

Definice 3.2.10

Náhodný proces W (t ) se nazýva stacionární, jsou-li jeho statistické charakteristiky časově invariantní.

Z toho vyplývá, že stacionární náhodný signál, např. šum, který je popsaný stacionární náhodnou funkcí W (t ) , lze popsat jejími statistickými charakteristikami, jako jsou například spektrální hustoty, korelační funkce, atd.

Signály s konečným spektrem a signály s konečným časem

Signál w(t ) se nazývá signál s konečným spektrem, je-li jeho Fourierova transformace W (ω) rovná nule pro ω > B , kde B je konečné kladné číslo, nazývané šířka základního pásma signálu. Je-li signál w(t ) ohraničený šířkou pásma B, nemá žádné harmonické složky nebo harmonické elementy o kmitočtech vyšších než B. Z toho vyplývá, že každý signál s konečným spektrem musí mít nekonečné trvání, neboli, signály s konečným spektrem nemohou být tzv. signály s konečným časem. (Signál w(t ) se nazývá signál s konečným časem, jestliže existuje kladné číslo T takové, že w(t ) = 0 pro všechna t < −T a t > T .) Z toho vyplývá, že každý signál s konečným časem musí mít nekonečné spektrum, nekonečnou šířku základního pásma. Navíc, protože všechny fyzikální signály jsou signály s konečným časem, musí mít všechny nekonečnou šířku základního pásma. Nicméně, pro každý v praxi se vyskytující rozumný deterministický signál s konečným časem w(t ) lze dokázat, že jeho Fourierova transformace W (ω) konverguje k nule když se ω → ∞ .

59

Lze proto vždy předpokládat že W (ω) ≈ 0 pro všechna ω > B kde se B zvolí dostatečně velké.

Definice 3.2.11

Šířka pásma B signálu w(t ) B = ω max − ω min

(3.2.46)

kde ω min je minimální úhlový kmitočet, na kterém má signál w(t )

nenulovou hodnotu spektra energie nebo spektra výkonu ω max je maximální úhlový kmitočet, na kterém má signál w(t ) nenulovou hodnotu spektra energie nebo spektra výkonu

Definice 3.2.12

Šířka pásma signálu v základním pásmu w(t ) B = ω max

(3.2.47)

kde ω max je maximální úhlový kmitočet, na kterém má signál w(t ) nenulovou hodnotu spektra energie nebo spektra výkonu

Definice 3.2.13

Signál v komunikačním pásmu w(t ) Signál v komunikačním pásmu w(t ) má nenulové hodnoty spektra energie nebo spektra výkonu pouze v pásmech, soustředěných kolem úhlových kmitočtů ± ω c , často platí ω c >> B

(3.2.48)

60

4

ANALÝZA SOUSTAV SE SPOJITÝM ČASEM

V této kapitole se budeme zabývat analýzou soustav se spojitým časem. Soustava se spojitým časem je soustava, jejíž vstupní signály a výstupní signály jsou signály se spojitým časem. Zde problematiku ještě zjednodušíme předpokladem, že i stavové veličiny soustavy, viz dále, jsou signály se spojitým časem.

4.1

Lineární časově invariantní soustava se soustředěnými parametry (LTIL soustava)

Stav soustavy

Cílem je nalézt výstupní signál y (t ) soustavy F , jejíž stav je ovlivněn vstupním signálem u (t ) . Popíšeme stav soustavy F množinou {xi (t )} vhodně zvolených veličin, proměnných

xi (t ) jejichž hodnoty lze měřit na soustavě. Proměnné xi (t ) se nazývají stavové veličiny soustavy, někdy též stavové proměnné soustavy. Stav soustavy F je v libovolném čase t i popsán hodnotami x1 (t i ), x 2 (t i ),L svých stavových veličin. Použijeme minimální množinu n stavových veličin a definujeme stav soustavy stavovým vektorem

x(t ) = [x1 (t ), x2 (t ),L , xn (t )]

(4.1.1)

soustavy F .

Definice 4.1.1

Řád n soustavy F Řád n soustavy F je roven rozměru n stavového vektoru x(t ) soustavy F .

Poznamenejme, že naše schopnost definovat stavové veličiny určité konkrétní soustavy často závisí na naší apriorní znalosti této soustavy. V dalším textu často použijeme pro stručnost následující zkrácenou terminologii

61

x(t ) − xi (t ) −

u (t ) − y (t ) −

stavový vektor soustavy F stavová veličina soustavy F

− −

stav stavová veličina

vstupní signál soustavy F výstupní signál soustavy F

− −

vstup výstup

V tomto textu omezujeme diskusi na LTIL SISO soustavy. Studujeme lineární soustavu H řádu n , u které lze n stavových veličin standardně vybrat tak, že za stavovou veličinu zvolíme její výstup y (t ) a n − 1 prvních derivací výstupu y (t ) . x1 (t ) = y (t ) , x 2 (t ) = y& (t ) , L , x n (t ) = y (n −1) (t )

(4.1.2 )

Tento algoritmus bývá charakteristický pro soustavy nazývané často pojmem black-box system, černá skříňka. Může však vést k požadavku na derivování vstupního signálu. Abychom se tomuto nežádoucímu derivování vyhnuli, uvedeme později standardní algoritmus dekompozice soustavy H na stavovou soustavu S a výstupní soustavu O .

Definice 4.1.2

Černá skříňka (Black-box soustava) Černá skříňka (Black-box soustava) se nazývá soustava, u které není známa její vnitřní struktura.

Černou skříňku popíšeme vztahy (4.1.2) nebo vztahy z (4.1.2) odvozenými.

Obrázek 4.1.1

Klidový, nulový, stav soustavy

Jako klidový, nulový, stav x(t i ) = 0 označíme stav soustavy, nacházející se v naprosto uvolněném stavu. V našem modelu soustava, nacházející se v klidovém, nulovém, stavu, nemění svůj stav, je-li vstup u (t ) = 0 t ≥ t i .

62

Odvodíme matematický model jednoduché soustavy F , která má pouze jeden vstup u (t ) a pouze jeden výstup y (t ) . Tento typ soustavy se nazývá single-input single-output neboli SISO soustava. Obecně budeme soustavu F studovat pro čas t ∈ (− ∞, ∞ ) . Budeme předpokládat znalost těch parametrů, které určují fyzikální vlastnosti soustavy. Například u RLC obvodu budeme předpokládat znalost hodnot prvků R, L a C . Navíc budeme předpokládat, že se charakteristiky a soustavy nemění s časem. Ani struktura RLC obvodu, ani hodnoty jejich prvků nejsou závislé na čase. Uvidíme, že takovéto soustavy se nazývají časově invariantní soustavy. Také budeme předpokládat, že se charakteristiky soustavy nemění s hodnotami stavových veličin. Ani struktura RLC obvodu, ani hodnoty jejich prvků nejsou závislé na napětí na prvcích nebo na proudu který jimi protéká. Uvidíme, že takovéto soustavy se nazývají lineární soustavy. Navíc budeme předpokládat, že soustava má pouze konečný počet stavových veličin. Uvidíme, že takovéto soustavy se nazývají soustavy se soustředěnými parametry (Lumped Systems) Soustava se soustředěnými parametry Definice 4.1.3

Soustava se soustředěnými parametry Soustava F se nazývá soustava se soustředěnými parametry jestliže její stav x(t ) má konečný počet stavových veličin x(t k )

⎫ ⎬ u (t ) t ≥ t k ⎭

y (t ) t ≥ t k

(4.1.3)

kde rozměr n vektoru x(t ) je konečné číslo. Poznamenejme, že soustavu se soustředěnými parametry F lze popsat konečným počtem parametrů. Soustava F se nazýva soustava s rozloženými, distribuovanými, parametry, jestliže její stav má, na rozdíl od modelu soustavy se soustředěnými parametry, nekonečně mnoho stavových veličin. Poznamenejme, že soustavu F s rozloženými, distribuovanými, parametry lze popsat, na rozdíl od modelu soustavy se soustředěnými parametry, nekonečným počtem 63

parametrů. V tomto textu naši diskusi omezíme na soustavy se soustředěnými parametry. Soustava se soustředěnými parametry F může být v libovolném čase t popsána hodnotami minimálního konečného počtu n stavový veličin x(t ) = x1 (t ), x 2 (t ),L , x n (t ) . Stav x(t i ) shrnuje minulé chování systému. V našem modelu shrnuje působení minulého vstupního signálu u (t ) t < t i na systém nacházející se v určitém počátečním stavu x(t ) t = −∞ , na současný výstup y (ti ) a tedy i na budoucí výstup y (t ) t > ti . V mnoha případech budeme soustavu studovat pro t ∈ [t 0 , ∞ ) . Potom budeme předpokládat, že stav x(t 0 ) soustavy v čase t = t 0 je znám a že ho lze určit, vypočítat, pro libovolný čas t > t 0 . Stav x(t 0 ) se nazývá počáteční stav soustavy. Budeme rovněž předpokládat, že vstup u (t ) je znám pro libovolný čas t ≥ t 0 a že výstup y (t ) lze určit, vypočítat, pro libovolný čas t ≥ t 0 . Definice 4.1.4

Fundamentální úloha analýzy soustavy Je dána soustava F , F ( x,u ,t ) = 0 . Je dán počáteční stav soustavy x(t 0 ) . Je dán vstupní signál u (t ) t ≥ t0 . Fundamentální úlohou analýzy soustavy je určit stav x(t ) t ≥ t 0 . x(t 0 ) u (t ) t ≥ t 0

Definice 4.1.5

⎫⎪ ⎬ x(t ) t ≥ t 0 ⎪⎭

(4.1.4)

Základní úloha analýzy soustavy Je dána soustava F , F ( x,u ,t ) = 0 . Je dán počáteční stav soustavy x(t 0 ) . Je dán vstupní signál u (t ) t ≥ t0 . Základní úlohou analýzy soustavy je určit výstup y (t ) t ≥ t0 .

64

x(t 0 ) u (t ) t ≥ t 0

⎫⎪ ⎬ ⎪⎭

y (t ) t ≥ t 0

(4.1.5)

Vzhledem k tomu, že jsme omezili téma naší diskuse LTIL SISO soustavami, liší se fundamentální úloha analýzy soustavy a základní úloha analýzy soustavy od sebe jen metodickými detaily. Pro další text si zvolíme za vodítko spíše základní úlohu analýzy soustavy. Víme, že počáteční stav soustavy x(t 0 ) nese informaci, která, společně se vstupním signálem u (t ) t ≥ t 0 jednoznačně definuje výstupní signál soustavy y (t ) t ≥ t 0 . Počáteční stav

soustavy x(t 0 ) popisuje systém v okamžiku, ve kterém zahajujeme jeho analýzu. x(t 0 ) je výsledkem působení minulého vstupu u (t ) t < t 0 , ale také výsledkem výchozího stavu x(t ) t = −∞ . Bez ztráty obecnosti budeme předpokládat x(t ) = 0 t = −∞ .

Obrázek 4.1.2

Základní úloha analýzy soustavy

Snad v každé fyzikální soustavě lze najít stavové veličiny, které nejsou závislé na vstupním signálu. V teorii soustav se takovéto soustavy nazývají soustavami které nejsou říditelné. Snad v každé fyzikální soustavě lze také najít stavové veličiny, na kterých nezávisí výstupní signál soustavy. V teorii soustav se takovéto soustavy nazývají soustavami které nejsou pozorovatelné. V naší diskusi se těmito modely soustav nebudeme zabývat. Stavové veličiny stavového vektoru tvoří minimální množinu stavových veličin která jednoznačně popisuje stav soustavy. Na druhé straně, výběr stavových veličin pro tento jednoznačný popis jednoznačný není. Například, u RLC obvodu, můžeme volit za stavové veličiny různé kombinace napětí a proudů v jednotlivých větvích obvodu.

65

Výstup může mít jiný fyzikální rozměr než vstup. Například poloha y (t ) tělesa může být regulována mechanickou silou u (t ) , úhlová rychlost, otáčky, y (t ) rotoru mohou být řízeny elektrickým proudem u (t ) , atp. V dalším textu, s výjimkou několika fyzikálních příkladů, bude všem fyzikálním veličinám přiřazen fyzikální rozměr 1 . V této studii se budeme zabývat základními obecnými metodami analýzy soustav. Tyto metody jsou nezávislé na fyzikální podstatě soustav. Přímá aplikace obecných postupů pro řešení skutečných technických problémů může vést k úspěchu jen u velmi jednoduchých úloh. V praktických a složitějších aplikačních úlohách jsou tyto obecné postupy užitečné pouze tehdy, jsou-li dostatečně doplněny znalostmi a detailními informacemi z fyzikálního a technického aplikačního prostředí, ve kterém se analýza a syntéza soustav provádí. Dynamická soustava

Definice 4.1.6

Dynamická soustava Soustava F se nazývá dynamická, závisí-li její stav x(t ) t > t k na stavu x(t k ) . x(t k )

⎫ ⎬ x(t ) t ≥ t k u (t ) t ≥ t k ⎭

(4.1.6)

Poznamenejme, že říditelná pozorovatelná nedynamická soustava má n = 0 stavových veličin. Definice 4.1.7

Ryze dynamická soustava Soustava F se nazývá ryze dynamická, je-li dynamická a je-li y (t k ) nezávislé na u (t k ) , kde u (t k ) je fyzikálně realizovatelný signál.

66

Časově invariantní soustava, time-invariantní soustava Definice 4.1.8

Time-invariantní soustava Soustava F se nazývá časově invariantní, time-invariantní (časově neměnná) jestliže pro každé x (t k )

⎫ ⎬ u (t ) t ≥ t k ⎭

y (t ) t ≥ t k

(4.1.7)

platí x(t k + τ ) = x(t k ) ⎫ ⎬ u (t − τ ) t ≥ t k + τ⎭

y (t − τ ) t ≥ t k + τ

(4.1.8)

V dalším textu omezíme naši diskusi na studium time-invariantních soustav. To nám umožňuje umístit čas t 0 na standardní hodnotu t 0 = 0 a studovat chování soustavy v čase t ∈ [0, ∞ ) .

x(0) = x0 u (t ) t ≥ 0

⎫ ⎬ ⎭

y (t ) t ≥ 0

(4.1.9)

Počáteční stav x(0 ) popisuje stav soustavy v čase t = 0 kdy zahajujeme analýzu soustavy. Základní úloha analýzy soustavy požaduje znalost vstupu u (t ) pouze pro t ≥ 0 . Popis soustavy vztahem (4.1.9) je tedy zcela dostatečný, vyhovující, nevzniká žádný problém, pokud je vstup u (t ) fyzikální. Fyzikální vstup u (t ) t ≥ 0 svým působením na soustavu nemění v čase t = 0 hodnotu počátečního stavu x(0) . V naší studii však budeme často na soustavu přivádět vstup, konstruovaný z významného teoretického modelu signálu známého jako Diracův impuls. Paul A. M. Dirac (1902-1984) definoval tzv. jednotkový impuls, neboli delta impulse δ(t ) , neboli Diracův impuls δ(t ) následovně:

67

Definice 4.1.9

Diracův impuls δ(t ) Diracův impuls δ(t ) je funkcionál splňující vztahy ∞

∫ δ(t )dt = 1

(4.1.10)

−∞

a δ(t ) = 0 for t ≠ 0 .

(4.1.11)

Pro tento text bude stačit představit si Diracův impuls δ(t ) jako velmi úzký sudý pravoúhlý impuls s jednotkovou plochou, s limitní šířkou ε → 0 , výškou

1 → ∞ , ležící v intervalu ε

(0 − ,0 + ) , viz Obrázek 4.1.4. Fourierova transformace Diracova impulsu je F{δ(t )} = 1 . Diracův impuls není fyzikální signál. V časovém intervalu (0 − , 0 + ) nese nekonečný výkon. Takovýto vstupní signál u (t ) může v časovém intervalu (0 − , 0 + ) změnit počáteční stav soustavy x(0 ) . Aby se zabránilo nejasnostem, která část vektoru x(0 ) je minulostí soustavy a která část x(0) je změněna vlivem vstupu u (t ) t = 0 , budeme předpokládat znalost počátečního stavu soustavy v čase t = 0 − . Odpovídající, již konečný, zápis úlohy analýzy soustavy je tedy následující,

x(0 − ) = x0 ⎫ ⎬ u (t ) t ≥ 0 ⎭

y (t ) t ≥ 0

(4.1.12)

viz Obrázek 4.1.3.

Obrázek 4.1.3

Konečný zápis úlohy analýzy soustavy

68

Časové intervaly, ve kterých se proměnné x(t ), y (t ) a u (t ) definují a studují jsou uvedeny na schématu na Obrázku 4.1.3 a nebudou vždy v dalším textu opakovaně uváděny. Poznamenejme, že časový interval [0 − ,0 + ] bude přesně diskutován pouze v jednoduchém příkladu Příklad 4.1.1 a v komentáři k rovnici (4.1.13). Příklad 4.1.1

Integrace Diracova impulsu Soustava sestává z ideálního integrátoru. Stav soustavy je popsán jejím výstupem, x(t ) = y (t ) . Integrátor se v čase t = 0 − nachází v nulovém, klidovém, počátečním stavu x(0 − ) = 0 . Vstup integrátoru je Diracův impuls, u (t ) = δ(t ) .

x(0 − ) = 0⎫ ⎬ δ(t ) ⎭

y (t ) t ≥ 0 =

t

∫ δ(t ) dt + 0 = η(t ) t ≥ 0

(4.1.13)

0−

Výstup integrátoru y (t ) je zobrazen na Obrázku 4.1.4.

Obrázek 4.1.4

Příklad 4.1.1 Výstup integrátoru y (t )

Time-invariance nabízí dvě hlavní výhody. První je, že nám umožňuje u určitých typů úloh najít jejich analytické řešení. Time-invariantní úlohy, úlohy popsané rovnicemi s konstantními koeficienty např., se řeší mnohem snadněji než úlohy, které nejsou time invariantní. Druhá výhoda je mnohem významnější. Uvažujme určitou soustavu. Time-invariantnost nám

69

umožňuje aplikovat výsledky, nebo zkušenosti, získané při řešení určité úlohy na soustavě, při řešení jiné podobné úlohy. Pokud bychom u naší soustavy nepředpokládali určitý stupeň time invariantnosti, nebyl by tento postup možný.

Lineární soustava Definice 4.1.10

Lineární soustava splňuje zákon superposice: Soustava F se nazývá lineární, nebo dále lineární soustava H , platí-li pro každé

x1 (t k )

⎫ ⎬ u1 (t ) t ≥ t k ⎭

y1 (t ) t ≥ t k

(4.1.14)

y 2 (t ) t ≥ t k

(4.1.15)

a pro každé

x 2 (t k )

⎫ ⎬ u 2 (t ) t ≥ t k ⎭ vztahy

a x1 (t k )

⎫ ⎬ ay1 (t ) t ≥ t k au1 (t ) t ≥ t k ⎭ (zákon homogenity)

(4.1.16)

a

x1 (t k ) + x 2 (t k )

⎫ ⎬ u1 (t ) + u 2 (t ) t ≥ t k ⎭ (zákon aditivity)

y1 (t ) + y 2 (t ) t ≥ t k

(4.1.17)

kde a je konečná reálná konstanta. V dalším textu omezíme naši diskusi na studium lineárních soustav. Definice 4.1.11

Odezva při nulovém vstupním signálu y a (t ) Odezva při nulovém vstupním signálu y a (t ) je odezva soustavy na počáteční stav x(0 − ) při nulovém vstupu u (t ) = 0 t ≥ 0 .

70

Definice 4.1.12

Odezva při nulovém počátečním stavu y b (t ) Odezva při nulovém počátečním stavu y b (t ) je odezva soustavy na vstup u (t ) t ≥ 0 při nulovém počátečním stavu x(0 − ) = 0 .

Linearita ospravedlňuje používání většiny algoritmů spojených s řešením základní úlohy analýzy soustav metodou dekomposice problému na sekvenci jednodušších kroků. Například, linearita umožňuje nalézt celkovou odezvu y (t ) soustavy H jako součet odezvy při nulovém vstupním signálu y a (t ) a odezvy při nulovém počátečním stavu yb (t ) , y (t ) = y a (t ) + y b (t ), umožňuje nalézt celkovou odezvu y (t ) soustavy H jako součet přechodné odezvy y aa (t ) a ustálené odezvy y bb (t ) , y (t ) = y aa (t ) + y bb (t ), atd. Odezvy lze najít per partes, odezvu y a (t ) lze vypočítat jako součet dílčích odezev y a1 (t ), y a 2 (t ), L, y ak (t ) na správně zvolené dílčí počáteční stavy x1 (0 − ), x 2 (0 − ),L , x k (0 − ) , ya (t ) = ya1(t ) + ya 2 (t ) + L + yak (t ) , kde x(0 − ) = x1 (0 − ) + x 2 (0 − ) + L + x k (0 − ) , odezvu y b (t ) lze vypočítat jako součet dílčích odezev yb1 (t ), yb2 (t ), L , ybl (t ) to na správně zvolené dílčí vstupy u1 (t ), u 2 (t ), L.ul (t ) ,

yb (t ) = yb1(t ) + yb 2 (t ) + L + ybl (t ) , kde u (t ) = u1(t ) + u2 (t ) + L + ul (t ) , atd. Symbol H je běžně vyhrazen pro označení lineární soustavy. Shrneme-li omezení která jsme si zavedli výše, omezujeme naši diskusi na lineární time-invariantní soustavy se soustředěnými parametry, neboli LTIL soustavy. V následujícím textu proto použijeme symbol H pro označení LTIL soustavy. Zdůrazněme, že jak model soustavy se soustředěnými parametry, tak model time-invariantní soustavy, tak model lineární soustavy jsou v praxi nejčastěji používaná vzájemně nezávislá zjednodušení fyzikální reality. Použitelnost každého modelu závisí na studované soustavě a na řešené úloze. Výše uvedené modely lze aplikovat v libovolných kombinacích a jejich použití není omezeno jinými modely. Například soustavy popisující vedení tepla jsou často popsány lineárními time-invariantními modely s distribuovanými parametry; elektronické spínače jsou studovány s využitím nelineárních time invariantních modelů se soustředěnými parametry, atp. LTIL model nám proto umožňuje zkonstruovat celistvou, jednotnou a inteligentní metodu analýzy soustav na úkor její přítulnosti k fyzikální realitě. 71

Kauzální soustava

Poslední co do pořadí, ale ne poslední co do významu pojem, který zde uvedeme je pojem kauzální soustava. Výstup y (t i ) kauzální soustavy nesmí záviset ani na vstupu u (t ) t > t i , ani na stavu x(t ) t > t i .

Příklad 4.1.2

Ideální dolnopropustný filtr Ideální dolnopropustný filtr H má pravoúhlou kmitočtovou charakteristiku ⎛ ω H(ω) = Π ⎜⎜ ⎝ 2ω max

⎞ ⎟⎟ ⎠

(4.1.18)

kde ω max je kritická frekvence filtru. Počáteční stav filtru je x(0 − ) = 0 . Vstup u (t ) = η(t ) , kde η(t ) je Heavisideův skok (viz dále) η(t ) = 0 t < 0, η(t ) = 1 2 t = 0, η(t ) = 1 t > 0

(4.1.19)

Lze ukázat, viz dále, že výstup je ω y (t ) = max π

t

⎛ω ⎞ sin c ⎜ max t ⎟ dt ⎝ π ⎠ −∞

∫

(4.1.20)

Výstup má nenulovou hodnotu v každém nenulovém časovém intervalu v intervalu

(− ∞, ∞ ) . Filtr tak generuje výstupní signál dříve, než byl buzen nenulovým vstupním signálem. Tento fakt je v rozporu s fyzikálním zákonem kauzality. Závěr je ihned zřejmý a je velmi jednoduchý: ideální dolnopropustný filtr nelze realizovat

72

Příklad 4.1.3

Non-LTIL soustava Výstup zesilovače A je vázán se vstupem rovnicí y (t ) = a u (t ) + b

(4.1.21)

kde a a b jsou nenulové konečné reálné konstanty. Rovnice (4.1.21) je rovnice přímky.

Obrázek 4.1.5

Příklad 4.1.3. Závislost mezi vstupem a výstupem zesilovače

Úloha: je zesilovač A LTIL soustava? Řešení: Soustava se soustředěnými parametry. Zesilovač A je nedynamická soustava. Má dva parametry a , b a n = 0 stavových veličin. Zesilovač A je soustava se soustředěnými parametry. Time-invariantní soustava. Ani struktura, ani parametry zesilovače A se v čase nemění. Soustava A je time-invariantní soustava. Lineární soustava. Pro vstup u1 (t ) = 1 je výstup y1 (t ) = a + b . Pro vstup u 2 (t ) = 2 u1 (t ) = 2 je výstup y 2 (t ) = 2a + b . Aby byla soustava lineární, musel by

být výstup y 2 (t ) roven y 2 (t ) = 2 y1 (t ) = 2a + 2b . Vidíme, že zesilovač A je soustava nelineární, což nemusí být, vzhledem k tvaru rovnice ( 2. 19 ) a přímky na Obrázku 2.5, na první pohled zřejmé.

73

5

ANALÝZA LTIL SOUSTAVY V ČASOVÉ OBLASTI

5.1

Časové odezvy LTIL soustavy

Uvažujme LTIL soustavu, která se nachází v nulovém počátečním stavu x(0 − ) = 0 . Na soustavu přivedeme vstup δ (t ) kde δ (t ) je Diracův impuls a měříme výstup y (t ) . Tento zvláštní výstup označíme h(t ) . Je to odezva LTIL soustavy na jednotkový impuls neboli impulsová charakteristika LTIL soustavy. Odezva soustavy na jednotkový impuls - impulsová charakteristika h (t )

Definice 5.1.1

Impulsová charakteristika h (t ) Impulsová charakteristika h(t ) je výstup LTIL soustavy generovaný vstupem δ (t ) při nulovém počátečním stavu soustavy x(0 − ) = 0 ⎫ ⎬ → h(t ) t ≥ 0 u (t ) = δ(t ) ⎭

(5.1.1)

Stručně:

δ (t ) → h(t ) Je-li soustava kauzální, potom h (t ) = 0 t < 0 . A naopak, tato podmínka je nutná a postačující podmínka pro to, aby byla LTIL soustava kauzální. Poznamenejme, že impulsová charakteristika ideálního dolnopropustného filtru H z Příkladu 4.1.2 je

h(t ) =

ω max ⎛ω ⎞ sinc ⎜ max t ⎟ , má nenulovou hodnotu v každém nenulovém časovém intervalu π ⎝ π ⎠

v intervalu (− ∞, ∞ ) .

Asymptoticky stabilní soustava Definice 5.1.2

Asymptoticky stabilní LTIL soustava ve smyslu Ljapunova LTIL soustava se nazývá asymptoticky stabilní ve smyslu Ljapunova, platí-li pro každé konečné x(0 − )

74

x(0 − ) = x0 u (t ) = 0 t ≥ 0

⎫ ⎬ ⎭

y (t ) → 0 t → ∞

(5.1.2)

kde ∞

∫ y(t ) dt < K , K < ∞

(5.1.3)

0−

Věta 5.1.1

Z Definice 5.1.2 vyplývá, že pro asymptoticky stabilní soustavu H platí vztahy h (t ) → 0 t → ∞

(5.1.4)

∞

∫ h(t ) dt < K , K < ∞

(5.1.5)

0−

V Odstavci 5.2 uvedeme popis LTIL soustavy obyčejnou lineární diferenciální rovnicí s konstantními koeficienty (5.1.6) y (n ) (t ) + an −1 y (n −1) (t ) + L + a0 y (t ) = bm u (m ) (t ) + bm −1u (m −1) (t ) + L + b0u (t )

, m≤n

Uvidíme, že impulsová charakteristika h(t ) asymptoticky stabilní ryze dynamické LTIL soustavy má tvar h(t ) = k1t i −1e Q1t + k 2 t i −1e Q2 t + L + k n t i −1e Qn t

(5.1.7)

kde Q1 , Q2 ,L, Qn jsou kořeny charakteristické rovnice k rovnici (5.1.6) . Za i se postupně dosadí i = 1,2, L , r kde r je násobnost kořenu. Konstanty k1 , k 2 ,L + k n se určí z nulového počátečního stavu soustavy.

Lineární členy LTIL soustavy

LTIL soustavu lze sestavit propojením lineárních členů. Může být překvapující, že přitom vystačíme jen se čtyřmi typy lineárních členů. Jsou to derivátor, integrátor, násobička

75

konstantou a sumátor. Navíc, fyzikální LTIL soustava derivátor obsahovat nemůže. Problematiku ilustruje příklad 5.2.4 v následujícím textu. Odezva při nulovém vstupním signálu y a (t )

Odezva při nulovém vstupním signálu y a (t ) byla definována Definicí 4.1.11 Leonhard Euler (1707 – 1783) ukázal, že odezva LTIL soustavy při nulovém vstupním signálu má tvar y a (t ) = K1t i −1e Q1t + K 2 t i −1e Q2 t + L + K n t i −1e Qn t

(5.1.8)

kde Q1 , Q2 ,L, Qn jsou kořeny charakteristické rovnice k rovnici (5.1.6) a odpovídají exponentům v (5.1.7). Za i se postupně dosadí i = 1,2, L , r kde r je násobnost kořenu. Konstanty K1 , K 2 ,L, K n se určí z počátečního stavu soustavy. Rovnice (5.1.8) vede ke Větě 5.1.2.

Věta 5.1.2

LTIL soustava popsaná rovnicí (5.1.6) je stabilní ve smyslu Ljapunova tehdy a jen tehdy, mají-li kořeny charakteristické rovnice k rovnici (5.1.6) zápornou nebo nulovou reálnou část a kořeny, které mají nulovou reálnou část jsou jednonásobné.

Věta 5.1.3

LTIL soustava popsaná rovnicí (5.1.6) je asymptoticky stabilní ve smyslu Ljapunova tehdy a jen tehdy, mají-li kořeny charakteristické rovnice k rovnici (5.1.6) zápornou reálnou část.

Poznámka:

Vztah (5.1.8) lze upravit použitím identity e jωt = cos ωt + j sin ωt ,

(5.1.9)

y a (t ) = ∑ Bk e −α k t cos(ω k t + ϕ k )

(5.1.10)

do tvaru

k

kde koeficienty Bk , α k , ω k a ϕ k jsou reálné. Vztah (5.1.9) je známý Eulerův vztah. 76

Odezva při nulovém počátečním stavu y b (t )

Odezva při nulovém počátečním stavu y b (t ) byla definována Definicí 4.1.12 Použijeme impulsovou charakteristiku h(t ) a odvodíme matematické rovnice, popisující vztah mezi výstupem a vstupem LTIL soustavy. Aproximujme vstupní signál u (t ) časovou posloupností Diracových impulsů tak, jak je uvedeno na Obrázku 5.1.6. Vzorkovací vlastnost Diracova impulsu říká, že ∞

∫ f (t )δ(t − τ)dt = f (τ)

(5.1.11)

−∞

Získáme u (t ) =

∞

∫ u(τ)δ(τ − t )dτ =

−∞

∞

∫ u(τ)δ(t − τ)dτ

(5.1.12)

−∞

Fyzikální interpretaci vztahu (5.1.12) ukazuje Obrázek 5.1.1, kde pro Δt → dt → 0 máme u (t ) ≈

∞

∞

∑ u(k Δt )δ(t − k Δt )Δt ≈ ∫ u(τ)δ(t − τ)dτ

k = −∞

Obrázek 5.1.1

(5.1.13)

−∞

Aproximace vstupu, vztah (5.1.12)

77

Použijeme notace z Definice 5.1.1, píšeme

δ (t ) → h(t )

(5.1.14)

Využijeme time invariance LTIL soustavy, máme δ(t − τ ) → h(t − τ )

(5.1.15)

Využijeme zákona homogenity, máme u (τ )δ(t − τ)dτ → u (τ )h(t − τ)dτ

(5.1.16)

Využijeme zákona aditivity, máme u (t ) → yb (t ) =

∞

∞

∫ u(τ)h(t − τ)dτ = ∫ u(t − τ)h(τ)dτ

−∞

(5.1.17)

−∞

Výraz na pravé straně vztahu (5.1.17) se nazývá konvoluční integrál nebo stručně konvoluce. U kauzálních soustav h (t ) = 0 t < 0 , je tedy h (t − τ ) = 0 t < τ a interval integrace lze omezit. Získáme u (t ) → yb (t ) =

t

t

∫ u(τ)h(t − τ)dτ = ∫ u(t − τ)h(τ)dτ .

0−

(5.1.18)

0−

Rovnice (5.1.17) a (5.1.18) jsou základní, popisují y b (t ) , odezvu LTIL soustavy při nulovém počátečním stavu. Zápis konvoluce se zkracuje: yb (t ) = u (t ) ∗ h(t ) ≡ u ∗ h (t ) yb (t ) = h(t ) ∗ u (t ) ≡ h ∗ u (t )

(5.1.19) (5.1.20)

Pro t ≥ 0 určíme podle principu superposice celkovou odezvu soustavy jako součet odezvy soustavy při nulovém vstupním signálu y a (t ) a odezvy soustavy při nulovém počátečním stavu yb (t ) : y (t ) = y a (t ) + yb (t )

(5.1.21)

Jedním ze základních závěrů, plynoucích z předchozích úvah je, že jak odezva při nulovém vstupním signálu y a (t ) , tak odezva při nulovém počátečním stavu yb (t ) , mají stejné vlastní dynamické charakteristiky, určené impulsovou charakteristikou h(t ) . V našem modelu soustavy je stav LTIL soustavy ovlivněn minulým vstupem. Výše uvedený závěr bylo tedy možno očekávat.

78

Konvoluce také ihned odpovídá na principielní otázku, jaká je ustálená odezva soustavy asymptoticky stabilní LTIL soustavy na harmonický vstup u (t ) = umax cos(ωu t + θu ) . Harmonický vstupní signál je definovaný na časovém intervalu (− ∞ , ∞ ) . Byl přiveden na soustavu nacházející se v konečném stavu x(t ) t = −∞ . Předpokládáme x(t ) = 0 t = −∞ , to zde ale není důležité. Rovnice ( 2. 30 ) ukazuje, že všechny přechodné děje na soustavě pro t → ∞ konvergují k nule. Rovnice ( 2. 30 ) rovněž ukazuje, že je výstup y (t ) pro

harmonický vstup u (t ) ohraničený. Pro t > t0 je y (t ) =

∞

∫ h(τ)u(t − τ)dτ . Ihned vidíme, že

−∞

y (t ) je harmonický signál o stejném úhlovém kmitočtu jako u (t ) . Pro t → ∞ je

(

y (t ) = ybb (t ) = y max cos ω y t + θ y

)

(5.1.22)

kde ω y = ωu . Lineární soustava může měnit amplitudu a počáteční fázi harmonického signálu, nemůže ale měnit jeho kmitočet. To je další základní závěr. Poznamenejme, že symbolem ybb (t ) označíme dále ustálenou odezvu LTIL soustavy. Příklad 5.1.1

Výpočet odezvy při nulovém počátečním stavu Impulsová charakteristika h(t ) = 0.4 e −100t cos(800t ) LTIL soustavy je aproximována sekvencí MATLAB delta=0.00005; t=0:delta:.05; h=.4*exp(-100*t).*cos(800*t);

viz Obrázek 5.1.2 A. Vstupní signál u (t ) = 8 u (t ) = −8

t ≤ 0.016 t > 0.016

je definován příkazem MATLAB u=8*sign(0.016-t);

viz Obrázek 5.1.2 B. 79

Sekvence příkazů MATLAB yb=delta*conv(h,u); subplot (313) plot(0:delta:.1,yb) axis([0 .05 -.01 .01]) ylabel('output signál yb(t)') grid xlabel('t (s)')

počítá a zobrazuje konvoluci vstupního signálu a impulsové charakteristiky soustavy. Generuje odezvu soustavy při nulovém počátečním stavu. Výsledek je zobrazen na Obrázku 5.1.2 C.

Obrázek 5.1.2

Příklad 5.1.1 Výpočet odezvy při nulovém počátečním stavu yb (t )

80

Příklad 5.1.2

Výpočet odezvy při nulovém počátečním stavu yb (t ) Uvažujme soustavu na Obrázku 5.1.3. Soustava se nachází v nulovém počátečním stavu. x(0 − ) = y (0 − ) = 0 . Pro jednoduchost zvolíme x (t ) = y (t ) . Vstupem soustavy je napětí u (t ) . Výstupem je

napětí na kondenzátoru. Impulsovou charakteristiku h(t ) lze určit nepřímým měřením jako h (t ) = 0.5 e −0.5t . Situaci ukazuje Obrázek 5.1.3.

Příklad 5.1.2

Obrázek 5.1.3

Soustava prvého řádu

Jakmile je h(t ) známo, určíme odezvu při nulovém počátečním stavu yb (t ) pro libovolný vstup u (t ) jako t

yb (t ) =

∫ u(τ)0.5 e

0−

− 0.5 (t − τ )

t

dτ = ∫ u (t − τ ) 0.5 e − 0.5τ dτ

(5.1.23)

0−

Poznamenejme, že v praxi nelze Diracův impuls δ (t ) generovat. V praktických aplikacích můžeme například generovat vstup blízký integrálu z Diracova impulsu, kterým je Heavisideův skok η (t ) a měřit odezvu na η (t ) , kterou je přechodová charakteristika

g (t ) =

t

∫ h(t ) dt . Potom nalezneme h(t ) jako

0−

h (t ) =

d g (t ) dt

(5.1.24)

81

Odezva na jednotkový skok - přechodová charakteristika g (t )

Přechodová charakteristika g (t )

Definice 5.1.3

Přechodová charakteristika g (t ) je výstup LTIL soustavy generovaný vstupem η(t ) při nulovém počátečním stavu soustavy. x(0 − ) = 0 ⎫ ⎬ → g (t ) t ≥ 0 u (t ) = η(t ) ⎭

nebo, stručně,

(5.1.25)

η(t ) → g (t )

Porovnejme: přechodová charakteristika g (t ) v Příkladu 5.1.2 má tvar tvar g (t ) = 1 − e −0.5t . 5.2

Popis LTIL soustavy diferenciální rovnicí

Konvoluce (5.1.18) popisuje odezvu při nulovém počátečním stavu yb (t ) libovolné LTIL soustavy. S využitím vztahu

d dt

t

∫

f (t , τ ) dτ =

0−

t

∫

0−

d f (t , τ ) dτ + f (t , τ ) τ = t derivujeme dt

(5.1.18), y&b (t ) =

t

∫

0−

d u (τ)h(t − τ)dτ + u (τ)h(t − τ) τ = t dt

(5.2.1)

Lze ukázat, že popisuje-li h(t ) dynamiku LTIL soustavy n -tého řádu, potom n -tou derivací konvoluce získáme obyčejnou lineární diferenciální rovnici n -tého řádu s konstantními koeficienty. Algoritmus odpovídající (5.2.1) budeme ilustrovat na Příkladu 5.2.1. Získáme diferenciální rovnici (5.2.2) popisující odezvu yb (t ) soustavy při nulovém počátečním stavu (5.2.2)

(n −1)

yb(n ) (t ) + a n −1 yb

(t ) + L + a0 yb (t ) = bm u (m ) (t ) + bm −1u (m −1) (t ) + L + b0u (t ) ,

m ≤ n,

yb (0 − ) = 0 y& b (0 − ) = 0 M

yb (n −1) (0 − ) = 0

82

Tuto je třeba kombinovat s rovnicí, popisující odezvu y a (t ) soustavy při nulovém vstupním signálu

(n −1)

y a(n ) (t ) + a n −1 y a

(t ) + L + a0 y a (t ) = 0

(5.2.3)

y a (0 − ) = y 0 y& a (0 − ) = y1 M

y a (n −1) (0 − ) = y n −1 Úplný popis úlohy ve tvaru (5.2.4) je dán komposicí vztahů (5.2.2) a ( 5.2.3) pro celkovou odezvu y (t ) = y a (t ) + y b (t ) :

(5.2.4) y (n ) (t ) + an −1 y (n −1) (t ) + L + a0 y (t ) = bm u (m ) (t ) + bm −1u (m −1) (t ) + L + b0u (t ) , m ≤ n, y (0 − ) = y0 y& (0 − ) = y1 M

y (n −1) (0 − ) = yn −1 Poznamenejme, že podmínka m ≤ n je podmínkou nutnou k tomu, aby byla LTIL soustava kauzální.

Příklad 5.2.1

Odvození diferenciální rovnice Výstup soustavy H na Obrázku 5.2.1 je pro y (0 − ) = 0 a pro libovolný vstup u (t ) , viz Příklad 5.1.2, y (t ) = yb (t ) =

t

∫ u(τ)0.5 e

− 0.5 (t − τ )dτ

(5.2.5)

0−

83

Příklad 5.2.1

Obrázek 5.2.1

d S využitím vztahu dt

t

t

0−

0−

∫ f (t, τ) dτ = ∫

Soustava prvého řádu

d f (t , τ ) dτ + f (t , τ ) τ = t derivujeme (5.2.5) , dt (5.2.6)

y& (t ) =

t

∫ u(τ) 0.5(− 0.5) e

− 0.5 (t − τ )

dτ + u (τ ) 0.5 e − 0.5 (t − τ ) τ = t

0−

t

y& (t ) = −0.5 ∫ 0.5 e − 0.5 (t − τ )u (τ ) dτ + 0.5u (τ )

(5.2.7)

0−

Dosazením z (5.2.7) získáme y& (t ) = −0.5 y (t ) + 0.5u (t ) y (0 − ) = 0

(5.2.8)

y& (t ) + 0.5 y (t ) = 0.5u (t ) y (0 − ) = 0

(5.2.9)

tj.

Získali jsme obyčejnou lineární diferenciální rovnici prvního řádu s konstantními koeficienty popisující úlohu typu x (0 − ) = 0

u (t ) t ≥ 0

⎫ ⎬ ⎭

y (t ) t ≥ 0

(5.2.10)

Diferenciální rovnice popisující úlohu typu

84

x (0 − ) = y0 ⎫ ⎬ u (t ) t ≥ 0 ⎭

y (t ) t ≥ 0 ,

(5.2.11)

má tvar y& (t ) + 0.5 y (t ) = 0.5u (t ) y (0 − ) = y0

(5.2.12)

Existují dva základní přístupy k odvození diferenciální rovnice popisující LTIL soustavu. První přístup je založen na měření. Znalost vnitřní stuktury soustavy se obvykle nepředpokládá. Na soustavu přivedeme známý vstup a měříme výstupní odezvu soustavy při nulovém počátečním stavu soustavy. Ve zvláštním případě přivedeme aproximaci Diracova impulsu a měříme aproximaci impulsové charakteristiky. Poté z měrných dat odvodíme diferenciální rovnici soustavy. Druhý přístup je založen na apriorní znalosti vnitřní stuktury soustavy. Pro odvození diferenciální rovnice využijeme fyzikální zákony, jako jsou například Newtonův nebo Kirchhoffův zákon.

Příklad 5.2.2

Mechanická soustava druhého řádu Uvažujme soustavu na Obrázku 5.2.2. Hmota m tělesa je umístěna v jeho těžišti. Těleso leží na podlaze a je spojeno pružinou se stěnou. Vstup soustavy je síla u (t ) která na těleso působí. Výstup soustavy je vzdálenost y (t ) mezi těžištěm tělesa a jeho rovnovážnou polohou.

(t) Obrázek 5.2.2

Příklad 5.2.2. Soustava druhého řádu

Síla, kterou vyvolává lineárně působící pružina je − k1 y (t ) , kde k1 je tzv. konstanta pružiny. Tření mezi tělesem a podlahou lze přibližně popsat pomocí tří složek: statického tření, Coulombova tření a viskozního tření, viz Obrázek 5.2.3.

85

Obrázek 5.2.3

Příklad 5.2.2 Statické, Coulombovo a viskozní tření

Protože jsme omezili studium na lineární soustavy, můžeme z těchto tří složek dále pracovat pouze s viskozním třením. Sílu viskozního tření modelujeme výrazem − k2 y& (t ) , kde k2 je koeficient viskozního tření.

Newtonův zákon říká, že přivedená síla jednak překonává sílu pružiny a sílu tření, jednak urychluje hmotné těleso. u (t ) − k1 y (t ) − k 2 y& (t ) − m&y&(t ) = 0

Obrázek 5.2.4

(5.2.13)

Příklad 5.2.2. Odvození rovnic (5.2.13), (5.2.16)

Získali jsme obyčejnou lineární diferenciální rovnici druhého řádu s konstantními koeficienty popisující LTIL soustavu druhého řádu. Abychom získali standardní tvar rovnice, podělíme rovnici konstantou m a dosadíme 1 k k a0 = 1 , a1 = 2 and b0 = . m m m Získáme 86

&y&(t ) + a1 y& (t ) + a0 y (t ) = b0 u (t )

(5.2.14)

Aby bylo možno rovnici pro konkrétní případ řešit, je třeba znát počáteční podmínky jejího řešení, popisující počáteční stav soustavy. Zde například informaci o vzdálenosti a rychlosti tělesa v čase t = 0 − . y (0 − ) = y& (0 − ) =

y0 y1

(5.2.15)

Poté získáme celkový popis soustavy a úlohy:

&y&(t ) + a1 y& (t ) + a0 y (t ) = b0 u (t ) y (0− ) = y& (0− ) =

(5.2.16)

y0 y1

Příklad 5.2.3

Popis soustav s odlišnou fyzikální podstatou diferenciální rovnicí stejného typu Síla studované teorie je v tom, že ji lze použít obecně. Zde demonstrujme, jak lze použít diferenciální rovnice stejného typu pro popis soustav se zcela odlišnou fyzikální podstatou. Zvolme si pro demonstraci tři rozdílné soustavy prvního řádu. Obecný tvar diferenciální rovnice je a1 y& (t ) + a0 y (t ) = b0 u (t )

y (0 − ) = y 0

(5.2.17)

Elektrická soustava Elektrický RC filtr je nakreslen na Obrázku 5.2.5. Při běžných pracovních podmínkách můžeme tento obvod považovat za LTIL soustavu. Vstupem je napětí ui přivedené ze zdroje napěťového signálu. Výstup je napětí uo které lze měřit na kondenzátoru C .

87

Obrázek 5.2.5

Elektrická soustava prvního řádu

Platí 1 1 uo (t ) = ui (t ) RC RC uo (0 − ) = uo0 u&o (t ) +

(5.2.18)

kde resistanci R lze vyjádřit jako u − uo R= i i

(5.2.19)

kapacitanci C lze vyjádřit jako C=

dQ duo

(5.2.20)

kde i je elektrický proud obvodem a Q je elektrický náboj na kondenzátoru, Q (t ) = ∫

t

−∞

i (t )dt

(5.2.21)

RC filtr na Obrázku 5.2.5 je snad nejznámější soustava prvního řádu. Použijeme ji v dalším textu v příkladech Soustava s proměnnou výškou hladiny Soustava s proměnnou výškou hladiny je načrtnuta na Obrázku 5.2.6. Studujeme-li pouze malé změny stavových veličin od jejich ustálených hodnot, můžeme soustavu uvažovat jako lineární. Zde předpokládejme, že se jak vstup, kterým je průtočné množství qi , tak výstup, kterým je průtočné množství qo mění jen málo a tak, že se následně i stavová veličina soustavy, výška hladiny h , mění jen v okolí své ustálené hodnoty.

88

Obrázek 5.2.6

Soustava s proměnnou výškou hladiny prvního řádu

1 1 qo (t ) = qi (t ) RC RC qo (0 − ) = qo0 q&o (t ) +

(5.2.22)

Odpor R pro průtočné množství je definován poměrem průtočného množství a odpovídající změny výšky hladiny v nádobě

1 dqo = R dh

(5.2.23)

Kapacitance C nádoby je definovaná poměrem změny množství kapaliny uchovávané v nádobě ke změně výšky hladiny v nádobě C=

dQ dh

(5.2.24)

kde Q je množství kapaliny v nádobě. Rovnice (5.2.22) má tvar rovnice (5.2.18) , která popisuje RC filtr.

Mechanická soustava Mechanickou analogií elektrického RC filtru je soustava na Obrázku 5.2.7. Při běžných pracovních podmínkách můžeme tuto soustavu považovat za LTIL soustavu. Vstup je výchylka xi závěsu pružiny. Výstup je výchylka xo horní části tlumiče.

89

Obrázek 5.2.7

Mechanická soustava prvního řádu

Výchylka xo závisí na xi podle rovnice (5.2.25) b x&o (t ) + xo (t ) = xi (t ) k xo (0 − ) = xo0

(5.2.25)

kde k je konstanta pružiny k=

dFs dxi − dxo

(5.2.26)

b je koeficient viskosního tlumení tlumiče b=

dFd dx&o

(5.2.27)

Fs je síla působící na pružinu, Fd je síla působící na tlumič. Dosadíme-li, čistě formálně, v rovnici (5.2.25), za

b výraz RC , získáme k

1 1 xo (t ) = xi (t ) RC RC xo (0 − ) = xo0 x&o (t ) +

(5.2.28)

Rovnice (5.2.28) má opět tvar rovnice (5.2.18) , která popisuje RC filtr.

90

LINEÁRNÍ time-invariantní DYNAMICKOU soustavu se

soustředěnými parametry

n-tého řádu SE SPOJITÝM ČASEM, SE

ZADANÝM POČÁTEČNÍM STAVEM Lze popsat

obyčejnou LINEÁRNÍ

DIFERENCIÁLNÍ rovnicí n-tého řádu

s konstantními koeficienty , SE ZADANÝMI POČÁTEČNÍMI PODMÍNKAMI

Poznámka o řešení lineární diferenciální rovnice s konstantními koeficienty

V tomto kurzu se setkáváme s obyčejnými lineárními diferenciálními rovnicemi s konstantními koeficienty. Tato rovnice popisuje LTIL soustavu, její řešení je výstupní signál y (t ) soustavy. Analytický tvar y (t ) je obvykle možné, a žádoucí, nalézt tehdy, je-li vstupní

signál u (t ) standardní periodický signál nebo jiný standardní, například impulsní nebo exponenciální signál. Nicméně i v těchto, a prakticky ve všech ostatních případech, použijeme pro výpočet y (t ) počítač a numerické aproximační metody. Analytické řešení lineární diferenciální rovnice s konstantními koeficienty lze obecně najít jako konečný výraz tehdy, je-li její pravá strana daná konečným součtem součinů funkcí z množiny funkcí, jejíž derivace jsou funkce stejné třídy. V našem případě je u (t ) reálný signál. Z toho vyplývá, že výše uvedená množina obsahuje funkce konst ; t n , n = 0,1,2,L ; e at ; cos(at ); sin (at ); cosh (at ); sinh (at )

(5.2.29)

Klasický algoritmus pro nalezení analytického řešení diferenciální rovnice v tomto případě využívá metodu neurčitých koeficientů, kterou se nejdříve určí libovolné řešení rovnice nehomogenní y P ( x ) vyhovující pravé straně diferenciální rovnice bez ohledu na počáteční podmínky rovnice. Toto libovolné řešení se nazývaná partikulární integrál. U složitějších úloh je obvykle obtížné i v případech, kdy pravá strana rovnice, vstupní signál soustavy,

91

splňuje výše uvedenou podmínku, určit y P ( x ) bez konzultace s literaturou, respektive s tabulkami. Takovéto tabulky využívá pro řešení diferenciálních rovnic i metoda Laplaceovy transformace, pracující původně v praxi s tabulkami dvojic Laplaceovy transformace. Výhodou metody Laplaceovy transformace navíc je, že zjednodušuje zápis základních typů LTIL soustav. Laplaceova transformace je podporována i prostředím MATLAB. Z hlediska rozvoje počítačových metod nalezení analytického řešení obyčejné lineární diferenciální rovnice s konstantními koeficienty metoda Laplaceovy transformace postupně zastarává. K disposici jsou stále modernější algoritmy, například v MATLAB Maple software. V tomto textu se Laplaceově transformaci nevěnujeme. V následujícím odstavci ukážeme, na rovnici typu ( 2. 56 ), řešení obyčejné lineární diferenciální rovnice s konstantními koeficienty analytickou metodou neurčitých koeficientů. V dalším textu potom dáme přednost numerickým metodám obsaženým v MATLAB Simulink Toolbox a budeme vždy hledat celkovou odezvu soustavy y (t ) jako součet odezvy soustavy při nulovém vstupním signálu y a (t ) a odezvy soustavy při nulovém počátečním stavu y b (t ) , y (t ) = y a (t ) + yb (t )

Metoda neurčitých koeficientů

Úkolem je najít analytické řešení diferenciální rovnice, popisující RC filtr, na který je přiveden harmonický sinusový signál. V čase t = 0 − je na výstupu obvodu signál y 0 . Zvolme úhlový kmitočet vstupního signálu ω = 1 . Rovnice má pak bezrozměrný tvar

y& (t ) +

1 1 y (t ) = sin (t ) RC RC y (0 − ) = y 0

(5.2.30)

92

Partikulární řešení rovnice nehomogenní y (t ) najdeme určením konstant, zde jedné konstanty, v obecném řešení rovnice nehomogenní y N (t ) tak, aby toto odpovídalo počátečním podmínkám, zde jedné počáteční podmínce, diferenciální rovnice (5.2.30). Celkové řešení y N (t ) najdeme jako součet obecného integrálu rovnice homogenní y O (t ) a libovolného partikulárního integrálu rovnice nehomogenní y P (t ) . Obecný integrál rovnice homogenní y O (t ) najdeme pomocí řešení charakteristické rovnice. Obsahuje konstanty, které se přenesou do y N (t ) . Libovolný partikulární integrál rovnice nehomogenní y P (t ) najdeme metodou neurčitých koeficientů. 1.

Určíme obecný integrál homogenní diferenciální rovnice y O (t ) y& O (t ) +

1 y O (t ) = 0 RC

(5.2.31)

Řešení hledáme ve tvaru y O = Ke Qt kde Q = −

(5.2.32)

1 je kořen charakteristické rovnice RC Q+

1 =0 RC

(5.2.33)

k diferenciální rovnici (5.2.30). Konstanta K je zatím libovolná a bude použita v závěru výpočtu k určení partikulárního integrálu rovnice nehomogenní vyhovujícího počáteční podmínce řešení.

2.

Určíme libovolný partikulární integrál diferenciální rovnice nehomogenní y P (t ) Libovolný partikulární integrál rovnice nehomogenní y P (t ) nemusí vyhovovat počátečním podmínce řešení diferenciální rovnice. Libovolný partikulární integrál rovnice nehomogenní y P (t ) zde najdeme metodou neurčitých koeficientů. Nehomogenní rovnice je y& (t ) +

1 1 y (t ) = sin (t ) RC RC

(5.2.34)

93

Její řešení y P (t ) odhadneme ve tvaru y P (t ) = A sin (t ) + B cos(t )

(5.2.35)

y& P (t ) = A cos(t ) − B sin (t )

(5.2.36)

odtud

kde A a B jsou neurčité koeficienty. Získáváme tedy A cos(t ) − B sin (t ) +

1 [A sin (t ) + B cos(t )] = 1 sin (t ) RC RC

(5.2.37)

Koeficienty A a B určíme porovnáním koeficientů u vzájemně lineárně nezávislých funkcí sinus a kosinus sin (t ) ⇒

1 1 A= RC RC 1 A+ B= 0 RC

−B+

cos(t ) ⇒

(5.2.38)

Odtud A=

1

1 + R 2C 2 RC B=− 1 + R 2C 2 3.

(5.2.39)

Obecné řešení rovnice (5.2.30) má tedy tvar y N (t ) = Ke Qt + A sin (t ) + B cos(t )

(5.2.40)

kde koeficienty A a B jsou dány rovnicí (5.2.39) a exponent Q vztahem (5.2.33) 4.

Abychom nalezli partikulární integrál y (t ) diferenciální rovnice (5.2.30), je třeba ještě určit konstantu K . Určí se snadno řešením rovnice (5.2.40) pro t = 0 a počáteční podmínku diferenciální rovnice: K = y0 − B

(5.2.41) 94

Je tedy y (t ) = ( y 0 − B )e Qt + A sin (t ) + B cos(t )

(5.2.42)

Rozklad celkové odezvy soustavy y (t ) na součet y (t ) = ya (t ) + yb (t )

Je mnoho způsobů, jak algoritmizovat řešení úloh analýzy LTIL soustavy. Mnoho z nich přímo podporuje i pochopení chování soustavy. První takovýto základní algoritmus hledá celkovou odezvu y (t ) jako součet odezvy při nulovém vstupním signálu y a (t ) a odezvy při nulovém počátečním stavu y b (t ) , y (t ) = y a (t ) + yb (t )

(5.2.43)

Druhý takovýto základní algoritmus hledá celkovou odezvu y (t ) jako součet přechodné odezvy y aa (t ) a ustálené odezvy y bb (t )

y (t ) = y aa (t ) + ybb (t )

(5.2.44)

V těchto skriptech je dále dána přednost algoritmu, který hledá celkovou odezvu y (t ) jako součet odezvy při nulovém vstupním signálu y a (t ) a odezvy při nulovém počátečním stavu y b (t ) , y (t ) = y a (t ) + yb (t )

(5.2.45)

Příklad 5.2.4

Odezva při nulovém vstupním signálu y a (t ) a odezva při nulovém počátečním stavu y b (t ) soustavy Uvažujme LTIL soustavu se spojitým časem, která je popsána rovnicí ( 2. 93 )

kde

&y&(t ) + y& (t ) + 2 y (t ) = u (t )

(5.2.46)

u (t ) = 10 sin(5t )

(5.2.47)

Počáteční stav soustavy určuje počáteční podmínku diferenciální rovnice 95

y (0 − ) = 1 y ' (0 − ) = 1

(5.2.48)

Rovnici (5.2.46) lze rozložit následovně:

&y&a (t ) + y& a (t ) + 2 y a (t ) = 0 y a (0 − ) = 1 y& a (0 − ) = 1 &y&b (t ) + y& b (t ) + 2 yb (t ) = 10 sin(5t ) y b (0 − ) = 0 y& b (0 − ) = 0 y (t ) = y a (t ) + yb (t ) .

(5.2.49) (5.2.50)

(5.2.51) (5.2.52)

(5.2.53)

Blokové schéma modelu rovnice (5.2.53) v MATLAB Simulink Toolbox je na Obrázku 5.2.8. Z rovnice (5.2.49) obdržíme

&y&a (t ) = − y& a (t ) − 2 y a (t )

(5.2.54)

Rovnice (5.2.54) je modelována sumátorem, viz Obrázek 5.2.8. Integrací &y&a (t ) prováděnou integrátorem z počáteční hodnoty y& a (0 − ) , se generuje derivace y& a (t ) . Následnou integrací y& a (t ) prováděnou druhým integrátorem z počáteční hodnoty y a (0 − ) , se generuje funkce y a (t ) , viz Obrázek 5.2.8.

y& a (t ) = y a (t ) =

t

∫ &y&a (t )dt + y& a (0 − )

0− t

(5.2.55)

∫ y& a (t )dt + y a (0 − )

0−

Podobné schéma snadno odvodíme pro yb (t ) . Výsledky simulace jsou uvedeny na Obrázku 5.2.9.

96

y'a (0− )

ya (0− )

Obrázek 5.2.8 Blokové schéma, realizující stukturu diferenciální rovnice (5.2.46)

97

Obrázek 5.2.9

Diferenciální rovnice (5.2.46). MATLAB Simulink Toolbox. Výsledky simulace.

Řešení úlohy v MATLAB Symbolic Toolbox

Mějme diferenciální rovnici druhého řádu a&y& + by& + cy = g (t ) y (0) = y 0 y& (0) = y1

(5.2.56)

Rovnice může být v MATLAB Symbolic Toolbox řešena například opět jako součet odezvy při nulovém vstupním signálu a&y&a + by& a + cy a = 0 y a (0) = y a 0 y& a (0) = y a1

(5.2.57)

a odezvy při nulovém počátečním stavu a&y&b + by& b + cy b = g (t ) y b (0) = 0 y& b (0) = 0

(5.2.59)

Jádro programu pro řešení y a (t ) v MATLAB Symbolic toolbox software vyplývá přímo z tvaru rovnice: (5.2.60) ya=dsolve('a*D2ya+b*Dya+c*ya=0','ya(0)=ya0','Dya(0)=ya1'); disp('ya='),pretty(simplify(ya))

MATLAB Symbolic toolbox generuje výsledky a zobrazuje je ve analytickém tvaru na obrazovce, ya= 2 1/2 1/2 (%2 b y0 + %2 (b - 4 a c) y0 + 2 %2 y1 a - %1 b y0 2 1/2 / 2 1/2 + %1 (b - 4 a c) y0 - 2 %1 y1 a) / (b - 4 a c) / (5.2.61) 2

1/2

98

(b + (b - 4 a c) ) t %1 := exp(- 1/2 -----------------------) a 2 1/2 (b - (b - 4 a c) ) t %2 := exp(- 1/2 -----------------------) a

K disposici je příkaz ezplot(.), zde, např. ezplot(’ya’[0,10]), pro grafické zobrazení výsledků. Obdobně jádro programu pro odezvu při nulovém počátečním stavu: (5.2.62) yb=dsolve('a*D2yb+b*Dyb+c*yb=u','yb(0)=0','Dyb(0)=0'); disp('yb='),pretty(simplify(yb))

MATLAB Symbolic Toolbox generuje v tomto případě řešení yb (t ) , které je funkcí g (t ) . Na čtenáři ponecháváme, aby si výše uvedené algoritmy a příkazy v MATLAB Simulink Toolbox a MATLAB Symbolic Toolbox naprogramoval a výpočet celkové odezvy y (t ) = y a (t ) + yb (t ) pro konkrétní hodnoty koeficientů rovnice, počáteční podmínky a pravou stranu rovnice, sám provedl.

6

ANALÝZA LTIL SOUSTAVY V KMITOČTOVÉ OBLASTI

6.1

Popis LTIL soustavy kmitočtovým přenosem

Odvození H (ω) Fourierovou transformací konvoluce

Spektrum Y (ω) výstupního signálu y (t ) soustavy lze určit pomocí Fourierovy transformace obou členů, vstupujících do konvoluce (5.1.17). Pokud tyto Fourierovy transformace existují, potom podle (3.1.49)

F {u (τ ) * h(τ )} = U (ω)H (ω)

(3.1.49)

(6.1.1) 99

máme Y (ω) = U (ω)H (ω)

(6.1.2)

Odtud H (ω) =

Y (ω) U (ω)

(6.1.3)

H (ω) = F{ h(t ) } se nazývá kmitočtový přenos, nebo též kmitočtová charakteristika, soustavy.

Definice 6.1.1

Kmitočtový přenos, kmitočtová charakteristika soustavy, H (ω) Kmitočtový přenos, kmitočtová charakteristika soustavy, H (ω) LTIL soustavy je definován, pokud existuje, jako H (ω) = F{h(t )}

(6.1.4)

nebo H (ω) =

Y (ω) U (ω)

(6.1.5)

Znamená to, že impulsová charakteristika a kmitočtový přenos jsou transformační dvojice Fourierovy transformace. Nutnou podmínkou je, aby existovala Fourierova transformace funkce h(t ) . Fourierova transformace funkce y (t ) existuje, jestliže

∞

∫ y(t ) dt < K , K < ∞ .

Kmitočtový

−∞

přenos H (ω) = F{h(t )} proto existuje pro každou asymptoticky stabilní LTIL soustavu. Kmitočtový přenos H (ω) je, obecně, komplexní veličina a může být zapsán v polárním tvaru H (ω) = H (ω) e jφ (ω)

(6.1.6)

kde amplitudový kmitočtový přenos, amplitudová kmitočtová charakteristika soustavy, H (ω) = Im{H (ω)}2 + Re{H (ω)}2

(6.1.7)

a fázový kmitočtový přenos, fázová kmitočtová charakteristika soustavy, ⎡ Im{H(ω)}⎤ φ (ω) = arctan ⎢ ⎥ ⎣ Re{H(ω)}⎦

(6.1.8)

100

Dále, protože h(t ) je reálná funkce času, je H (ω) sudá funkce kmitočtu a φ (ω) je lichá funkce kmitočtu. Odvození H (ω) Fourierovou transformací diferenciální rovnice

Kmitočtový přenos H (ω) LTIL soustavy lze rovněž najít Fourierovou transformací diferenciální rovnice, která soustavu popisuje. Pro názornost zvolíme zvláštní případ LTIL soustavy n-tého řádu (5.2.4) , m = 1, n = 2 ,

&y&(t ) + a1 y& (t ) + a0 y (t ) = b1u& (t ) + b0u (t )

(6.1.9)

y (0 − ) = y 0 y& (0 − ) = y1 Pro nalezení H (ω) transformujeme rovnici Fourierovou transformací. Integraci rovnice provádíme v časovém intervalu (− ∞ , ∞ ) . Pro integraci je tedy požadována znalost u (t ) a

y(t ) i pro t < 0 . Například volba u (t ) = 0 t < 0 vede k y (t ) = 0 t < 0 a splňuje tento požadavek. Důsledkem je nulový počáteční stav soustavy, tj. nulové počáteční podmínky řešení odpovídající rovnice (6.1.9), y (0 − ) = 0 , y& (0 − ) = 0 . Pro periodický vstupní signál můžeme počátek posunout do t = −∞ a počáteční podmínky neuvažovat. Transformací

&y&(t ) + a1 y& (t ) + a0 y (t ) = b1u& (t ) + b0 u (t )

(6.1.10)

obdržíme

( jω)2 Y (ω) + a1 ( jω)Y (ω) + a0Y (ω) = b1 ( jω)U (ω) + b0U (ω) (6.1.11) což lze přerovnat do tvaru ⎤ ⎡ ⎢( jω)2 + a ( jω) + a ⎥Y (ω) = [b ( jω) + b ]U (ω) 1 0 1 0 ⎥ ⎢ ⎦ ⎣

(6.1.12)

neboli H (ω) =

b1 ( jω) + b0 Y (ω) = U (ω) ( jω)2 + a1 ( jω) + a0

(6.1.13)

101

Dělíme čitatel a jmenovatel konstantou a0 a obdržíme standardní tvar H (ω) ,

H(ω) =

β1( jω) + β0 Y(ω) = U(ω) α2( jω)2 + α1( jω) + 1

(6.1.14)

Přenos LTIL soustavy lze zjistit měřením použitím testovacího harmonického vstupního signálu o známé amplitudě a kmitočtu. Například, jestliže u (t ) = u max cos(ω0 t )

(6.1.15)

potom ustálený výstup soustavy je y (t ) = u max H (ω 0 ) cos[ω 0 t + φ (ω 0 ) ]

(6.1.16)

kde amplituda a fáze výstupu může být určena pomocí nějakého vektorového měřícího přístroje, z nichž nejstarší a nejjednodušší je osciloskop.

Příklad 6.1.1

Trajektorie H (ω) v komplexní rovině Soustava H prvního řádu je dána rovnicí y& (t ) + 2 y (t ) = 3u (t ) .

(6.1.17)

Úkolem je určit průběh kmitočtového přenosu H (ω) , průběh, tvar, kmitočtové charakteristiky H (ω) soustavy (6.1.17) v komplexní rovině, neboli tzv. trajektorii H (ω) v komplexní rovině.

Fourierovou transformací rovnice (6.1.17) získáme jω Y (ω) + 2Y (ω) = 3U (ω)

(6.1.18)

Odtud H (ω) =

3 1.5 = jω + 2 0.5 jω + 1

(6.1.19)

Nyní použijeme pěkný algoritmus: rozložíme H (ω) na reálnou a imaginární část a

102

ukážeme, že hledaný průběh kmitočtové charakteristiky H (ω) =

b0 K = a 0 jω + 1 j ω T + 1

(6.1.20)

je v komplexní rovině tvořen kružnicí. 1.5 0.5 jω + 1 (1 − 0.5 jω) = 1.5 (1 + 0.5 jω)(1 − 0.5 jω) 1 − 0.5 jω = 1.5 1 − 0.5 jω + 0.5 jω + 0.25ω 2 1 − 0.5 jω = 1.5 1 + 0.25ω 2

H (ω) =

(6.1.21)

Nyní ukážeme, že pro výraz (6.1.22), tento výraz je blízký, ale není totožný s (6.1.21), platí ⎛ 1 ⎜⎜ − 0.5 ⎝ 1 + 0.25ω 2

2

2

⎞ ⎛ 0.5ω ⎞ ⎟⎟ + ⎜⎜ ⎟ = 0.5 2 2⎟ ⎠ ⎝ 1 + 0.25ω ⎠

(6.1.22)

Skutečně 2

2 2⎞ ⎛ + ω 1 1 0 . 25 0 . 5 ω ⎛ ⎞ ⎜ ⎟ +⎜ − 0.5 ⎟ 2⎟ 2 ⎜ 1 + 0.25ω2 + ω 1 0 . 25 1 + 0 . 25 ω ⎝ ⎠ ⎝ ⎠

=

(

0.25 − 0.125ω2 + 0.125ω2

= 0.25

(1 + 0.25ω2 )2

2

(

1 + 0.5ω + 0.25ω

(1 + 0.25ω2 )2

)

2 2

)2 + 0.25ω2 (6.1.23)

= 0 .5 2

Tím je potvrzeno, že součet čtverců nad reálnou a imaginární částí výrazu (6.1.22) je roven konstantě, z čehož podle Pythagorovy věty vyplývá, že průběh kmitočtové charakteristiky (6.1.19)

103

H (ω) =

1.5 0.5 jω + 1

(6.1.19)

(6.1.24)

je v komplexní rovině tvořen kružnicí zobrazenou na Obrázku 6.1.1. Kmitočtovou charakteristiku (6.1.24) lze snadno vypočítat programem MATLAB omega=-100:0.01:100; H=1.5./(j*omega*0.5+1); X=real(H); Y=imag(H); plot(X,Y) xlabel('Re[H(omega)]') ylabel('Im[H(omega)]') axis([-2 2 -2 2]) axis square grid

Obrázek 6.1.1

Příklad 6.1.1 Trajektorie H (ω) =

1.5 v komplexní rovině 0.5 jω + 1

Výsledek lze zobecnit pro každou stabilní LTIL soustavu prvního řádu, viz Obrázek 6.1.2 .

104

Obrázek 6.1.2

Trajektorie H (ω) =

K v komplexní rovině jω T + 1

MATLAB nabízí mnoho dalších prostředků pro výpočet a grafický záznam kmitočtového přenosu H (ω) . Zde je jeden takovýto jednoduchý program, využívající pro určení amplitudové kmitočtové charakteristiky H (ω) a fázové kmitočtové charakteristiky ⎡ Im{H(ω)}⎤ φ (ω) = arctan ⎢ ⎥ kmitočtového přenosu H (ω) MATLAB funkci bode( ). ⎣ Re{H(ω)}⎦ num=[1.5]; den=[0.5 1]; omega=0:0.01:20; [mag,phase]=bode(num,den,omega); subplot(211) plot(omega,mag) ylabel('mag[H(omega)]') xlabel('omega') axis([0 20 0 1.5]) subplot(212) plot(omega,phase) ylabel('phase[H(omega)],degrees') xlabel('omega') axis([0 20 -90 0])

105

Obrázek 6.1.3, Bodeho diagram, ukazuje výslednou amplitudovou a fázovou charakteristiku kmitočtového přenosu H (ω) =

Obrázek 6.1.3

6.2

1.5 . Obrázek 6.1.3 odpovídá Obrázku 6.1.1. 0.5 jω + 1

Příklad 6.1.1. Kmitočtové charakteristiky soustavy H (ω) =

1.5 . 0.5 jω + 1

Transformace signálu LTIL soustavou

Lineární time-invariantní soustava s konstantními koeficienty může měnit spektrum vstupního signálu. Na daném kmitočtu může měnit amplitudu a počáteční fázi harmonické složky spektra, nemůže měnit úhlový kmitočet této složky. Příklad 6.2.1

Zkreslení signálu soustavou Vraťme se k časovému sledu pravoúhlých impulsů popsanému v Příkladu 3.2.1, vztah (3.2.33).

106

w(t ) = 0.5 +

1 π 2 cos t π 10

+

−1 π 2 cos 3 t 3π 10

+

1 π 2 cos 5 t 5π 10 (3.2.33)

+ L (6.2.1)

Předpokládejme, že je tento signál zkreslován soustavou, která ho přenáší ze svého vstupu na svůj výstup. Pro iluslustraci zkreslení signálu dynamickou soustavou zvolme nejjednodušší případ, kdy signál prochází LTIL soustavou prvního řádu. Takovouto soustavou může být například RC filtr popsaný rovnicí (5.2.18), jehož vstupním signálem je napětí získávané z napěťového zdroje a výstupním signálem napětí na zatěžovacím odporu o nekonečně velké rezistenci. Pro nalezení kmitočtového přenosu RC filtru (5.2.18) předpokládáme nulový počáteční stav soustavy, y (0 − ) = 0 . Rovnici (5.2.18) transformujeme Fourierovou transformací. Protože je transformovatelná, získáme 1 RCjω + 1 K = Tjω + 1

H(ω) =

(6.2.2)

Abychom mohli v dalším textu uvést kvantitativní výsledky, zvolme konkrétní jednoduchou soustavu prvního řádu. Zvolme K = 1 a T = 1 , neboli RC = 1 .

Získáme (6.2.3) H (ω) =

−ω 1 1 = + j = Re{H (ω)} + j Im{H (ω)} = H (ω) e jφ(ω) 2 jω + 1 1 + ω 2 1+ ω

Vypočítejme kvantitativní hodnoty částí H (ω) podle (6.2.3) pro úhlové kmitočty, odpovídající harmonickým složkám časového sledu pravoúhlých impulsů. Tabulka 6.2.1 shrnuje výsledky pro konstantní složku a pro prvních pět harmonických složek.

107

Index

Úhlový

m

kmitočet

0 1 2 3 4 5

0 π/10 2π/10 3π/10 4π/10 5π/10

Tabulka 6.2.1

Re{H (ω)} Im{H (ω)}

1 0.2884 0.0920 0.0431 0.0247 0.0160

0 -0.4530 -0.2890 -0.2031 -0.1552 -0.1253

H (ω)

1 0.5370 0.3033 0.2076 0.1572 0.1263

φ(ω)

0 -1.0039 -1.2626 -1.3617 -1.4130 -1.4438

Příklad 6.2.1 Zesílení harmonických složek vstupního sledu pravoúhlých impulsů (6.2.1) soustavou (6.2.3)

Je-li tedy vstupním signálem zvoleného RC filtru signál u (t ) daný časovým sledem pravoúhlých impulsů (6.2.1), je výstupním signálem RC filtru signál y (t ) časový sled přibližně exponenciálně tvarovaných impulsů (6.2.4) y(t ) = 0.5 +

0.5370 ⎛π ⎞ 0.2076 ⎛ π ⎞ 0.1263 ⎛ π ⎞ 2 cos⎜ 5 t − 1.4438⎟ + L 2 cos⎜ 3 t − 1.3617⎟ + 2 cos⎜ t − 1.0039⎟ − π 5π 3π ⎝ 10 ⎠ ⎝ 10 ⎠ ⎝ 10 ⎠

Časový průběh výstupního signálu y (5) (t ) vypočítaný ze vztahu (6.2.4) dosazením pro konstantní složku a prvních pět harmonických složek je zobrazený, spolu se vstupním signálem u (t ) , na Obrázku 6.2.1.

Obrázek 6.2.1

Příklad 6.2.1

Časový průběh výstupního signálu y (5) (t )

Amplitudová kmitočtová charakteristika H (ω) se občas udává v decibelech, označuje se jako H (ω) dB a je definovaná vztahem 108

H (ω) dB = 20 log10 H (ω)

(6.2.5)

Poznamenejme, že H (ω) dB < 0 dB pro

H (ω) dB = 0 dB pro H (ω) dB > 0 dB pro

H (ω) < 1 H (ω) = 1

(6.2.6)

H (ω) > 1

Grafy H (ω) dB a φ(ω) , kde se ω vynáší v logaritmické stupnici, se nazývají Bodeho diagramy soustavy. Bodeho diagramy se často stanovují experimentálně měřením ustálené odezvy soustavy na budící harmonický vstupní signál u (t ) = U cos(ω0t ) s tím, že se měření provádí pro různé hodnoty ω0 . V MATLABu lze Bodeho diagramy vypočítat přímo z přenosu H (s ) nebo z kmitočtového přenosu H (ω) pomocí příkazu bode. Pouhý příkaz bode kreslí automaticky Bodeho diagramy. Pro další úpravu Bodeho diagramů lze pak použít různá options, viz Příklad 6.2.2 níže.

Příklad 6.2.2

Bodeho diagramy

Uvažujme soustavu druhého řádu s kmitočtovým přenosem H (ω) =

200( jω + 5) ( jω + 10)( jω + 50)

(6.2.7)

Program v MATLABu num=200*[1 5]; den=conv([1 10],[1 50]); omega=1:0.1:1000; [mag,phase]=bode(num,den,omega); magnitude=20*log10(mag); subplot(211) semilogx(omega,magnitude) title('Bode Diagram') axis([1 1000 -20 20]) ylabel('magnitude(dB)') subplot(212) semilogx(omega,phase) axis([1 1000 -90 45]) ylabel('phase(deg)') grid xlabel('omega(rad/s)')

109

generuje Bodeho diagramy soustavy (6.2.7). Diagramy jsou zobrazené na Obrázku 6.2.2.

Obrázek 6.2.2

Příklad 6.2.2. Bodeho diagramy soustavy (6.2.7)

LITERATURA 1. Edward W. Kamen., Bonnie S. Heck.: Fundamentals of Signals and Systems. Prentice Hall Inc., New Jersey 2000. 2. Chi-Tsong Chen: System and Signal Analysis. Saunders College Publishing, 1994. 3. Nevřiva P.: Analýza signálů a soustav. BEN Prana, Praha 2000. 4. MATLAB documentaticon, Math-Works Inc. Natick 2000 - 2004.

110