Digitální zpracování signálů Fourierova transformace, FFT Frekvenční analýza
3. přednáška
Jean Baptiste Joseph Fourier (1768 - 1830) Základy experimentální mechaniky
1
Frekvenční analýza Proč se frekvenční analýza provádí:
B
C a(t)
A(f)
A A
B
CD
E
čas
E
frekvence
vibrace
D
3. přednáška
Jednotlivé jevy jsou v časové oblasti promíchány
Základy experimentální mechaniky
Jednotlivé jevy jsou ve frekvenční oblasti od sebe odděleny
2/32
Digitální zpracování signálů Fourierova transformace (pro spojitou funkci): funkce x(t), periodická v čase T, může být vyjádřena jako nekonečná posloupnost:
a0 ∞ 2πnt 2πnt x ( t ) = + ∑ a n ⋅ cos + b ⋅ sin n 2 n =1 T T kde an a bn mohou být vypočteny ze znalosti x(t) pomocí vztahů
2 T 2πnt an = ⋅ ∫ x( t ) ⋅ cos dt 0 T T
2 T 2πnt bn = ⋅ ∫ x( t ) ⋅ sin dt 0 T T 3. přednáška
Základy experimentální mechaniky
2π = ω T
3/32
Digitální zpracování signálů Fourierova transformace (pro diskretizovanou funkci): funkce x(t) je diskretizována a trvá konečný čas T, je definována na množině N jednotlivých časových okamžiků tk:
a0 N / 2 2πnt k 2πnt k x k (= x(t k )) = + ∑ an ⋅ cos + bn ⋅ sin 2 n=1 T T
; k = 1,N
Koeficienty an a bn jsou Fourierovy neboli spektrální koeficienty funkce x(t) a často jsou zobrazovány ve tvaru amplitudy a fáze:
4 ⋅ c n ⋅ c − n = a 2n + b 2n
b φn = arctg − n an
tedy:
2πnt k c n ⋅ cos + ϕ n ; k = 1, N T n =− N / 2 N/2
x k (= x (t k )) = ∑
3. přednáška
Základy experimentální mechaniky
4/32
Diskrétní Fourierova transformace - DFT Vstupní signál je A/D převodníkem digitalizován a zaznamenán jako množina N diskrétních hodnot s pravidelnými časovými rozestupy v intervalu T. Předpokládá se, že vzorek v čase T je periodický. Je vypočtena konečná Fourierova řada (transformace) jako odhad požadované Fourierovy transformace. Platí základní vztah mezi délkou vzorku T, počtem diskrétních hodnot N, vzorkovací (digitalizační) frekvencí fs a rozsahem a rozlišením frekvenčního spektra (fmax, ∆f):
f 1 N fmax = s = ⋅ 2 2 T fS 1 ∆f = = N T 3. přednáška
fs … vzorkovací frekvence fmax … Nyquistova frekvence ∆f … frekvenční rozlišení rozsah získaného spektra je <0;fmax>
Základy experimentální mechaniky
5/32
Diskrétní Fourierova transformace - DFT Základní rovnice pro určení spektrálního složení: sin(2π / T ) x1 0.5 cos(2π / T ) x 0.5 cos(4π / T ) sin(4π / T ) 2 x 3 = 0.5 cos(6π / T ) sin(6π / T ) M 0.5 M M x N 0.5 cos(2Nπ / T ) sin(2Nπ / T )
K a0 K a1 K ⋅ b1 K M K M
≡
{xk } = [C]⋅ {an }
K určení neznámých spektrálních (Fourierových) koeficientů obsažených v {an} tedy použijeme:
{an } = [C]−1{x k } Nejpoužívanějším algoritmem výpočtu spektrální analýzy je rychlá Fourierova transformace - FFT. Tato metoda vyžaduje, aby N bylo celočíselnou mocninou 2. 3. přednáška
Základy experimentální mechaniky
6/32
Specifické rysy DFT Digitální Fourierova analýza má mnoho rysů, které mohou vést k chybným výsledkům, pokud nejsou správně ošetřeny. Jsou důsledkem – diskretizace – nutnosti omezit délku časového signálu Je třeba uvážit: – aliasing – chyby únikem – vliv oken – filtrování – frekvenční lupu – průměrování
3. přednáška
Základy experimentální mechaniky
7/32
Aliasing Tento jev plyne z diskretizace původně spojitého časového signálu. Při malé vzorkovací frekvenci je přítomnost vysokých frekvencí v původním signálu při tomto diskretizačním procesu špatně interpretována. Tyto vysoké frekvence se ve spektru objeví jako nízké frekvence, nebo spíše budou od skutečných nízkofrekvenčních složek nerozpoznatelné.
nízkofrekvenční signál vysokofrekvenční signál
Spektrum získané pomocí DFT je zkreslené, i když výpočet je proveden přesně!
3. přednáška
Základy experimentální mechaniky
8/32
Použití anti-aliasingového filtru Nejvyšší frekvence, která může být ve spektru obsažena, je ωs/2. Vyšší frekvence jsou “zrcadleny“ do nižších frekvencí.
skutečné spektrum signálu
spektrum získané z DFT
3. přednáška
Základy experimentální mechaniky
9/32
Použití anti-aliasingového filtru Anti-aliasingový filtr podrobí původní časový signál nízkopásmovému filtru s ostrou sestupnou hranou. nefiltrovaný signál
anti-aliasingový filtr
filtrovaný signál Protože filtry mají konečný sklon sestupné hrany, odstraňují se i spektrální měření ve frekvenčním rozsahu blízkém Nyquistově frekvenci ωs/2. Proto při 2048 bodové transformaci není výsledkem úplné 1024 čárové spektrum, ale typicky se zobrazuje pouze prvních 800 čar.
Anti-aliasingová opatření tvoří nedílnou součást analyzátoru! 3. přednáška
Základy experimentální mechaniky
10/32
Chyba únikem - leakage, použití váhových oken Když signál není periodický, energie "unikne" do mnoha spektrálních čar blízkých skutečné frekvenci a spektrum je rozprostřeno přes několik čar: periodický signál
neperiodický signál
a(t)
b(t)
čas
čas
T
obdélníkové okno (žádné vážení)
T
A(f)
B(f)
frekvence
3. přednáška
Základy experimentální mechaniky
frekvence 11/32
Chyba únikem - leakage, použití váhových oken periodický signál
neperiodický signál
a(t)
b(t)
čas
čas
T
T
Hanningovo okno
2πt 1 − cos T
A(f)
B(f)
frekvence
3. přednáška
Základy experimentální mechaniky
frekvence 12/32
Chyba únikem - leakage Vztah mezi časovým omezením signálu a chybou únikem ve spektru: spojitý signál a(t)
A(f) čas ×
· W(f)
w(t)
*
*
frekvence
= =
frekvence
čas =
A(f)*W(f)
a(t)⋅w(t) čas
frekvence digitální signál - naměřená data a(t)
DFT
A(f)
čas
∆f = časové omezení
3. přednáška
únik
Základy experimentální mechaniky
1 T
frekvence
13/32
Chyba únikem - leakage Minimalizace chyby únikem: •
Změnou délky trvání měřeného vzorku tak, aby vyhověla základní periodicitě signálu, např. změnou doby měření T ⇒ úplné odstranění chyby, málokdy realizovatelné
•
Zvětšení délky trvání doby měření T, takže frekvenční rozlišení je jemnější ⇒ vliv chyby se zmenší
•
Uzavření signálu do oken - "okenní transformace"
3. přednáška
Základy experimentální mechaniky
14/32
Průměrování Dosavadní poznámky k DFT se týkaly deterministických dat. Vibrační data jsou náhodné signály
⇒
– nestačí vypočíst Fourierovu transformaci - ta ani pro náhodný proces neexistuje – potřebujeme odhady spektrálních hustot a korelačních funkcí – (tyto jsou vypočteny z Fourierovy transformace) – je nutné provést proces průměrování, který zahrne několik jednotlivých časových záznamů – počet požadovaných průměrů ovlivňuje: » požadovaná statistická spolehlivost » míra šumu v signálu
3. přednáška
Základy experimentální mechaniky
15/32
Typy průměrování
bez překrytí
s překrytím
lineární exponenciální s držením špičky
3. přednáška
doba zpracování
Základy experimentální mechaniky
16/32
Typy průměrování
3. přednáška
Základy experimentální mechaniky
17/32
…a frekvenční analýza ještě jednou jinak
K frekvenčním složkám budeme přistupovat jako k rotujícím vektorům – fázorům Na obrázku je vektor F = a + ib zobrazený v komplexní rovině
Re
F = a + ib = |F| eiϕ
F* = a - ib = |F| e-iϕ
a
F = a2 + b2
b ϕ = arctg a
b
|F| ϕ
a = |F| cos ϕ b = |F| sin ϕ
F = F (cos ϕ + i sin ϕ)
-ϕ
e iϕ = cos ϕ + i sin ϕ
Im -ib
F = F ⋅ e iϕ i2b = -b
3. přednáška
F1 ⋅ F2 = F1 ⋅ e iϕ1 ⋅ F2 ⋅ e iϕ2 = F1 ⋅ F2 ⋅ e i(ϕ1 + ϕ2 )
Základy experimentální mechaniky
18/32
Základní předpoklady
3. přednáška
Základem frekvenční analýzy je Fourierova transformace Předpokládá, že signál se skládá z (∞) mnoha (ko)sinusových složek různých frekvencí Každá složka má svou frekvenci, amplitudu a počáteční fázi
Základy experimentální mechaniky
19/32
Zobrazení typické složky
A cos(2πft + φ)
A cos θ = 3. přednáška
(
A iθ e + e −iθ 2
)
Jako součet dvou vektorů s amplitudou A/2 rotujících v opačném smyslu
kde
θ = 2πft + φ
Základy experimentální mechaniky
20/32
Rozklad periodické funkce do Fourierovy řady g(t) je periodická funkce, tzn. g(t) = g(t+nT) Může být vyjádřena jako součet sinusových složek (rotujících vektorů) na frekvencích k·f1, kde f1 = 1/T k … celé číslo včetně nuly a záporných čísel k-tou složku dostaneme z integrálu: T/2
1 − i 2 πf k t ( ) G (f k ) = g t e dt ∫ T −T / 2
3. přednáška
kde fk = k ·f1 ,tj. k-tá harmonická f1
Základy experimentální mechaniky
21/32
Rozklad periodické funkce do Fourierovy řady
Re
Re
φk Im
Im
To znamená, že když signál g(t) obsahuje složku, která rotuje s frekvencí fk, tak násobení jednotkovým vektorem e − i 2 πfk t
(který rotuje s frekvencí –fk) anuluje
rotaci této složky signálu a její integrací v čase dostaneme konečnou hodnotu. Všechny složky na jiných frekvencích stále rotují i po násobení
e −i2 πfk t
a proto jejich integrál za časovou periodu je nulový.
3. přednáška
Základy experimentální mechaniky
22/32
Rozklad periodické funkce do Fourierovy řady T/2
Pomocí integrálu G (f k ) = 1 g(t ) e −i 2 πf k t dt T −T∫/ 2
složky rotující na všech frekvencích fk. Tím taky fázové úhly každé z nich “zamrznou” v poloze, ve které byly v nulovém čase (kdy e −i2 πfk t = 1)
extrahujeme ze signálu g(t)
Skutečnou polohu každého vektoru v libovolném čase lze získat násobením jeho počáteční hodnoty G(fk) opačně rotujícím jednotkovým vektorem i 2 πfk t
e
Celkový signál g(t) bude vektorový součet všech těchto vektorů v jejich okamžitých polohách, tj.:
g(t ) =
3. přednáška
∞
∑ G(f )e
k = −∞
i2 πfk t
k
Základy experimentální mechaniky
23/32
3D zobrazení spektra
Řada komplexních hodnot G(fk) je nazývána spektrum složek signálu g(t) Protože každá z nich má amplitudu a fázi (nebo reálnou a imaginární složku), je k úplnému zobrazení třeba 3D
komplexní spektrum 3. přednáška
Základy experimentální mechaniky
24/32
Vlastnosti spektra
Signál g(t), který je v čase periodický, má diskrétní spektrum, jehož složky mají frekvence, které jsou vždy celočíselnými násobky základní frekvence f1. Vektor se základní frekvencí f1 se otočí o 360º jednou za periodu T. Vektor se základní frekvencí fk se otočí o 360º k-krát za periodu T. ⇒ Po uplynutí periody T se všechny vektory vrátí do své výchozí pozice. Funkce g(t) je reálná funkce, protože každé složce s frekvencí fk odpovídá složka s frekvencí –fk, která má stejnou amplitudu, ale opačnou fázi (tedy stejnou reálnou složku a opačnou imaginární složku). Imaginární složky na všech frekvencích se vyruší a výsledek je vždy reálný. Spektrum reálné funkce je sudá funkce: G(fk) = G*(-fk) Stejnosměrná složka je vždy reálná – fázový úhel je 0 nebo ± π
3. přednáška
Základy experimentální mechaniky
25/32
Výkonové spektrum Okamžitý výkon časového signálu g(t) je roven [g(t)]2 a průměrný výkon za jednu periodu je dán integrací okamžité hodnoty signálu v průběhu periody: T
Pprůměrný
1 2 = ∫ {g(t )} dt T0
Pro typickou složku A k cos(2πfk t + φ k ) to vede na: T
Pprůměrný
T
1 A k2 1 1 A k2 2 2 = ∫ A k cos (2πfk t + φk )dt = − cos [2(2πfk t + φk )] dt = T0 T ∫0 2 2 2
protože integrál sinové části frekvence 2fk za periodu T je nulový. To je známý výsledek pro střední hodnotu sinusovky s amplitudou Ak a vede na hodnotu A efektivní hodnota (RMS ) = k 2
3. přednáška
Základy experimentální mechaniky
26/32
Výkonové spektrum
Výkon na každé frekvenci je dán přímo druhou mocninou amplitudy složky Fourierova spektra:
Amplituda složky G(fk) je Ak/2, kde Ak je amplituda k-té sinusovky, takže její druhá mocnina je Ak2/4.
Jelikož amplitudové spektrum je sudá funkce, platí to podobně pro zápornou frekvenci a celkový výkon na frekvenci fk bude dán součtem těchto složek, tedy Ak2/2, což je stejný výsledek, jaký jsme dostali v časové oblasti.
Celkový výkon může být získán buď integrací časového signálu nebo sečtením druhých mocnin amplitud všech frekvenčních složek (Parsevalův teorém).
Protože ve výkonovém spektru už není obsažena informace o fázi, není možné z něj získat zpět původní časový signál.
3. přednáška
Základy experimentální mechaniky
27/32
Fourierova transformace
Předešlé výsledky se týkaly periodických signálů. To je možné zobecnit na případ, kdy T → ∞ V tomto případě se frekvenční rozlišení 1/T mezi harmonickými blíží nule a G(f) se stává spojitou funkcí frekvence. G(f ) =
∞
∫ g(t )e
−i 2 πft
dt
−∞
∞
g(t ) = ∫ G(f )e i2 πft dt
… přímá Fourierova transformace
… zpětná Fourierova transformace
transformační pár
−∞
Funkce je spojitá v časové i frekvenční oblasti.
3. přednáška
Základy experimentální mechaniky
28/32
Vzorkovaná časová funkce
Funkce reprezentovaná časovými vzorky v ekvidistantních časových krocích (diskrétní v časové oblasti, spojitá ve frekvenční oblasti) Je to opačný případ k Fourierově řadě (spojitá v časové oblasti, diskrétní ve frekvenční oblasti). V důsledku symetrie Fourierova transformačního páru je spektrum periodické s periodou rovnou vzorkovací frekvenci fs fs = 1/∆t Transformační pár pro funkci reprezentovanou časovými vzorky má tvar:
∞
G(f ) =
∑
g(t n )e −i2 πft n
n = −∞
g(t n ) =
3. přednáška
1 fs
fs / 2
∫
G(f )ei2πft n dt
, kde tn = n· ∆t
… n-tý časový vzorek
− fs / 2
Základy experimentální mechaniky
29/32
Diskrétní Fourierova transformace
Funkce je vzorkovaná v časové i frekvenční oblasti ⇒ signál i spektrum jsou implicitně periodické
Transformační pár má tvar: −i 1 N−1 G(k ) = ∑ g(n) e N n=0
N−1
g(n) = ∑ G(k ) e
i
2 πkn N
2 πkn N
… přímá DFT
… zpětná DFT
k =0
3. přednáška
Základy experimentální mechaniky
30/32
Diskrétní Fourierova transformace 1 Agn N Gk … vektor reprezentující N komplexních frekvenčních složek
Vztah pro přímou DFT lze zjednodušeně napsat: kde
Gk =
g n … vektor reprezentující N časových vzorků 1/N … jednoduché měřítko A … čtvercová matice jednotkových vektorů e ↑ G0 G1 ↑ G ↑ 2 G3 1 ↑ = G4 8 ↑ G5 ↑ ↑ G6 G ↑ 7
↑
↑
↑ ↑ ↑
↑
→ ↓ → ↓ ← ↑ → ← ↓ ↓ ↑ ↓ ↑ ↓
← ↓ → ↑
→ ↓ ← ← ↓ → ↑ ← ↓ ← ↓ →
3. přednáška
↑ g0 g1 ← g2 ⋅ g3 ↓ g4 g5 → g6 g7
řady matice …
−i
2 πkn N
frekvence k = 0, 1, 2 …7
sloupce matice … časové okamžiky n = 0, 1, 2 …7
Základy experimentální mechaniky
31/32
Ekvivalence otáčení v kladném a záporném smyslu pro diskrétní funkce Re -1/8 otáčky
Im
+7/8 otáčky
Fourier ova transformace
3. přednáška
Základy experimentální mechaniky
32/32
Děkuji za pozornost !