Fourierova transformace
Jean Baptiste Joseph Fourier (1768-1830)
Jeho obdivovatel (nedatováno)
Opáčko harmonických signálů Spojitý harmonický signál
= x ( t ) C1 cos (ω1t + ϕ1 )
()
4 cos ( 3t − 0.5 )
x t Příklad:= 4
C1cos(φ1)
3
C1 amplituda
2
ω1
1
ϕ1 f1 T1
[ rad / s ] počáteční fáze [ rad ] úhlová frekvence
skutečná frekvence 2π f1 = ω1 1 2π T = = základní perioda 1 f1 ω1
0 -1 -2 -3 -4 -0.5
0
0.5
1
1.5
2
2.5
3
3.5
4
T1=1/f1=2π/ω1
C1 = 4, ω1 = 3, ϕ1 = −0.5 Chyták: Proč jsou tam indexy 1 ? Cvičení: vyrobte několik harmonických signálů s různou úhlovou frekvencí a fází. Pozorujte data, ať vám „zalezou pod kůži“.
Diskrétní harmonický signál
= x [ n ] C1 cos (ω1n + ϕ1 )
C1 amplituda
ω1 ϕ1
úhlová frekvence [ rad ] počáteční fáze [ rad ]
2π
N1 základní perioda: pozor !!! N1 ≠ ω1 Musí platit že
cos (ω1 ( n + N1 ) ) = cos (ω1n )
Tedy rozdíl argumentů funkce cos musí být násobek základní periody funkce cos, což je 2π.
ω1 ( n + N1 ) − ω1n = ω1 N1 = k 2π
Příklad
x [ n ] = cos ( 8π n / 31) , tedy ω1 = 8π / 31
Dosadím do vzorečku
8π N1 = k 2π , tedy 4 N1 = 31k a mám 31
řešení
4 = k , N1 = 31 je základní perioda
1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 -5
0
5
10
15
N1=31
20
25
30
35
Proč je potřeba vědět hodně o harmonických funkcích ? Protože Fourierova dekompozice rozdělí signál na hromadu sinusovek a cosinusovek. Jde to vůbec ? • U diskrétních funkcí vždycky • U spojitých funkcí nejdou „ostré rohy“, ale můžeme se k nim hooodně přiblížit Takže to jde. A k čemu to je ? Proč vůbec rozděluji signál na hromadu harmonických funkcí ? Protože mi to o signálu hodně řekne !
Fourierova transformace a signály Signálů je spousta druhů (to už víte) a podle toho existují různé druhy Fourierovy transformace:
Cože ? Tolik různých Fourierů ? (a to ještě existují další verze, např. reálná DFT, komplexní DFT, FFT, …).
Reálně jsme schopni spočítat pouze DFT
Reálná DFT - základní myšlenka (stále stejný trik jako u konvoluce) Signál rozložím na hromadu cosinusovek a sinusovek. • Kolik ? Je-li vstupní signál N (od 0 do N-1), pak na N/2+1 sinusovek a N/2+1 cosinusovek. • Jakých ? Tyto harmonické funkce mají postupně vyšší a vyšší frekvence, a to tak šikovně, že se do daného signálu „vejdou“ právě 0x, 1x, 2x, …. U harmonických funkcí s takto danou frekvencí se mění pouze amplituda. • K čemu to ? Amplitudy nám říkají, „jak moc“ dané frekvence je v původním signálu • Co je vlastně výsledkem ? Amplitudy harmonických funkcí (frekvence jsou dané) • Sedí to ? Měl jsem N vzorků, teď mám N+2 vzorků ! Jak to ??? (první a poslední sinusovka má vždy amplitudu 0, takže ve skutečnosti to sedí) Příklad: Mám signál o 16 vzorcích (0-15). Rozložím ho tedy na 9 cosinusovek a 9 sinusovek. Ty mají stejný počet prvků jako původní signál. Když všechny harmonické sečtu, dostanu zpátky původní signál. Chyták: jakou jednotku má x-ová osa u harmonických funkcí ?
cos
sin
Původní signál (časová oblast)
Reálná / komplexní DFT Reálná DFT – pochopí i dítě ReX a ImX jsou jen označení ! Z N vzorků dostanu N+2 vzorků Co je na vodorovné ose ? Co je na svislé ose ?
Komplexní DFT – tohle je v Matlabu Všechny prvky jsou komplexní č. Z N vzorků dostanu N vzorků Matematicky elegantní Poznámka – obrázek není úplně správně, poznáte proč?
Dopředná DFT (reálná) Jak to spočítám ? N −1
ReX [ k ] = ∑ x [ j ] cos ( 2π kj / N ) j =0
N −1
ImX [ k ] = −∑ x [ j ] sin ( 2π kj / N ) j =0
Slovně: každý vzorek ve frekvenční oblasti najdu tak, že vynásobím časový signál sinusovkou nebo kosinusovou kterou hledám a všechno sečtu. Podivnost: mínus ve druhém vztahu – jen kvůli kompatibilitě s komplexní DFT.
Cvičení: naprogramujte reálnou DFT
Co jsme to dostali ? Konstanta casova oblast - konstanta 21 20.5 20 19.5 19
0
200 400 cas [sec] ReX
600 -11
10
15000 10000
ImX
x 10
5
5000 0
0 -5000
0
100
200
300
Co je to za podivnosti v ImX ?
-5
0
100
200
300
Sinusovka
casova oblast - sinusovka
casova oblast - sinusovka 1
1
0.5
0.5
0
0
-0.5
-0.5
-1
0
0.5
1.5 1 cas [sec] ReX
-1
2
0
0.5
ImX 100
8 6
0
4
1 1.5 cas [sec] ReX
2 ImX
40
50
30
0
20
-50
10
-100
0
-150
-100 2 -200
0 -2
0
50
100
150
200
-300
0
50
100
150
200
-10
0
50
100
150
200
-200
0
50
100
150
200
Dvě sinusovky casova oblast - dve sinusovky 2 1 0 -1 -2
0
0.5
1 1.5 cas [sec] ReX
2 ImX
40
50
30
0
20
-50
10
-100
0
-150
-10
0
50
100
150
200
-200
0
50
100
150
200
Co je na vodorovné ose ? • Číslo vzorku (sample number) [-] o to mají rádi programátoři o rozsah 0-N/2 • Podíl vzorkovací frekvence (fraction of sampling rate) f [-] o Nejvyšší frekvence je polovina vzorkovací frekvence (Nyquistův teorém) o Rozsah 0-0,5 • Vlastní frekvence (natural frequency) ω [rad] o Rozsah 0-π • Frekvence [Hz] o Nejlépe pochopitelné (čísla co dávají konkrétní smysl) o Vztaženo ke konkrétní vzorkovací frekvenci o Rozsah DC-polovina vzorkovací frekvence Příklad: Jak napsat kosinusovou ?
c [ n ] = cos ( 2π kn / N ) c [ n ] = cos ( 2π fn ) c [ n ] = cos (ω n )
Polární zobrazení Citace: když se budete snažit pochopit pravoúhlé zobrazení, vybuchne vám hlava Jak to zobrazit jinak ? Dvojici v pravoúhlých souřadnicích (aktuální vzorek ve frekvenční oblasti) si představím jako vektor. Převedením do polárních souřadnic: z dvojice (Re,Im) dostanu opět dvojici: • Magnituda (modul, absolutní hodnota, …) • Fáze (argument, úhel, …)
Polární zobrazení Převod pravoúhlé / polární MagX = [k ]
Převod polární / pravoúhlé
ReX [ k ] + ImX [ k ] 2
2
ImX [ k ] FazeX [ k ] = arc tan ReX [ k ]
ReX [ k ] = MagX [ k ] cos ( FazeX [ k ]) ImX [ k ] = MagX [ k ] sin ( FazeX [ k ])
Proč je polární zobrazení tak boží ? Hned vidím, jaké frekvence tam jsou a jaké ne, jak se mi mění fáze, je to „human readable“.
Cvičení: naprogramujte si převodník z pravoúhlého zobrazení DFT na polární
Poznámky k fázi 1. bacha na dělení nulou (Im/Re) 2. bacha znaménka u arctan (jinak bude fáze jen v rozmezí
−
π π
, 2 2
if (ReX[k] < 0) AND (ImX[k] < 0) THEN Faze[k]=Faze[k] - PI if (ReX[k] < 0) AND (ImX[k] >= 0) THEN Faze[k]=Faze[k] + PI 3. malé hodnoty magnitudy → nesmyslné hodnoty fáze
)
4. Nespojitost fáze Sinusovky jsou periodické, a fázový posuv o ϕ je stejný jako ϕ + 2π , ϕ + 4π ,... Takže možností kolik fáze vlastně přesně je, je nekonečno. Při výpočtu se vybere vždy ta nejmenší. Pro představu může pomoct rozbalená fáze na obrázku
5. Magnituda je vždy kladná Na obrázku máme příklad. Reálná část je pěkně hladká, imaginární je nulová. V polárním zobrazení ale máme ostré rohy, fuj fuj. Je tomu tak proto, že magnituda je vždy kladná (z definice). Matematicky je to ook, ale interpretovat to někdy může být obtížné. 6. Rychlohops ve fázi mezi −π a π protože fáze pro −π a π je stejná, zaokrouhlovací chyby mohou způsobit hopsání fáze, viz. obr.
Inverzní DFT (reálná) Mám hotovou dopřednou DFT. Jak se dostanu zpět na původní signál ?
= x[ j]
N /2
N /2
∑ ReX [ k ] cos ( 2π kj / N ) + ∑ ImX [ k ] sin ( 2π kj / N )
k 0= k 0 =
Kde
ReX [ k ] =
ReX [ k ]
ImX [ k ] = − Ale (!)
ReX [ 0] = ReX [ N / 2] = Proč ??? (další stránka)
N /2 ImX [ k ] N /2
ReX [ 0] N ReX [ N / 2] N
Cvičení: Naprogramujte inverzní DFT. Vyzkoušejte dopřednou a následně zpětnou transformaci.
Proč musíme čarovat s prvním a posledním koeficientem ?
Vše je jasné z obrázku . Šířka pásma pro jednotlivé vzorky není stejná. Co to vůbec je šířka pásma ? Jaký „rozsah frekvencí“ pokrývá náš jediný vzorek.
Jak se to dělá v Matlabu? Pozor na rozdíl mezi reálnou DFT (to co klovete) a komplexní DFT (Matlab). clear all % generování dat t = 0:0.01:1.99; % tedy vzorkujeme s periodou 0.01 sec, (100Hz) % vyrobim si dilci signaly o frekvencich 4 a 9 Hz y1 = 2*sin(t*4*2*pi); y2 = 4 * sin(t*9*2*pi); y = y1 + y2; % sectu je a mam dve namixovane sinusovky Y = fft(y); % provedu FFT % spoctu si magnitudu a fazi - prevedu do polarniho zobrazeni magnituda = abs(Y); faze = phase(Y); % ted to vsechno nakreslime subplot(3,1,1); plot(t,y); ylabel('amplituda'); xlabel('cas [sec]'); title('puvodni signal v casove oblasti'); subplot(3,1,2); plot(magnituda); ylabel('amplituda'); xlabel('vzorek [-]'); title('signal ve frekvencni oblasti'); subplot(3,1,3); % vezmeme jen pulku frekvenci! pulmagnituda = magnituda(1:length(magnituda)/2); % na x osu vyrobime data od 0 do pulky vzorkovacky f = linspace(0,50,length(pulmagnituda)); plot(f,pulmagnituda); ylabel('amplituda'); xlabel('frekvence [Hz]'); title('signal ve frekvencni oblasti');
