Digitalna obrada signala, Vladimir Petrović
Ovaj notebook će predstaviti osnovni alat za frakvencijsku analizu signala u digitalnoj obradi signala - diskretnu Furijeovu transformaciju. Najpre ćemo definisati diskretnu Furijeovu transformaciju, ispitaćemo neke njene osobine, a zatim i primeniti je na signalima iz realnog sveta.
USE_WIDGETS = False
def importEssentialLibs(USE_WIDGETS):
import numpy as np
if USE_WIDGETS:
%matplotlib widget
else:
%matplotlib inline
import matplotlib as mpl
mpl.rc('text', usetex = True)
mpl.rc('font', family = 'serif', size = 18)
import matplotlib.pyplot as plt
return np, mpl, plt
Kao što je ranije rečeno, svaki diskretni signal se može predstaviti kao linearna kombinacija nekih bazisnih signala. Najjednostavniji bazisni signali su vremenski pomereni jedinični impulsi, pa smo videli ranije da se svaki signal može predstaviti kao suma skaliranih jediničnih impulsa (linearna kombinacija):
$$ x[n] = \sum_{m=-\infty}^{+\infty} x[m]\delta[n-m]. $$Međutim, bazisni signali ne moraju nužno da budu jedinični impulsi. Dodatno, ne mora ih imati prebrojivo mnogo, kao u prethodnom primeru. Sve dok je moguće signal zapisati kao linearnu kombinaciju nekih drugih signala, ti drugi signali se mogu smatrati bazisnim signalima. Primetite da koristimo izraz ,,bazisni signal" kao ekvivalentan izrazu ,,bazisni vektor". Dimenzije vektora su beskonačne. Ako je signal periodičan, onda se on može predstaviti linearnom kombinacijom konačnog broja kompleksnih sinusoida čije su relativne učestanosti celobrojni umnožak osnovne relativne učestanosti periodičnog signala. Broj kompleksnih sinusoida odgovara broju odbiraka u jednoj periodi signala, kao u sledećem izrazu
$$ x[n] = \sum_{k=\langle N \rangle} X[k]e^{j k \Omega_0 n}, $$gde je $\Omega_0 = 2\pi/N$, a $k=\langle N \rangle$ znači da se indeks $k$ kreće unutar bilo kog niza od $N$ uzastopnih celih brojeva. Koeficijenti uz bazisne signale (vektore) $e^{j k \Omega_0 n}$ su projekcije vektora $x[n]$ na te bazisne vektore:
$$ X[k] = \frac{1}{N} \sum_{n=n_0}^{n_0 + N - 1} x[n]e^{-j k \Omega_0 n}. $$Prethodni izraz nam je poznat kao razvoj u Furijeov red (Discrete Time Fourier Series - DTFS). S obzirom na to da je $e^{j k \Omega_0 n} = e^{j k n 2\pi/N} \cdot 1 = e^{j k n 2\pi/N}e^{j k 2\pi n} = e^{j (k+N) n 2\pi/N}$, jasno je da je koeficijent $X[k]$ isti kao i koeficijent $X[k+N]$ jer oba množe isti bazisni signal $e^{j k \Omega_0 n} = e^{j (k+N) n 2\pi/N}$ u linearnoj kombinaciji kompleksnih sinusoida u predstavi signala $x[n]$. Dakle, Furijeov red diskretnih signala je periodičan sa periodom $N$.
Ako signal nije periodičan, može se formalno reći da mu perioda $N$ teži beskonačnosti, pa u tom slučaju relativna kružna učestanost bazisnih signala $\Omega_k = k \Omega_0 = k \cdot 2\pi/N$ postaje kontinualna promenljiva $\Omega$, pa se signal ne može predstaviti linearnom kombinacijom konačnog broja kompleksnih sinusoida (bazisnih signala), ali se može predstaviti linearnom kombinacijom beskonačnog broja bazisnih signala:
$$ x[n] = \frac{1}{2\pi} \int_{2\pi} X(j\Omega)e^{-j \Omega n} d\Omega. $$Koeficijenti uz bazisne vektore su, kao i ranije, projekcije vektora $x[n]$ na te bazisne vektore:
$$ X(j\Omega) = \sum_{n=-\infty}^{+\infty} x[n]e^{-j \Omega n}. $$Prethodni izraz nam je poznat kao Furijeova transformacija diskretnog signala (Discrete Time Fourier Transform - DTFT).
Koeficijenti Furijeovog reda i koeficijenti Furijeove transformacije su periodični sa periodom $N$ i $2\pi$, respektivno. S obzirom na to da indeks $k$ odgovara relativnoj kružnoj učestanosti $k \cdot 2\pi/N$, jasno je da je i Furijeov red periodičan sa periodom $2\pi$ ako indekse preslikamo u relativne kružne učestanosti.
Primetimo da su do sada analizirani signali bili beskonačnog trajanja, što je neprihvatljivo u slučaju da želimo da računar uradi frekvencijsku analizu nepoznatog signala, jer računar ne može da radi sa beskonačnim brojem odbiraka odjedanput. Dakle, u digitalnoj obradi signala, signal beskonačnog ili veoma velikog trajanja će na ulazu u digitalni sistem za frekvencijsku analizu morati da se podeli na signale konačnog trajanja (signale koji imaju konačan broj nenultih odbiraka). Zatim će svaki od pojedinačnih delova signala biti analiziran nezavisno.
Odsecanjem delova signala koji je možda i bio periodičan, napravili smo aperiodične signale, pa bi stoga intuitivno bilo da koristimo Furijeovu transformaciju za frevencijsku analizu. Međutim, računar ne može da sačuva beskonačan broj koeficijenata Furijeove transformacije u konačne memorijske resurse, pa korišćenje Furijeove transformacije ne dolazi u obzir. Ipak, možemo izračunati Furijeovu transformaciju signala samo na nekim relativnim učestanostima i taj ograničen broj odbiraka Furijeove transformacije smestiti u memoriju. Zatim ostaje da se nadamo da ćemo iz tih odbiraka moći da izvučemo neke zaključke tokom analize signala.
Diskretna Furijeova transformacija (Discrete Fourier Transform - DFT) nije ništa drugo do Furijeova transformacija izračunata u onoliko ekvidistantnih tačaka koliko je i odbiraka diskretnog signala.
$$ X[k] = X(j\Omega_k) = X(j k \frac{2\pi}{N}) = \sum_{n=0}^{N-1} x[n]e^{-j k \frac{2\pi}{N} n},\, k \in [0, N-1]. $$Ponovimo izraz sa Furijeov red periodičnog diskretnog signala $x_p[n]$:
$$ X_p[k] = \frac{1}{N} \sum_{n=n_0}^{n_0 + N} x_p[n]e^{-j k \Omega_0 n} = \frac{1}{N} \sum_{n=n_0}^{n_0 + N} x_p[n]e^{-j k \frac{2\pi}{N} n}. $$Ako usvojimo da je $n_0 = 0$, dobijamo potpuno isti izraz kao u slučaju da računamo DFT jedne periode signala $x_p[n]$, nedostaje samo faktor skaliranja. Obrnuto gledano, diskretna Furijeova transformacija vremenski ograničenog signala $x[n]$ je, praktično, Furijeov red signala koji predstavlja beskonačno periodično produženje tog signala $x[n]$.
Iz koeficijenata diskretne Furijeove transformacije se može dobiti originalni signal (ponovo linearna kombinacija bazisnih vektora) na sledeći način:
$$ x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k]e^{j k \frac{2\pi}{N} n}. $$Ovaj postupak nazivamo inverznom diskretnom Furijeovom transformacijom (Inverse Discrete Fourier Transform - IDFT).
Da zaključimo, svaki diskretan signal konačnog trajanja koji možemo predstaviti vektorom konačne dužine može da se predstavi linearnom kombinacijom bazisnih vektora $e^{j k \frac{2\pi}{N} n}$, za $k = 0, 1, 2, ..., N-1$. Prikažimo ove bazisne vektore u Pajton programu:
np, mpl, plt = importEssentialLibs(USE_WIDGETS)
import matplotlib.colors as colors
from matplotlib.patches import FancyArrow
colorsList = list(colors.TABLEAU_COLORS.values())
def basisVectorsDft(N = 32):
# vektor a (dimenzija Nx1) - niz brojeva on 0 do N-1
a = np.expand_dims(np.arange(N), 0)
# matrica baziskinh vektora, broadcasting pomaze
U = np.exp(2j * (np.pi / N) * (a.T * a))
return U
N = 16
v = np.ones((N,1))
U = basisVectorsDft(N)
fig, ax0 = plt.subplots(figsize=(4,4))
ax0.set_aspect(1)
a = 2*np.pi/N*np.arange(N)
for j, i in enumerate(a):
ax0.add_patch( FancyArrow(0, 0, np.cos(i), np.sin(i), width=0.02,
length_includes_head=True, color=colorsList[j % 8]))
ax0.text(1, 0.1, '0')
ax0.text(0.1, 1, r'$\frac{\pi}{2}$', fontsize = 22)
ax0.text(-1, 0.1, r'$\pi$')
ax0.text(0.1, -1.2, r'$\frac{3\pi}{2}$', fontsize = 22)
ax0.set_xlim(-1.3, 1.3);
ax0.set_ylim(-1.3, 1.3);
ax0.set_xticks([])
ax0.set_yticks([])
ax0.set_title('Relativna kružna učestanost $\Omega$')
mpl.rc('font', family = 'serif', size = 12)
fig, ax = plt.subplots(8, 2 ,figsize=(15,15))
plt.subplots_adjust(hspace = 1)
for j in range(2):
for i in range(8):
ax[i][j].plot(U.real[:,i+j*8], '-o', color=colorsList[i])
ax[i][j].plot(U.imag[:,i+j*8], '--o', color=colorsList[i], alpha=0.3)
ax[i][j].set_ylim(-1.2, 1.2);
ax[i][j].title.set_text(r'$\Omega_{%d}=%d\times\frac{2\pi}{16}$'%(i+j*8,i+j*8))
ax[i][j].set_xlim(0, 15);
if i == 7:
ax[i][j].set_xlabel(r'$n$', fontsize = 18)
ax[i][j].set_xticks(np.arange(16))
Primer 1 Napišimo funkciju za DFT po definiciji i nacrtajmo rezultate izračunavanja za više test signala.
np, mpl, plt = importEssentialLibs(USE_WIDGETS)
import matplotlib.colors as colors
def dftDef(x):
N = len(x)
X = np.zeros(N, dtype = complex)
n = np.arange(N)
WN = np.exp(-2j*np.pi/N*n)
for k in range(N):
X[k] = np.sum(x*(WN**k))
return X
def plotSignalAndDft(x, X):
n = np.arange(len(x))
fig, ax = plt.subplots(2,1)
plt.subplots_adjust(hspace=0.5)
ax[0].plot(n, x.real, 'o-', label = 'real', color = 'tab:blue')
ax[0].plot(n, x.imag, 'o--', label = 'imag', color = 'tab:red')
ax[0].set_ylabel('$x[n]$')
ax[0].set_ylim(np.min([x.real, x.imag])-0.1, np.max([x.real, x.imag])+0.1)
ax[0].set_xlabel('$n$')
ax[0].set_xlim(0, N*1.3)
ax[0].legend(loc="upper right");
k = n
ax[1].stem(k, X.real, label = 'real')
markerline, stemlines, baseline = ax[1].stem(k, X.imag, label = 'imag')
plt.setp(markerline, markerfacecolor = 'tab:red', markeredgecolor = 'tab:red')
plt.setp(stemlines, color = 'tab:red')
ax[1].set_ylabel('$X[k]$')
ax[1].set_ylim(np.min([X.real, X.imag])-1, np.max([X.real, X.imag])+1)
ax[1].set_xlabel('$k$')
ax[1].set_xlim(0, N*1.3)
ax[1].legend(loc="upper right");
N = 32
n = np.arange(N)
# relativna ucestanost
F = 1/8
#pocetna faza
Phi = 0 #np.pi/3#, 2*np.pi/3
x = np.exp(1j*(2*np.pi*F*n + Phi))
X = dftDef(x)
plotSignalAndDft(x, X)
# relativne ucestanosti
F1 = 1/8
F2 = 3/8
#pocetne faze
Phi1 = 0 # np.pi/3, 2*np.pi/3
Phi2 = np.pi/3
x = np.exp(1j*(2*np.pi*F1*n + Phi1)) + np.exp(1j*(2*np.pi*F2*n + Phi2))
X = dftDef(x)
plotSignalAndDft(x, X)
x = np.cos(2*np.pi*F1*n + Phi1) + np.cos(2*np.pi*F2*n + Phi2)
X = dftDef(x)
plotSignalAndDft(x, X)
Vežba: Napišite funkciju za inverznu diskretnu Furijeovu transformaciju i proverite njenu funkcionalnost.
Često se izraz $e^{-j \frac{2\pi}{N}}$ predstavlja sa $W_N$ radi lakšeg zapisa. Na taj način, izraz za diskretnu Furijeovu transformaciju postaje:
$$ X[k] = \sum_{n=0}^{N-1} x[n]e^{-j k \frac{2\pi}{N} n} = \sum_{n=0}^{N-1} x[n]W^{kn}. $$Ako sve pojedinačne vektore $W^{kn}$ (sa $N$ vrsta i jednom kolonom) poređamo u matricu dobijamo sledeću matricu:
$$ \mathbf{W}_N = \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & W^{1 \cdot 1} & W^{1 \cdot 2} & \cdots & W^{1 \cdot (N-1)} \\ 1 & W^{2 \cdot 1} & W^{2 \cdot 2} & \cdots & W^{2 \cdot (N-1)} \\ 1 & W^{3 \cdot 1} & W^{3 \cdot 2} & \cdots & W^{3 \cdot (N-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & W^{(N-1) \cdot 1} & W^{(N-1) \cdot 2} & \cdots & W^{(N-1) \cdot (N-1)} \\ \end{bmatrix} $$Vektor koeficijenata diskretne Furijeove transformacije sada lako možemo dobiti množenjem ove matrice i vektora signala: $$ \mathbf{X} = \mathbf{W}_N \mathbf{x}. $$
Matrica $\mathbf{W}_N$ ima lepu osobinu da je $\mathbf{W}_N^*\mathbf{W}_N = N\mathbf{I}$, pa se inverzna disretna Furijeova transformacija, takođe lako piše u matričnom obliku:
$$ \mathbf{W}_N^*\mathbf{X} = \mathbf{W}_N^*\mathbf{W}_N \mathbf{x} = N \mathbf{x} \implies \mathbf{x} = \frac{1}{N} \mathbf{W}_N^*\mathbf{X}. $$Primer 2 Napišimo funkciju za DFT u matričnoj formi i uporedimo rezultate. Matricu $\mathbf{W}_N$ možemo i izračunati samo jedanput, ako ćemo uvek raditi sa istom dužinom signala $N$. Ovim se značajno ubrzava ponovno izračunavanje DFT-a. Zbog toga ćemo u primeru izračunavanje DFT-a razdvojiti na dve funkcije: 1) izračunavanje matrice $\mathbf{W}_N$ i 2) izračunavanje DFT-a. Napomena: Ovaj primer se nadovezuje na Primer 1, pa je podrazumevano da je pre pokretanja Primera 2 pokrenut Primer 1 u kome su definisane funkcije i uključene biblioteke koje se koriste i u Primeru 2.
def calculateDftsW(N):
# vektor a (dimenzija Nx1) - niz brojeva on 0 do N-1
a = np.expand_dims(np.arange(N), 0)
# matrica W, broadcasting pomaze
W = np.exp(-2j * (np.pi / N) * (a.T * a))
return W
def dftDefMatrixForm(W, x):
X = np.squeeze(np.dot(W, x))
return X
N = 32
W = calculateDftsW(N)
n = np.arange(N)
# relativna ucestanost
F = 1/8
#pocetna faza
Phi = 0 # np.pi/3, 2*np.pi/3
x = np.exp(1j*(2*np.pi*F*n + Phi))
X = dftDefMatrixForm(W, x)
plotSignalAndDft(x, X)
# relativne ucestanosti
F1 = 1/8
F2 = 3/8
#pocetne faze
Phi1 = 0 # np.pi/3, 2*np.pi/3
Phi2 = np.pi/3
x = np.exp(1j*(2*np.pi*F1*n + Phi1)) + np.exp(1j*(2*np.pi*F2*n + Phi2))
X = dftDefMatrixForm(W, x)
plotSignalAndDft(x, X)
x = np.cos(2*np.pi*F1*n + Phi1) + np.cos(2*np.pi*F2*n + Phi2)
X = dftDefMatrixForm(W, x)
plotSignalAndDft(x, X)
Vežba: Napišite funkciju za inverznu diskretnu Furijeovu transformaciju u matričnoj formi i proverite njenu funkcionalnost na prethodnim primerima.
Diskretna Furijeova transformacija se može izračunati značajno brže i efikasnije nego što se to postiže izračunavanjem po definiciji koje smo videli u prethodnim primerima. Postoji čitav set algoritama koji ubrzavaju izračunavanje DFT-a mnogo puta. O njima više na predavanjima. Programski paketi numpy i scipy imaju implementirane ove algoritme i njih ćemo koristiti ubuduće. Ako se za izračunavanje diskretne Furijeove transformacije koriste brzi algoritmi, često se takva diskretna Furijeova transformacija naziva brzom Furijeovom transformacijom (Fast Fourier Transform - FFT).
Funkcije vezane za FFT ćemo koristiti iz paketa scipy iz biblioteke fft.
import scipy.fft as fft
Funkcija za izračunavanje FFT-a se poziva na sledeći način:
X = fft.fft(x)
dok se inverzna diskretna Furijeova transformacija se poziva na sledeći način:
x = fft.ifft(X)
Primer 3: Na signale iz primera 1 i 2 primeniti ugrađene funkcije za brzu Furijeovu transformaciju, uporediti rezultate. Izračunati inverznu Furijeovu transformaciju i uporediti originalni signal i signal dobijen IDFT-om.
import scipy.fft as fft
N = 32
# relativne ucestanosti
F1 = 1/8
F2 = 3/8
#pocetne faze
Phi1 = 0 # np.pi/3, 2*np.pi/3
Phi2 = np.pi/3
x = np.cos(2*np.pi*F1*n + Phi1) + np.cos(2*np.pi*F2*n + Phi2)
X = fft.fft(x)
plotSignalAndDft(x, X)
xhat = fft.ifft(X)
diff = np.sum(np.real(x - xhat))
print(f'Razlika signala je: %d' % (diff))
Razlika signala je: 0
Kao i kod Furijeove transformacije i kod diskretne Furijeove transformacije je često potrebno da izračunamo amplitudsku i faznu karakteristiku. Tačnije, kod DFT-a ćemo izračunati samo odbirke amplidutdske i fazne karakteristike jer uz pomoć DFT-a dobijamo samo odbirke frekvencijske karakteristike. Koeficijenti su kompleksni brojevi koji se mogu zapisati u Ojlerovom obliku:
$$ X[k] = |X[k]| e^{j \angle(X[k])}, $$gde su $|X[k]|$ odbirci amplitudske, a $\angle(X[k])$ odbirci fazne karakteristike.
np, mpl, plt = importEssentialLibs(USE_WIDGETS)
import scipy.fft as fft
n = np.arange(32)
# relativne ucestanosti
F1 = 1/8
F2 = 3/8
#pocetne faze
Phi1 = 2*np.pi/3 # np.pi/3, 2*np.pi/3
Phi2 = np.pi/3
x = np.cos(2*np.pi*F1*n + Phi1) + np.cos(2*np.pi*F2*n + Phi2)
X = fft.fft(x)
fig = plt.figure(figsize = [6, 3])
X = fft.fft(x);
plt.stem(np.abs(X));
plt.subplots_adjust(bottom=0.25, left=0.15)
plt.xlabel(r'$k$')
plt.ylabel(r'$|X[k]|$');
fig = plt.figure(figsize = [6, 3])
plt.stem(np.angle(X));
plt.subplots_adjust(bottom=0.25, left=0.15)
plt.xlabel(r'$k$')
plt.ylabel(r'$\angle(X[k])$');
Primetimo jednu zbunjujuću stvar. Fazna karakteristika ima vrlo značajne vrednosti za gotovo sve koeficijente $k$. Međutim, ovo nije teoretski očekivano. Zašto je onda alat dao pogrešne rezultate? Reč je o konačnoj tačnosti izračunavanja. Koeficijenti DFT-a izračunati na računaru iako izuzetno male vrednosti, nisu tačno 0. Proverimo ovo:
print(X)
[-1.55431223e-15-0.00000000e+00j 4.14018138e-16+2.48398958e-15j 2.83430252e-15+2.54524974e-16j -1.25130829e-14+1.17522219e-15j -8.00000000e+00+1.38564065e+01j -1.00271283e-14-1.05969226e-14j -6.80594692e-15-9.73035361e-15j 8.12534275e-15-1.05483352e-14j 1.26565425e-14+4.44089210e-16j 5.67555423e-15-3.51142572e-15j 6.74553348e-15+1.31875814e-14j -1.60045167e-14+5.72145760e-15j 8.00000000e+00+1.38564065e+01j -1.26979923e-14+6.44907815e-15j -8.54704881e-15-1.32428553e-14j 1.30469877e-14-3.29998201e-15j 1.11022302e-15-0.00000000e+00j 1.30469877e-14+3.29998201e-15j -8.54704881e-15+1.32428553e-14j -1.26979923e-14-6.44907815e-15j 8.00000000e+00-1.38564065e+01j -1.60045167e-14-5.72145760e-15j 6.74553348e-15-1.31875814e-14j 5.67555423e-15+3.51142572e-15j 1.26565425e-14-4.44089210e-16j 8.12534275e-15+1.05483352e-14j -6.80594692e-15+9.73035361e-15j -1.00271283e-14+1.05969226e-14j -8.00000000e+00-1.38564065e+01j -1.25130829e-14-1.17522219e-15j 2.83430252e-15-2.54524974e-16j 4.14018138e-16-2.48398958e-15j]
Dakle, svega 4 koeficijenta imaju vrednosti realnog ili imaginarnog dela veće od $10^{-13}$. S obzirom na to, rezultate dobijene za faznu karakteristiku u onim tačkama u kojima je amplitudska karakteristika izuzetno mala ne možemo smatrati pouzdanim, zapravo moramo ih smatrati netačnim. Zato ćemo sve odbirke fazne karakteristike u onim tačkama gde je amplitudska karakteristika izuzetno mala, npr. manja od $10^{-13}$, postaviti na 0.
Xp = np.angle(X) # originalna fazna karakteristika
Xa = np.abs(X) # amplitudska karakteristika
Xp = Xp*(Xa > 1e-13) # ispravljena fazna karakteristika
fig = plt.figure(figsize = [6, 3])
plt.stem(Xp);
plt.subplots_adjust(bottom=0.25, left=0.15)
plt.xlabel(r'$k$')
plt.ylabel(r'$\angle(X[k])$');
DFT koeficijent $X[k]$ sa indeksom $k$ govori o doprinosu bazisnog signala oblika kompleksne sinusoide čija je učestanost
$$ \Omega_k = \frac{2\pi}{N}k. $$Što je veća vrednost ovog koeficijenta, to je veći udeo sinusoidalne komponente na učestanosti $\Omega_k$ u signalu $x[n]$.
Kompleksna sinusodida, tj. eksponencijalna funkcija sa imaginarnim argumentom, ima osobinu rotacione simetrije kružne učestanosti:
$$ e^{j\Omega n} = e^{j\Omega n} \cdot 1 = e^{j\Omega n} \cdot e^{-j2\pi n} = e^{j(\Omega - 2\pi) n} = e^{-j(2\pi - \Omega) n}, $$što znači da je pozitivna učestanost $\Omega$ koja se nalazi u opsegu od $\pi$ do $2\pi$ ekvivalentna negativnoj učestanosti $\Omega - 2\pi$. Ovo znači da polovina DFT koeficijenata odgovara negativnim učestanostima. S obzirom na sve navedeno, kada se fokusiramo na fizičko značenje koeficijenata DFT-a, ima smisla prikazivati koeficijente za opseg učestanosti od $-\pi$ do $\pi$. Dakle, koeficijente koji odgovaraju negativnim učestanostima $\Omega - 2\pi$, a to su koeficijenti počev od $k = N/2$ do $k = N-1$, treba premestiti na početak sekvence koeficijenata. Na ovaj način koeficijent za $k=0$ postaje centralni koeficijent i predstavlja doprinost jednosmerne komponente signala.
Dakle, izmenom redosleda koeficijenata smo samo iskoristili svojstvo periodičnosti Furijeove transformacije (koja je periodična sa periodom $2\pi$). U sekvenci koeficijenata sada prvi koeficijenti odgovaraju negativnim učestanostima, a druga polovina pozitivnim. Ovo je posebno zgodno kada radimo sa signalima čiji su svi odbirci realni brojevi jer ti signali imaju simetričan spektar oko učestanosti 0.
Primetite da je koeficijent počev od koga smatramo da je učestanost negativna $k=N/2$. Međutim, šta ako je $N$ neparan broj? Pošto je koeficijent za $k=0$ centralni, sekvence neparne dužine će biti preuređene u simetrične sekvence sa $(N-1)/2$ tačaka levo i desno od centralnog koeficijenta. Sekvence sa parne dužine će biti asimetrične sa jednim koeficijentom više na levoj strani.
Možemo napisati funkciju koja menja redosled koeficijenata:
import numpy as np
def dft_shift(X):
N = len(X)
if (N % 2 == 0):
return np.concatenate((X[(N//2):], X[:(N//2)])) # // integer div
else:
return np.concatenate((X[((N+1)//2):], X[:((N+1)//2)]))
Nacrtajmo amplitudsku karakteristiku jednog proizvoljnog signala u originalnom rasporedu.
np, mpl, plt = importEssentialLibs(USE_WIDGETS)
import scipy.fft as fft
x = np.arange(0, 1.02, 0.02) - 0.5
fig = plt.figure(figsize = [6, 3])
plt.subplots_adjust(bottom=0.25, left=0.15)
plt.stem(x);
plt.xlabel(r'$n$')
plt.ylabel(r'$x[n]$')
fig = plt.figure(figsize = [6, 3])
X = fft.fft(x);
plt.stem(abs(X));
plt.subplots_adjust(bottom=0.25, left=0.15)
plt.xlabel(r'$k$')
plt.ylabel(r'$|X[k]|$');
Sada možemo obrnuti redosled koeficijenata DFT-a. Za te potrebe možemo korisiti funkciju koju smo napisali, ali i funkciju iz scipy biblioteke:
X = fft.fftshift(X)
Videćemo da je rezultat potpuno isti.
fig = plt.figure(figsize = [6, 3])
plt.stem(abs(dft_shift(X)))
plt.subplots_adjust(bottom=0.25, left=0.15)
plt.xlabel(r'$k$')
plt.ylabel(r'$|X[k]|$')
fig = plt.figure(figsize = [6, 3])
plt.stem(abs(fft.fftshift(X)))
plt.subplots_adjust(bottom=0.25, left=0.15)
plt.xlabel(r'$k$')
plt.ylabel(r'$|X[k]|$');
Bilo bi mnogo razumljivije kada bismo negativne učestanosti predstavili negativnim koeficijentima, pa stoga moramo promeniti x-osu stem
naredbe u skladu sa dužinom DFT-a.
N = len(X)
k = np.arange(N) - (N//2)
fig = plt.figure(figsize = [6, 3])
plt.stem(k, abs(fft.fftshift(X)))
plt.subplots_adjust(bottom=0.25, left=0.15)
plt.xlabel(r'$k$')
plt.ylabel(r'$|X[k]|$');
Ranije smo videli da koeficijenti $X[k]$ pokazuju doprinos kompleksne sinusoide relativne kružne učestanosti $\Omega_k =k \cdot 2\pi/N$. Ovo, naravno važi i kada smo koeficijente preuredili tako da indeksi odgovaraju i negativnim i pozitivnim učestanostima. Stoga na x-osi možemo crtati relativne kružne učestanosti, a ne indekse.
N = len(X)
k = np.arange(N) - (N//2)
Wk = 2*np.pi/N*k
fig = plt.figure(figsize = [6, 3])
plt.stem(Wk, abs(fft.fftshift(X)))
plt.subplots_adjust(bottom=0.25, left=0.15)
plt.xticks([-np.pi, -np.pi/2, 0, np.pi/2, np.pi])
plt.xlabel(r'$\Omega$')
plt.ylabel(r'$|X[k]|$');
Ukoliko radimo frekvencijsku analizu nekog signala koji je dobijen odabiranjem analognog signala, onda na x-osi ima mnogo više smisla prikazati analogne učestanosti koje odgovaraju relativnim učestanostima. Podsetimo se, ako je učestanost odabiranja $f_s$, onda je relativna učestanost $F=f/f_s$, a relativna kružna učestanost $\Omega = 2\pi F = 2\pi f/f_s$. Stoga možemo mapirati relativne kružne učestanosti u učestanosti analognog signala na sledeći način:
$$ f = fs\cdot F = \frac{f_s\cdot \Omega}{2\pi} = \frac{f_s \cdot k 2\pi/N}{2\pi} = f_s\frac{k}{N} $$Na osnovu prethodne analize, možemo uraditi frekvencijsku analizu jednog realnog signala. To će biti snimljeni zvuk sa klavira. Za učitavanje ćemo koristiti IPython i scipy.
np, mpl, plt = importEssentialLibs(USE_WIDGETS)
import scipy.fft as fft
import IPython
from scipy.io import wavfile
fs, x = wavfile.read("audio/01_guitar_E2.wav")
IPython.display.display(IPython.display.Audio(x, rate = fs))
N = len(x)
X = fft.fft(x)
k = np.arange(N) - N//2
f = fs*k/N
fig = plt.figure(figsize = [6, 5])
plt.plot(f, abs(fft.fftshift(X)))
plt.subplots_adjust(bottom=0.25, left=0.15)
plt.xlabel(r'$f$ [Hz]')
plt.ylabel(r'$|X[k]|$');
Iz spektra signala vidimo da su sve značajne frekvencijske komponente manje od 2000 Hz, pa nema mnogo smisla crtati ceo opseg učestanosti. Dodatno, signal je realan, pa je spektar simetričan oko nule. Iz tog razloga dovoljno je nacrtati samo jednu stranu spektra. Nacrtajmo samo deo spektra. Sledeća ćelija koristi ipywidgets biblioteku za učitavanje imena fajla koji se učitava. Za gitaru, spektar se crta do učestanosti 2000 Hz, dok se za ostale instrumente crta do učestanosti 6000 Hz. Dodatno, na odvojenoj slici se crta i vremenski oblik dela učitanog signala.
np, mpl, plt = importEssentialLibs(USE_WIDGETS)
import scipy.fft as fft
import IPython
from scipy.io import wavfile
import ipywidgets as widgets
def loadSignalAndPlotSpectra(filename):
fs, x = wavfile.read("audio/"+filename)
IPython.display.display(IPython.display.Audio(x, rate = fs))
N = len(x)
X = fft.fft(x)
if "guitar" in filename:
fMaxShow = 2000
else:
fMaxShow = 6000
Nmax = fMaxShow*N//fs
f = fs*np.arange(Nmax)/N
# only part of the spectra
X = X[:len(f)];
# Plot spectra
fig = plt.figure(figsize = [6, 4])
plt.plot(f, abs(X))
plt.subplots_adjust(bottom=0.15, left=0.15)
plt.xlabel(r'$f$ [Hz]')
plt.ylabel(r'$|X[k]|$');
plt.show()
# Plot signal
fig = plt.figure(figsize = [6, 4])
n = np.arange(N//10, N//10 + 2000)
t = n/fs
plt.plot(t, x[n])
plt.subplots_adjust(bottom=0.15, left=0.15)
plt.xlabel(r'$t$ [s]')
plt.ylabel(r'$x(t)$');
widgets.interact(loadSignalAndPlotSpectra, filename=["01_guitar_E2.wav",
"02_guitar_A2.wav",
"03_guitar_D3.wav",
"04_guitar_G3.wav",
"05_guitar_H3.wav",
"06_guitar_E4.wav",
"07_clarinet_A.wav",
"08_violin_H.wav"]);
Do sada smo u primerima birali relativnu kružnu učestanost signala tako da je ona na $\Omega_k = k \cdot 2\pi/N$, tj. učestanost koja je i učestanost nekod od bazisnih signala. Tako smo uvek dobijali da u $N$ odbiraka signala imamo ceo broj perioda kompleksnih sinusoida i u spektru smo dobijali nenulte komponente samo na indeksima koji odgovaraju relativnim učestanostima $\Omega_k$. Međutim, šta se dešava ako je relativna kružna učestanost kompleksne sinusoide takva da nije celobrojni umnožak $2\pi/N$? Ne postoji indeks $k$ koji odgovara ovoj učestanosti. Da li to znači da ne možemo prikazati DFT ovog signala? Apsolutno ne, štaviše, probajmo u Pajtonu.
np, mpl, plt = importEssentialLibs(USE_WIDGETS)
import scipy.fft as fft
N = 16
n = np.arange(N)
# relativne ucestanosti
F = 2.42/N
#pocetne faze
Phi = 0 # np.pi/3, 2*np.pi/3
x = np.exp(1j*(2*np.pi*F*n + Phi))
X = fft.fft(x)
fig = plt.figure(figsize = [6, 3])
X = fft.fft(x);
plt.stem(abs(X));
plt.subplots_adjust(bottom=0.25, left=0.15)
plt.xlabel(r'$k$')
plt.ylabel(r'$|X[k]|$');
fig = plt.figure(figsize = [6, 3])
X = fft.fft(x);
plt.stem(np.angle(X));
plt.subplots_adjust(bottom=0.25, left=0.15)
plt.xlabel(r'$k$')
plt.ylabel(r'$\angle(X[k])$');
Primećujemo da sada umesto jednog nenultog koeficijenta DFT-a, imamo sve koeficijente koji su različiti od nule! Ovo je posledica Parsevalove teoreme koja kaže da se energija signala u vremenskom domenu mora održati i u frekvencijskom, a koja važi i za diskretnu Furijeovu transformaciju. Da li ovaj spektar ima nekog smisla? Pogledajmo dobijeni rezultat. Koeficijenti koji su najveći su koeficijenti koji odgovaraju učestanostima $\Omega_k$ koje su najbliže ulaznoj učestanosti $\Omega$. Dakle, nije baš da nemaju smisla. Mi u opštem slučaju ne možemo da predstavimo kompleksnu sinusoidu na učestanosti $\Omega$ samo jednim od bazisnih vektora DFT-a, ali možemo kao linearnu kombinaciju više njih. Tako, možemo da zaključimo da ulazna sinusoida nije tačno na učestanosti nekog od bazisnih signala, ali možemo da zaključimo u blizini kojih učestanosti jeste.
Dodatno, setimo se šta je DFT. DFT je Furijeova transformacija signala odabirana u ekvidistantnim tačkama. Signal ima svega $N$ nenultih odbiraka. Svi ostali odbirci su nula. Možemo nacrtati Furijeovu transformaciju ovog signala.
np, mpl, plt = importEssentialLibs(USE_WIDGETS)
import scipy.fft as fft
def fourierTransform(x, N):
# Ova funkcija izračunava Furijeovu transformaciju u N ekvidistantnih tačaka
if len(x) < N:
x = np.append(x, np.zeros(N - len(x)))
X = np.zeros(N, dtype = 'complex')
n = np.arange(N)
for k in range(N):
Omega_k = 2*np.pi/N*k
X[k] = np.sum(x*np.exp(-1j*Omega_k*n))
return X
N = 16
n = np.arange(N)
# relativne ucestanosti
F = 3.5/N
#pocetne faze
Phi = 0 # np.pi/3, 2*np.pi/3
x = np.exp(1j*(2*np.pi*F*n + Phi))
NFT = 1024;
X = fourierTransform(x, NFT)
Omega = np.arange(NFT)/NFT*2*np.pi
fig = plt.figure(figsize = [6, 3])
plt.plot(Omega, abs(X));
plt.subplots_adjust(bottom=0.25, left=0.15)
plt.xlabel(r'$\Omega$')
plt.ylabel(r'$|X(e^{j\Omega})|$');
plt.xlim([0, 2*np.pi])
plt.ylim([0, max(abs(X))*1.05]);
Kao što smo videli, DFT signala dužine $N$ je vektor odbiraka Furijeove transformacije u $N$ ekvidistantnih tačaka. Napišitmo to u Pajtonu.
np, mpl, plt = importEssentialLibs(USE_WIDGETS)
def fourierTransform(x, N):
# Ova funkcija izračunava Furijeovu transformaciju u N ekvidistantnih tačaka
if len(x) < N:
x = np.append(x, np.zeros(N - len(x)))
X = np.zeros(N, dtype = 'complex')
n = np.arange(N)
for k in range(N):
Omega_k = 2*np.pi/N*k
X[k] = np.sum(x*np.exp(-1j*Omega_k*n))
return X
N = 16
n = np.arange(N)
# relativne ucestanosti
F = 3.1/N
#pocetne faze
Phi = 0 # np.pi/3, 2*np.pi/3
x = np.exp(1j*(2*np.pi*F*n + Phi))
NFT = 1024;
X = fourierTransform(x, NFT)
Omega = np.arange(NFT)/NFT*2*np.pi
Xk = X[::NFT//N]
Omega_k = Omega[::NFT//N]
fig = plt.figure(figsize = [6, 3])
plt.plot(Omega, abs(X))
plt.stem(Omega_k, abs(Xk), 'r', markerfmt='ro')
plt.subplots_adjust(bottom=0.25, left=0.15)
plt.xlabel(r'$\Omega$')
plt.ylabel(r'$|X(e^{j\Omega})|, \, |X[k]|$')
plt.xlim([0, 2*np.pi])
plt.ylim([0, max(abs(X))*1.05]);
Dakle, ništa se nije promenilo, DFT je i dalje Furijeova tranformacija izračunata u $N$ tačaka. Odabiranjem Furijeove transformacije je dobijen vektor koji smo dobili i ranije. Zaključujemo da ćemo imati jednu frekvencijsku komponentu od jedne kompleksne sinusoide samo ako je ta kompleksna sinusoida na učestanosti jednog od bazisnih vektora, tj. ako smo uzeli ceo broj perioda ove sinusoide. Ako ovo nije slučaj, javiće se takozvano curenje spektra. Signal moramo predstaviti koeficijentima koji su nam na raspolaganju. Korišćenjem sledećeg koda, možemo malo jasnije razumeti ovu pojavu. Dodatno, setimo se da je DTF praktično Furijeov red periodičnog produženja signala čiji DFT računamo. Sledeći kod crta i 3 periode tog periodičnog produženja. Šta uočavate?
np, mpl, plt = importEssentialLibs(USE_WIDGETS)
import ipywidgets as widgets
import scipy.fft as fft
fig, axs = plt.subplots(3, 1, figsize = [10, 10])
plt.subplots_adjust(bottom=0.25, hspace = 0.5)
def plotFTandDFT(F0, N):
n = np.arange(N)
x = np.exp(1j*(2*np.pi*F0*n))
NFT = 1024;
X = fourierTransform(x, NFT)
Omega = np.arange(NFT)/NFT*2*np.pi
F = Omega/2/np.pi
Xk = fft.fft(x)
Omega_k = np.arange(N)/N*2*np.pi
F_k = Omega_k/2/np.pi
print(f'k is: {F0*N}')
#Clear figures
axs[0].cla()
axs[1].cla()
axs[2].cla()
# Plot signal x
ax = axs[0]
ax.plot(x.real, 'o-', label='Real')
ax.plot(x.imag, 'o--', label = 'Imag')
ax.set_xlabel(r'$n$')
ax.set_ylabel(r'$x[n]$')
ax.legend(loc = 'upper right')
# Plot periodic expansion
xp = np.append(x, x)
xp = np.append(xp, x)
ax = axs[1]
ax.plot(xp.real, 'o-', label='Real')
ax.plot(xp.imag, 'o--', label = 'Imag')
ax.set_xlabel(r'$n$')
ax.set_ylabel(r'$x_p[n]$')
ax.legend(loc = 'upper right')
# Plot spectra
ax = axs[2]
ax.plot(F, abs(X))
ax.stem(F_k, abs(Xk), 'r', markerfmt='ro')
ax.set_xlabel(r'$F$')
ax.set_ylabel(r'$|X(e^{j\Omega})|, \, |X[k]|$')
ax.set_xlim([0, 1])
ax.set_ylim([0, max(abs(X))*1.05]);
widgets.interact(plotFTandDFT, F0 = (0, 1, 0.01), N = (10, 31, 5));
Videli smo da za različite kombinacije parametara N i F dobijamo da je mnogo veća verovatnoća da će do curenja spektra doći, nego da neće. Zamislite sada da nemamo samo jednu sinusoidalnu komponentu, već da u signalu imamo mnogo spektralnih komponenti na različitim učestanostima. U prirodi je ovo najčešći slučaj. U tim situacijama je gotovo nemoguće odabrati razumno malo $N$ da do curenja spektra ne dođe. Dodatno, mi DFT-om radimo frekvencijsku analizu signala. Stoga, najčešće ne znamo kakav je signal na ulazu i ne možemo unapred znati kakvo $N$ da odaberemo. Dakle, da budem slobodan da se izrazim: curenje spektra nam ne gine. :) Samo treba da vidimo kako da se sa njim borimo.
Ako posmatramo signal u svega 8 tačaka, u opštem slučaju je teško da ćemo iz njega moći da zaključimo mnogo šta. Koeficijenti DFT-a takvog signala će biti izračunati u malom broju tačaka i velika je verovatnoća da ćemo dobiti curenje spektra. Iz koeficijenata ćemo moći da zaključimo da se neka frekvencijska komponenta nalazi najverovatnije između dve učestanosti susednih bazisnih signala. Susedni bazisni signali su na razmaku od $2\pi/N$. Dakle, ne možemo odrediti učestanost bilo koje spektralne komponente sa rezolucijom manjom od $\Delta \Omega = 2\pi/N$. Ova vrednost se naziva frekvencijskom rezolucijom DFT-a.
Da bismo uradili bolju estimaciju učestanosti, moramo uzeti više odbiraka signala, što je lako uočiti na poslednjem primeru sa slajderima za $F$ i $N$.
S obzirom na to da je signal veće dužine takav da će frekvencijska rezolucija DFT-a biti veća, postavlja se pitanje da li se produženjem signala nulama može postići sličan efekat odabiranju u većem broju tačaka. Ako signal dužine $N_1$ produžimo nulama do dužine $N_2$, gde je $N_2 > N_1$, to nas može navesti da pretpostavimo i da će frekvencijska rezolucija biti $\Delta \Omega_2 = 2\pi/N_2 < \Delta \Omega_1 = 2\pi/N_1$. Ponovimo prethodnu analizu na signalu koji je produžen nulama tako da je broj tačaka duplo veći.
np, mpl, plt = importEssentialLibs(USE_WIDGETS)
import ipywidgets as widgets
import scipy.fft as fft
fig, axs = plt.subplots(3, 1, figsize = [10, 10])
plt.subplots_adjust(bottom=0.25, hspace = 0.5)
def plotFTandDFT(F, N):
n = np.arange(N)
x = np.exp(1j*(2*np.pi*F*n))
# PRODUZENJE NULAMA
x = np.append(x, np.zeros(N))
N = len(x) # NOVO N
NFT = 1024;
X = fourierTransform(x, NFT)
Omega = np.arange(NFT)/NFT*2*np.pi
Xk = fft.fft(x)
Omega_k = np.arange(N)/N*2*np.pi
print(f'k is: {F*N}')
#Clear figures
axs[0].cla()
axs[1].cla()
axs[2].cla()
# Plot signal x
ax = axs[0]
ax.plot(x.real, 'o-', label='Real')
ax.plot(x.imag, 'o--', label = 'Imag')
ax.set_xlabel(r'$n$')
ax.set_ylabel(r'$x[n]$')
ax.legend(loc = 'upper right')
# Plot periodic expansion
xp = np.append(x, x)
xp = np.append(xp, x)
ax = axs[1]
ax.plot(xp.real, 'o-', label='Real')
ax.plot(xp.imag, 'o--', label = 'Imag')
ax.set_xlabel(r'$n$')
ax.set_ylabel(r'$x_p[n]$')
ax.legend(loc = 'upper right')
# Plot spectra
ax = axs[2]
ax.plot(Omega, abs(X))
ax.stem(Omega_k, abs(Xk), 'r', markerfmt='ro')
ax.set_xlabel(r'$\Omega$')
ax.set_ylabel(r'$|X(e^{j\Omega})|, \, |X[k]|$')
ax.set_xlim([0, 2*np.pi])
ax.set_ylim([0, max(abs(X))*1.05]);
widgets.interact(plotFTandDFT, F = (0, 1, 0.01), N = (10, 31, 5));
Možemo primetiti da je jedino što smo postigli da smo Furijeovu transformaciju istog signala odabirali u većem broju tačaka. Prividno, frekvencijska rezolucija jeste poboljšana, ali periodično produženje signala više nikada nije periodična kompleksna sinusoida jer između nenultih odbiraka imamo čitav niz nula. Za konkretan primer, JEDNE kompleksne sinusoide, možemo dobiti bolju informaciju o učestanosti na kojoj je najveća frekvencijska komponenta. Međutim, i ovo je IZUZETNO VAŽNO, čim se u spektru pojavi makar još jedna frekvencijska komponenta (što je uvek slučaj u praktičnim signalima), Furijeova transofrmacija druge kompleksne sinusoide odabirane u $N$ tačaka će se sabrati sa Furijeovom transformacijom prve kompleksne sinusoide. Rezultujuča Furijeova transformacija je ona koja se odabira i tada se može desiti da se u zbiru pomere maksimumi Furijeove transformacije! Dakle, dopunjavanjem nulama se samo dobija precizniji prikaz Furijeove transformacije signala odabranog u $N$ tačaka i čiji su odbirci 0 za sve $n$ koji nisu između $0$ i $N-1$. Da bi se povećala Frekvencijska rezolucija, mora se uzeti veći broj odbiraka signala. Dakle, frekvencijsku rezoluciju određujemo pri ODABIRANJU.
Sledeći primer treba da pokaže da, iako smo dopunili nulama zbir dve kompleksne sinusoide odabrane u 10 tačaka, uvek može da se desi da najveći koeficijent ne odgovara učestanosti koja je najbliža pravoj učestanosti (prave učestanosti su obležene crnim isprekidanim linijama). Na primer, odaberite slučaj $F_1 = 0.24$ i $F_2 = 0.14$.
np, mpl, plt = importEssentialLibs(USE_WIDGETS)
import ipywidgets as widgets
import scipy.fft as fft
fig = plt.figure(figsize = [10, 5])
def plotFTandDFT(F1, F2):
plt.clf()
N = 10
n = np.arange(N)
x1 = np.exp(1j*(2*np.pi*F1*n))
x2 = np.exp(1j*(2*np.pi*F2*n))
# PRODUZENJE NULAMA
x1 = np.append(x1, np.zeros(N))
x2 = np.append(x2, np.zeros(N))
x = x1 + x2
N = len(x) # NOVO N
NFT = 1024;
X = fourierTransform(x, NFT)
X1 = fourierTransform(x1, NFT)
X2 = fourierTransform(x2, NFT)
Omega = np.arange(NFT)/NFT*2*np.pi
Xk = fft.fft(x)
Omega_k = np.arange(N)/N*2*np.pi
print(f'k is: {F*N}')
# Plot spectra
plt.plot(Omega, abs(X), label = '$|X(e^{j\Omega})|$')
plt.plot(Omega, abs(X1), '--', label = '$|X_1(e^{j\Omega})|$')
plt.plot(Omega, abs(X2), '--', label = '$|X_2(e^{j\Omega})|$')
plt.stem(Omega_k, abs(Xk), 'r', markerfmt='ro', label = '$|X[k]|$')
plt.legend(loc = 'upper right', ncol = 2)
# stem original freqs
if F1 == F2:
plt.stem(np.array([F1])*2*np.pi, np.array([N]), 'k--', markerfmt='ks')
else:
plt.stem(np.array([F1, F2])*2*np.pi, np.array([N/2, N/2]), 'k--', markerfmt='ks')
plt.subplots_adjust(bottom=0.25, left=0.15)
plt.xlabel(r'$\Omega$')
plt.ylabel(r'$|X(e^{j\Omega})|, \, |X[k]|$')
plt.xlim([0, 2*np.pi])
plt.ylim([0, N]);
widgets.interact(plotFTandDFT, F1 = (0, 0.3, 0.02), F2 = (0, 0.3, 0.01));
<Figure size 720x360 with 0 Axes>
ZADATAK 1: Generisati proizvoljan realan signal kod koga nije zadavoljeno signhrono odabiranje. Izračunati DFT ovog signala, izdvojiti njegov realni i imaginarni deo, izračunati njegov fazni i amplitudski spektar, nacrtati ih, a zatim proveriti sledeće relacije:
$$ X_r[k] = X_r[N - k] $$$$ X_i[k] = -X_i[N - k] $$$$ |X[k]| = |X[N - k]| $$$$ \angle (X[k]) = -\angle (X[N - k]) $$
ZADATAK 2: Generisati signal $x[n]=\cos(2\pi F_0 n) + \cos(2\pi F_1 n)$ za proizvoljne $F_0$ i $F_1$, u proizvoljnom broju tačaka $N$ i naći njegovu DFT. U kakvoj vezi treba da budu $F_0$, $F_1$ i $N$ da ne bi došlo do curenja spektra?