Primena diskretne Furijeove transformacije

Digitalna obrada signala, Vladimir Petrović

U ovom notebook-u ćemo primeniti DFT u izračunavanju konvolucije dva signala i u dodatnim primerima frekvencijske analize signala.

Brzo izračunavanje konvolucije

Periodična konvolucija periodičnih signala

Ranije smo videli da se konvolucija aperiodičnih signala može izračunati odbirak po odbirak korišćenjem sume:

$$ y[n] = (x*h)[n] = \sum_{m=-\infty}^{+\infty} x[m]h[n-m]. $$

Međutim, ovakva formulacija je neupotrebljiva za periodične signale jer je njihova dužina beskonačna, pa bi konvoluciona suma mogla da divergira. Zbog toga se za periodične signale definiše periodična konvolucija o kojoj će biti reči u ovom odeljku.

Kod periodičnih signala $\overset{\sim}{x}[n]$ smo mogli lako definisati razvoj u Furijeov red:

$$ X[k] = \frac{1}{N} \sum_{n=n_0}^{n_0 + N - 1} \overset{\sim}{x}[n]e^{-j k \Omega_0 n}. $$

Neka je $X_1[k]$ Furijeov red signala $\overset{\sim}{x_1}[n]$ i $X_2[k]$ Furijeov red signala $\overset{\sim}{x_2}[n]$ i neka su signali $\overset{\sim}{x_1}[n]$ i $\overset{\sim}{x_2}[n]$ periodični sa istim osnovnim periodom. Neka je $Y[k] = X_1[k]X_2[k]$ Furijeov red signala $\overset{\sim}{y}[n]$. Tada je signal $\overset{\sim}{y}[n]$:

\begin{align*} \overset{\sim}{y}[n] &= \sum_{k=\langle N \rangle} X_1[k]X_2[k] e^{j k \Omega_0 n} \\ &= \sum_{k=\langle N \rangle} \left( \frac{1}{N} \sum_{m=n_0}^{n_0 + N - 1} \overset{\sim}{x_1}[m]e^{-j k \Omega_0 m} \right) X_2[k] e^{j k \Omega_0 n} \\ &= \frac{1}{N} \sum_{m=n_0}^{n_0 + N - 1} \overset{\sim}{x_1}[m] \sum_{k=\langle N \rangle} X_2[k] e^{j k \Omega_0 (n-m)} \\ &= \frac{1}{N} \sum_{m=n_0}^{n_0 + N - 1} \overset{\sim}{x_1}[m] \overset{\sim}{x_2}[n-m] \end{align*}

Izraz $\sum_{m=n_0}^{n_0 + N - 1} \overset{\sim}{x_1}[m] \overset{\sim}{x_2}[n-m]$ nazivamo periodičnom konvolucijom. Iz prethodne analize zaključujemo da periodičnoj konvoluciji dva signala u vremenskom domenu odgovara proizvod koeficijenata Furijeovog reda pomnožen periodom $N$.

Cirkularna konvolucija

Ranije smo videli da su koeficijenti diskretne Furijeove transformacije signala $x[n]$ isti kao i koeficijenti Furijeovog reda periodičnog produženja signala $x[n]$, $x_p[n]$, samo pomnoženi periodom $N$. Ovo navodi na zaključak i da će proizvod DFT koeficijenata dva signala dati nešto slično periodičnoj konvoluciji periodičnih produženja ta dva signala u vremenskom domenu. Nađimo inverznu diskretnu Furijeovu transformaciju proizvoda koeficijenata DFT-a dve sekvence:

\begin{align*} y[n] &= \frac{1}{N} \sum_{k=0}^{N-1} X_1[k]X_2[k] e^{j k \frac{2\pi}{N} n} = \frac{1}{N} \sum_{k=0}^{N-1} \left(\sum_{m=0}^{N-1} x_1[m]e^{-j k \frac{2\pi}{N} m} \right)X_2[k] e^{j k \frac{2\pi}{N} n} \\ &= \frac{1}{N} \sum_{m=0}^{N-1} x_1[m] \sum_{k=0}^{N-1} X_2[k] e^{j k \frac{2\pi}{N} (n-m)} = \frac{1}{N} \sum_{m=0}^{N-1} x_1[m] \sum_{k=0}^{N-1} NX_{p2}[k] e^{j k \frac{2\pi}{N} (n-m)} \\ &= \sum_{m=0}^{N-1} x_1[m] \sum_{k=0}^{N-1} X_{p2}[k] e^{j k \frac{2\pi}{N} (n-m)} = \sum_{m=0}^{N-1} x_1[m] x_{p2}[n-m] = \sum_{m=0}^{N-1} x_1[m] x_2\langle n-m \rangle_N\\ \end{align*}

gde su sa $X_{p2}[k]$ označeni koeficijenti Furijeovog reda periodičnog produženja signala $x_2[n]$. Sa $x\langle n-m \rangle_N$ označen cirkularni pomeraj sekvence konačne dužine $x[n]$. Na jednom periodu, ovaj pomeraj je ekvivalentan linearnom pomeraju periodičnog produženja te sekvence. Izraz $\sum_{m=0}^{N - 1} x_1[m] x_2\langle n-m \rangle_N = x_1[n] \otimes x_2[n]$ nazivamo cirkularnom konvolucijom dva signala.

