4. KRUHOVÁ KONVOLUCE, RYCHLÁ FOURIEROVA TRANSFORMACE (FFT) A SPEKTRÁLNÍ ANALÝZA SIGNÁLŮ
Kruhová (cyklická) konvoluce Rychlá Fourierova transformace Aplikace DFT na analogové signály, frekvenční analýza periodických analogových signálů s využitím DFT
383ZS – P4
1
KRUHOVÁ (CYKLICKÁ) KONVOLUCE Kruhová konvoluce dvou periodických signálů xP(n) a h xP(n) s periodou N je
y P ( n ) = x P ( n ) ⊗ hP ( n ) =
N −1
N −1
i =0
i =0
∑ xP (i) ⋅ hP (n − i) = ∑ xP (n − i) ⋅ hP (i)
Pozn.: Konvoluce signálů diskrétních v čase je
y ( n) = x ( n) * h ( n) =
∞
∞
i = −∞
i = −∞
∑ x (i ) ⋅ h ( n − i ) = ∑ x ( n − i ) ⋅ h ( i )
Pro odlišení od kruhové konvoluce se též nazývá „neperiodická“ (někdy, ale nevhodně, také „lineární“) konvoluce. Pro DFT platí
yP (n) = xP (n) ⊗ hP (n) ↔ X P (k ) ⋅ H P (k ) = YP (k ) Grafický postup nalezeni kruhové konvoluce (KK). (KK je periodická se stejnou periodou jako vstupní signály (délka základního intervalu KK je stejná jako délka vstupních signálů, tj. N).
383ZS – P4
2
Dvě možnosti nalezení průběhu KK v základním intervalu
1) y P ( n) = xP ( n) ⊗ hP ( n) 2) xP ( n) ↔ X P ( k ), hP ( n) ↔ H P ( k ),
DFT
YP ( k ) ↔ X P ( k ) ⋅ H P ( k ) YP ( k ) ↔ y P ( n)
IDFT
Druhý způsob – rychlejší, pokud se použije FFT (tzv. rychlá konvoluce).
UŽITÍ KK PRO NALEZENÍ (NEPERIODICKÉ) KONVOLUCE KK definována pouze pro periodické posloupnosti, lze ji ale užít i pro nalezení (neperiodické konvoluce). Je nutno zařídit, aby oba typy konvoluce měly stejný průběh v základním intervalu, tj. pro body 0 až N -1. Postup závisí na délce vstupních signálů.
383ZS – P4
3
A) KONVOLUCE SIGNÁLŮ KONEČNÉ DÉLKY
x(n ), N1 , h( n), N 2 , x(n ) * h(n), Zvolit: N = N1 + N 2 − 1
N1 + N 2 − 1
Každý ze signálů x(n) a h(n) se doplní nulami na délku N1+N2-1 nebo větší, v praxi většinou na mocninu dvou. Pomocí algoritmu rychlé konvoluce (užití DFT a IDFT) se najde KK a její průběh v intervalu délky N1+N2 -1 je neperiodická konvoluce.
383ZS – P4
4
383ZS – P4
5
B) KONVOLUCE SIGNÁLU KONEČNÉ DÉLKY SE SIGNÁLEM “NEKONEČNĚ” DLOUHÝM “Nekonečně dlouhý signál” x(n) - např.: vzorkovaná řeč, vzorkovaný signál ze snímače. Signál konečné délky h(n) má délku N1 a je to např. impulsní odezva FIR filtru.
y ( n) = x( n) * h ( n) =
∞
∑ x (i ) ⋅ h ( n − i )
i = −∞
Přímý výpočet: je možný, lze ho urychlit pomocí DFT a IDFT, ale první hodnota y(n) se najde až po skončení „ ∞ “ dlouhého signálu x(n). Řešení – segmentace signálu x(n). Rozdělení x(n) na M segmentů, každý o délce L. Konvoluce je součtem M dílčích konvolucí:
y (n ) =
∞
∞
i = −∞
i = −∞
∑ x1 (i) ⋅h(n − i) + ∑ x2 (i) ⋅h(n − i) + L = y1 (n ) + y1 (n ) + L + yM (n )
Postupný výpočet dílčích konvolucí (pomocí kruhové konvoluce a DFT a IDFT). Každá má délku L+N1 -1, takže u každého segmentu je překryv N1 -1 hodnot. N1 je délka signálu h(n). Pro výpočet dílčích rychlých konvolucí stačí jen jednou vypočíst H(k), použije se ve všech výpočtech. Do výsledného součtu se musí zahrnout i překryvy dílčích konvolucí. Jde o metodu s částečným překrýváním výsledků – „overlap-add method“. V MATLABu = funkce fftfilt(h,x,y) 383ZS – P4
6
RYCHLÁ FOURIEROVA TRANSFORMACE (FFT – FAST FOURIER TRANSFORM) Z definice DFT a IDFT – nalezeni každého bodu: N (komplexních) násobení a (N -1) sčítání. Celkem: N2 násobení a N.(N-1) sčítání, tj. počet operací je řádu N2. Pro velká N – výpočet příliš dlouhý. Řešení: modifikace algoritmu pro výpočet DFT a IDFT tak, aby počet operací klesl. Prakticky použitelný algoritmus – z r.1965 (J.W.Cooley, J.W. Tukey, An algorithm for machine calculation of complex Fourier series. Math. Computation, 4 stránky)
Princip: Výpočet řady DFT kratší délky a kombinování jejich výsledků. Délka signálu pro základní algoritmy FFT je N = 2M, M - přirozené. Počet operací je řádu N.M čili N.log2N. +1 V teorii FFT se užívá označení exp(− j 2π / N ) = W N , takže FFT je dána vztahy
X (k ) =
N −1
∑
n=0
383ZS – P4
x ( n) ⋅ WN+ kn
1 x ( n) = N
N −1
− kn X ( k ) ⋅ W ∑ N
n=0
7
NEJBĚŽNĚJŠÍ ALGORITMUS: DIT (DECIMATION IN TIME, FFT S VÝBĚREM VE VSTUPNÍ POSLOUPNOSTI) N – bodová DFT signálu x(n) o délce N = 2M se vypočte kombinováním výsledků dvou N/2 bodových DFT, počítaných ze sudých vzorků signálu (x1(n)) a lichých vzorků signálu (x2(n)). Počet operací klesne z N2 na 2 . (N/2)2 =N2/2, čili o 50 %. Tento postup se opakuje při výpočtu obou (N/2) bodových DFT, ty se počítají pomocí (N/4) bodových DFT atd. Protože je N = 2M , je těchto kroků M. Poslední krok je kombinování výsledků M 2-bodových DFT.
APLIKACE DFT NA ANALOGOVÉ SIGNÁLY, DFT SPEKTRÁLNÍ ANALÝZA PERIODICKÝCH SIGNÁLŮ DFT (a FFT) – definovány pro signál diskrétní v čase. V praxi se používá i pro přibližné zjištění spektra analogových signálů (spojitých v čase i amplitudě). V důsledku konečného počtu vzorků N v časové i frekvenční oblasti, nevhodné vzorkovací frekvence nebo nesplnění podmínky navzorkování celého počtu period analogového signálu vznikají chyby – chyba useknutím, aliasing, leakage. Správné spektrum najdeme aplikací spojité Fourierovy Transformace: H (ω) = FT {h(t )} . 383ZS – P4
8
Princip nalezení spektra signálu h(t) spojitého v čase pomocí DFT: 1) Na signál aplikujeme obdélníkové okno délky TM (tj. pracujeme s úsekem signálu délky TM ): hw (t ) = h(t ) ⋅ w(t ) 2) Signál hw (t ) ovzorkujeme a aplikujeme na něj DFT: H P ( k ) = DFT {hw (n)} Spektrum H (ω) se ale liší od spektra H P (k ) amplitudou i tvarem. Může vzniknout chyba useknutím a aliasing. Protože je fyzikálně nemožné, aby h(t ) i H (ω) byly oba konečné délky, nelze zabránit oběma chybám současně. Projeví se alespoň jedna, v praxi často obě. Navíc se v reálném systému projeví chyby použitého AD převodníku (amplitudová a časová). Podrobnější výklad - viz následující obrázek.
383ZS – P4
9
383ZS – P4
10
UŽITÍ DFT PRO FREKVENČNÍ NALÝZU PERIODICKÉHO SIGNÁLU SPOJITÉHO V ČASE Periodičnost x(n) a X(k), možnost prosakování energie ve spektru (leakage). K tomu nedojde, je-li splněna podmínka koherentního vzorkování N ⋅ T = m ⋅ TSIG , kterou můžeme vyjádřit f SIG f T ⋅N = = SIG = m, m − celé kladné ( počet navzorkova ných period ) TSIG fVZ / N ∆f FFT Pak pro reálný signál x(n) dostaneme nezkreslené spektrum v intervalu f ∈ 〈 0, fVZ / 2〉 , který zobrazují FFT analyzátory. Pokud ale při vzorkování platí f SIG f T ⋅N 1 = = SIG = m + α, α < , α≠0 TSIG fVZ / N ∆f FFT 2 dojde při výpočtu DFT k prosakování energie do okolních frekvenčních binů (Spektrum se „rozmaže“, čili místo jediné čáry v základním intervalu 0 až N -1odpovídající koherentně ovzorkované sinusovce resp. exponenciále s imaginárním exponentem se objeví čar několik. Celková energie signálu se nezmění, ale rozdělí se do širšího frekvenčního pásma).
383ZS – P4
11
Spektrum komplexní exponenciály pro m=3, α=0 a α=0.3) (v obr. je α označeno jako a)
Pokud nelze zajistit koherentní vzorkování, lze zvýšit rozlišovací schopnost ve spektru zmenšením vzdálenosti sousedních spektrálních čar DFT spektra. Tyto čáry leží na tzv. mřížce DFT (DFT grid), čili na frekvencích v násobcích hodnoty ∆f = 1 (Hz ) . Hodnota ∆f je tzv. „frekvenční bin“. Zhuštění DFT mřížky (při NTVZ vyjádření frekvence na vodorovné ose v Hz, nikoliv v bezrozměrné hodnotě k !) se dosáhne zmenšením ∆f , čili zvětšením doby odebírání vzorků TM = N ⋅ TVZ . To lze provést zvětšením délky DFT (čili zvětšením počtu vzorků odebraného signálu) při zachování vzorkovací frekvence nebo zvětšením vzorkovacího intervalu TVZ , čili snížením vzorkovací frekvence při zachování počtu vzorků. 383ZS – P4
12
Nelze-li zajistit koherentní vzorkování, potlačuje se v praxi leakage okénkováním signálu. Použití vhodného okna sníží postranní čáry spektra harmonického signálu, ale změní energii signálu (a tedy např. efektivní hodnotu sinusového napětí měřeného signálu). Nezajímá-li nás jen tvar spektra, ale chceme-li zachovat energii signálu, je nutno spektrum násobit opravným koeficientem závislým na tvaru použitého okna. Dokonalejší metodou je použití interpolované DFT, kde se pomocí interpolace určí přesná frekvence signálu (která neleží na DFT mřížce) a na základě znalosti průběhu spektra okna se určí správná výška spektrální čáry. Spektrum komplexní exponenciály pro m=3, α=0 a α=0.3) při užití okna Hann (v obr. je α označeno jako a)
383ZS – P4
13
Je-li splněna podmínka koherentního vzorkování (je navzorkován celý počet period analogového signálu), k leakage nedojde. Pak je pro spektrální analýzu nejlepší nepoužít žádné okno (tj. použít okno obdélníkové, protože pracujeme s konečnou délkou signálu). Okno Hann a jeho spektrum. Pokud nebyl navzorkován celý počet period, pak při spektrální analýze použijeme okno, čili analyzujeme signál xW ( n) = x ( n) ⋅ w( n). Příkladem často používaného okna je okno Hann (někdy nazývané „hanning“ resp. „Hanning“), jehož spektrum je tvořeno superposicí tří spekter okna obdélníkového. Okno může být položeno souměrně kolem počátku souřadnic (pak je sudé a má sudé a reálné spektrum). Takové okno Hann značíme wHN 1 (n). Okno může ale začínat v počátku souřadnic. Pak je posunutým sudým oknem, zde oknem wHN 1 (n). Označíme je wHN 2 (n). Toto okno se používá v praxi číslicového zpracování signálu (protože fyzikální signály se generují pro kladné časy). Jeho amplitudové spektrum je totožné se spektrem sudého okna, ale jeho fázové spektrum je nenulové (okno má lineární fázi).
383ZS – P4
14
Pro okna Hann platí
wHN 1 (n) =
N N 1 2π n , − ≤ n < 1 + cos N 2 2 2
wHN 2 ( n) =
1 2π n , 0 ≤ n < N 1 − cos N 2
FTD spektrum okna wHN1 (n) je
1 1 WHN1 (e jθ ) = WO (e jθ ) + WO (e 2 4
2π j θ − N
1 ) + WO (e 4
2π j θ + N
) jθ
FTD spektrum sudého okna Hann je tedy součtem tří spekter obdélníkového okna WO (e ) . jθ
Spektrum WO (e ) je sudé a násobené koeficientem 0.5, zbývající dvě spektra jsou násobena 2π koeficienty 0,25 a posunuta proti prostřednímu spektru o hodnoty ∆θ = ± N ., které při fVZ vyjádření frekvence v Hz odpovídají hodnotě frekvenčního binu ∆f = ± N . Protože DFT
2π
spektrum odpovídá vzorkovanému spektru FTD v hodnotách ∆θ = ± k ⋅ N , je DFT spektrum sudého okna Hann tvořeno jen třemi vzorky a to (N/4) pro k=±1 a N/2 pro k=0. (V obrázku je θ označeno jako q a π jako p.)
383ZS – P4
15
FTD spektrum okna Hann:
DFT spektrum okna Hann je FTD spektrum okna Hann vzorkované v bodech k.2π/N, tj. 3 spektrální čáry, délky 1/4, 1/2 a 1/4 pro wHN1 a –1/4, 1/2 a –1/4 pro wH2.
383ZS – P4
16