Linearna konvolucija uz pomoć DFT-a

S obzirom na to da je diskretnu Furijeovu transformaciju moguće izračunati izuzetno brzo, postavlja se pitanje da li je moguće iskoristiti je za izračunavanje linearne konvolucije. Poznato je da će konvolucija dva signala u vremenskom domenu, u frekvencijskom domenu (Furijeovoj transformaciji) dati proizvod Furijeovih transformacija pojedinačnih signala:

$$ \mathcal{F} \left\lbrace x[n] * h[n] \right\rbrace (e^{j \Omega}) = X(e^{j \Omega})H(e^{j \Omega}). $$

Posmatrajući ovaj izraz lako je zaključiti da je konvolucija dva signala $$ x[n] * h[n] = \mathcal{F}^{-1} \left\lbrace X(e^{j \Omega})H(e^{j \Omega}) \right\rbrace. $$

Međutim, kao što smo videli, ovo ne važi za diskretnu Furijeovu transformaciju u opštem slučaju. Rezultat inverzne DFT proizvoda koeficijenata DFT-a dva signala je cirkularna konvolucija. DFT se može iskoristiti za izračunavanje linearne konvolucije samo onda ako je linearna konvolucija jednaka cirkularnoj. $$ \sum_{m=0}^{N - 1} x_1[m] x_2\langle n-m \rangle_N = x_1[n] \otimes x_2[n] \overset{?}{=} x_1[n] * x_2[n] = \sum_{m=0}^{N - 1} x_1[m] x_2[n-m] $$

Posmatrajući izraze za linearnu i cirkularnu konvoluciju možemo zaključiti da su oni jednaki samo onda kada proizvod dva signala od kojih je jedan cirkularno pomeren daje isti rezultat kao proizvod dva signala gde je jedan linearno pomeren. Ovo se može postići jedino dopunjavanjem nulama. Na kraju, postavlja se pitanje koliko je nula potrebno dopuniti. S obzirom na to da je dužina konvolucije $N_y = N_x + N_h - 1$, gde su $N_x$ i $N_h$ dužine signala koji se konvoluiraju, onda je jasno da broj koeficijenata DFT-a iz kojih se inverznom transformacijom dobija signal dužine $N_y$ mora biti baš $N_y$. Zbog toga se ulazne sekvence dopunjavaju nulama do dužine $N_y$.

Frekvencijska analiza nestacionarnih signala

Do sada smo frekvencijsku analizu radili nad signalima konačnog trajanja koji nisu menjali frekvencijske komponente u vremenu. Učitajmo jedan signal kod koga ovo nije slučaj i nacrtajmo spektar ovog signala.

Kao što vidimo, iz spektra se može malo šta zaključiti. Sve frekvencijske komponente su pomešane i nemoguće je razlučiti kada su pojedine frekvencijske komponente nastale i nestale. Zbog toga se u ovakvim primerima signal mora izdeliti na signale manjeg trajanja, a zatim ispitati frekvencijski sadržaj pojedinačnih delova. Na ovaj način se može dobiti onoliko vektora DFT-a koliko je delova signala. Ti vektori se mogu složiti u matricu, a onda se ta matrica može prikazati kao slika. Tako dobijena slika se naziva spektrogram.

Spektrogram se lako može dobiti korišćenjem ugrađene funkcije iz paketa Scipy: signal.spectrogram(). Ova funkcija vraća podatke o vremenskim trenucima, učestanostima i vrednostima spektra u tim vremenskim trenucima i na tim učestanostima. Da bi se spektrogram prikazao, najjednostavnije je koristiti funkciju pcolormesh(). Međutim ova funkcija može da bude izuzetno spora ako se koristi za crtanje velikih matrica. Kao obilazno rešenje se može koristiti funkcija matshow() ili imshow(). Ove funkcije su malo komplikovanije za korišćenje kada se koriste za crtanje spektrograma. Zbog toga se posle crtanja spektra mora podesiti odnos širine i visine slike. U narednim primerima se za to koristi funkcija forceAspect() definisana u narednoj ćeliji. Prilikom pozivanja funkcija matshow() ili imshow() moraju definisati ključne tačke za obeležavanje osa. To se postiže podešavanjem parametra extent. Takođe, podrazumevano je da je koordinatni početak u gornjem levom uglu, s obzirom na to da se ove funkcije najčešće koriste za prikaz slika. Ako želimo da koordinatni početak bude u donjem levom uglu, moramo podesiti i parametar origin na vrednost 'lower'.

Nacrtajmo sada spektrogram signala iz prethodnog primera.

Primećujemo da je sada moguće uočiti da tonovi muzičkog signala postaju sve viši i viši do sredine trajanja signala, a zatim postaju sve niži i niži. Pokažimo još jedan primer, ovoga puta govornog signala:

Više o različitim mapama boja možete naći ovde.