Ctrl+K

DFT in osnove naključnih procesov

Hide code cell source 
import numpy as np
from scipy import stats, signal
import matplotlib.pyplot as plt

DFT in osnove naključnih procesov#

Domača naloga#

Naloga

Z uporabo generatorja signalov in zajemnega sistema Arduino zajemite naključni signal s parametri (srednjo vrednostjo \(m\), amplitudo \(A\)) iz podatkov naloge. Zajem podatkov izvedite tako, da boste pri pretvorbi zajetega signala v frekvenčno domeno dobili spekter na podanem frekvenčnem območju \([0, f_k]\) ter z ločljivostjo v frekvenčni domeni \(\Delta f\) iz podatkov naloge.

LabView program za zajem signalov lahko prenesete v obliki zip arhiva.

Pripravite kratko poročilo v okolju Jupyter Notebook (od 3 do 10 celic s kodo), v katerem naj bodo razvidni podatki naloge ter prikaz in kratek komentar naslednjih lastnosti zajetega signala:

  • zajeto frekvenčno območje in frekvenčna ločljivost,

  • ocena gostote porazdelitve verjetnosti vrednosti \(p(x)\) zajetega signala \(x(t)\),

  • enostranska ter dvostranska avtospektralna gostota moči, \(G_{XX}\) in \(S_{XX}\),

  • prvi štirje centralni statistični momenti (\(m_0\), \(m_1\), \(m_2\), \(m_3\)) zajetega signala.

Dodatek:

  • Ovrednotite in komentirajte stacionarnost ter ergodičnost zgoraj generiranega naključnega signala (po potrebi zajemite dodatne signale, ki jih za to potrebujete).

  • Numerično generirajte in komentirajte primera nestacionarnega gaussovega ter stacionarnega ne-gaussovega naključnega signala (glejte na primer pyExSi).

Poročilo oddajte v .pdf obliki (glejte navodila za oddajo domačih nalog).

Diskretna Fourierova transformacija#

Računamo Fourierovo transformacijo diskretnega signala, \(x_n = x(n\,\Delta t)\). Kadar imamo končno mnogo vzorčenih podatkov (\(N\)) in gre \(n=[0,1,\dots,N-1]\), uporabimo diskretno Fourierovo transformacijo (DFT).

Note

Diskretna Fourierova transformacija:

\[ X_k = \sum_{n=0}^{N-1} x_n\,e^{-\mathrm{i}\,2\pi\,k\,n/N}. \]

Velja \(X_k=X(k/(N\,\Delta t))\).

Ker je DFT periodična z \(1/\Delta t\) (\(X_k=X_{k+N}\)), je treba izračunati samo \(N\) členov.

Modul numpy.fft izračuna \(X_k\) za diskretne frekvenčne točke pri \(k \in [-N/2, \dots, 0,\dots\,N/2-1]\) (v primeru sodega števila točk \(N\)).

Primer 1

Poglejmo diskretizacijo frekvenčnega vektorja pri DFT signala, vzorčenega pri časih: \(t = [0, 1, 2, \dots, 9] s\):

t = np.arange(0, 10)
delta_t = t[1] - t[0]
n = len(t)
delta_t, n

(1, 10)

np.fft.fftfreq(len(t), t[1]-t[0])

array([ 0. ,  0.1,  0.2,  0.3,  0.4, -0.5, -0.4, -0.3, -0.2, -0.1])

np.fft.rfftfreq(len(t), t[1]-t[0])

array([0. , 0.1, 0.2, 0.3, 0.4, 0.5])

Osnove naključnih procesov#

Procese, ki jih ne moremo obravnavati kot deterministične, pogosto modeliramo kot naključne procese.

Za te je značilno, da njihovih dejanskih vrednosti v času \(x(t)\) ne moremo natančno predvideti, lahko pa z določeno verjetnostjo sklepamo o porazdelitvi njihovih vrednosti. Tako definiramo funkcijo gostote porazdelitve verjetnost, \(p(x)\), ki podaja verjetnost, da se bo v vzorcu naključne spremenljivke \(x(t)\) pojavila določena vrednost \(x\).

Lastnosti takih signalov je pogosto smiselno obravnavati v frekvenčni domeni.

Primer 2

Normalno porazdeljena naključna spremenljivka \(x\) z \(N=1000\) vzorci, srednjo vrednostjo \(\mu = 5\) in standardno deviacijo \(\sigma = 2\) predstavlja realizacijo N ponovitev naključnega dogodka.

Ocenimo prva dva centralna statistična momenta funkcije gostote verjetnosti (\(n = [1, 2]\)):

\[ \sigma^n(t)= E[(x(t)-\mu_x)^n] = \int_{-\infty}^{+\infty} (x-\mu_x)^n\, p(x)\, \textrm{d} x \]
N = 1000
mu = 5
sigma = 2

x = np.random.randn(N)*sigma + mu

x_center = x - np.mean(x)

s_1 = np.mean(x_center)
s_1

-4.050093593832571e-16

s_2 = np.mean((x_center)**2)
s_2

3.9616020921751622

Bolj splošen generator naključnih vrednosti v Numpy:

rng = np.random.default_rng()
y = rng.normal(loc=mu, scale=sigma, size=N)
Hide code cell source 
hist, bins, _ = plt.hist(x, bins=50, alpha=0.5, density=True, label='$x$');
hist, bins, _ = plt.hist(y, bins=50, alpha=0.5, density=True, label='$y$');
plt.plot(bins, stats.norm.pdf(bins, mu, sigma), 'k', lw=2, label='$p(x)$')
plt.legend()
plt.xlabel('$x$')
plt.ylabel('$p(x)$ [/]');
plt.show()
../_images/263113ccdb1eca540ca94c4619869f569511294c8e48219fb8a07045057be154.png

Porazdelitev verjetnosti ni nujno normalna (Gaussova):

z = rng.weibull(a=2.5, size=N)
Hide code cell source 
plt.plot(x, alpha=0.75, label='$x$ (Gauss)')
plt.plot(y, alpha=0.75, label='$y$ (Gauss)')
plt.plot(z, alpha=0.75, label='$z$ (Weibull)')
plt.legend()
plt.xlabel('vzorec $i$ [/]')
plt.ylabel('$x_i$ [/]');
../_images/820a29dfab24a53e0e4fbde811832f22003e1972cac71eda14c632fbc17bb5cb.png
Hide code cell source 
hist, bins, _ = plt.hist(x, bins=50, alpha=0.5, density=True, label='$x$');
hist, bins, _ = plt.hist(y, bins=50, alpha=0.5, density=True, label='$y$');
hist_w, bins_w, _ = plt.hist(z, bins=50, alpha=0.5, density=True, label='$z$');
plt.plot(bins, stats.norm.pdf(bins, mu, sigma), 'k', lw=2, label='$p(x)$ (normal)')
plt.plot(bins_w, stats.weibull_min.pdf(bins_w, 2.5), 'r', lw=2, label='$p(z)$ (Weibull)')
plt.legend()
plt.xlabel('$x$')
plt.ylabel('$p(x)$ [/]');
plt.show()
../_images/e0240cd361a154098cf83b3d53662dd266c64f94e102ec84c4ca4639c56514a5.png

Note

Naključni proces \(x(t)\) je stacionaren, če so srednja vrednost in kovariančne funkcije procesa časovno neodvisne.

Naključni proces \(x(t)\) je ergodičen, če so statistične lastnosti nižjega reda ene same realizacije procesa v času (\(x(t)\)) enake statističnim lastnostim ansambla več realizacij dogodka \(x_k(t)\). Pravimo, da je posamezna realizacija procesa reprezentativna.

Primer 3

Opazujemo naključni proces \(x(t)\) z normalno porazdelitvijo verjetnosti (\(\mu = -4\), \(\sigma = 1\)), definiran v 25 časovnih točkah \(t = [0, 1, \dots, 24]\). Posnamemo ansambel \(N=50\) realizacij naključnega procesa.

Opazujmo razvoj srednje vrednosti in variance signala v času, \(\mu_x(t)\) in \(\sigma^2_x(t)\).

t = np.arange(25)
N = 50
mu = -4
sigma = 1

ansambel = []
for ponovitev in range(N):
    ansambel.append(rng.normal(loc=mu, scale=sigma, size=len(t)))
ansambel = np.array(ansambel)
Hide code cell source 
for p in ansambel:
    plt.plot(t, p, '.')
plt.xlabel('t')
plt.ylabel('x(t) (N realizacij)');
../_images/a7700e76e45059a78ff93df142417d758b5e2322940478f9eed08cc26aaa6a1d.png 
mu_t = np.mean(ansambel, axis=0)
var_t = np.mean((ansambel-np.mean(ansambel, axis=0))**2, axis=0)
Hide code cell source 
fig, ax = plt.subplots(2, 1)
ax[0].plot(t, mu_t)
ax[0].set_ylim(mu-sigma, mu+sigma)
ax[0].set_xlabel('t')
ax[0].set_ylabel('$\\mu_x(t)$');

ax[1].plot(t, var_t)
ax[1].set_ylim(sigma-sigma, sigma+sigma)
ax[1].set_xlabel('t')
ax[1].set_ylabel('$\\sigma^2_x(t)$');
../_images/4b34d3fb03204b08666f95631c93e7377cd918f42f44efe5b5a5a1d5ee3eb5c4.png

Ker se srednja vrednost in varianca naključnih realizacij v isti točki procesa s časom ne spreminjata, lahko sklepamo, da gre za stacionaren proces.

mu_ansambel = np.mean(ansambel, axis=1)
var_ansambel = np.mean((ansambel-np.mean(ansambel, axis=0))**2, axis=1)
Hide code cell source 
fig, ax = plt.subplots(2, 1)
x_ansambel = np.arange(len(mu_ansambel))
ax[0].plot(x_ansambel, mu_ansambel, 'k')
ax[0].set_ylim(mu-sigma, mu+sigma)
ax[0].set_ylabel('$\\mu_x(k)$');

ax[1].plot(x_ansambel, var_ansambel, 'k')
ax[1].set_ylim(sigma-sigma, sigma+sigma)
ax[1].set_xlabel('Realizacija procesa $k$')
ax[1].set_ylabel('$\\sigma^2_x(k)$');
../_images/ce1a50ab07857c7363af939ed6aaf2ae087f9874028107fa603f7beb57146ba0.png

Ker se srednja vrednost in varianca naključnega procesa \(x(t)\) pri večkratni realizaciji ne spreminjata, in sta enaki \(\mu_x(t)\) in \(\sigma^2_x(t)\) pri posameznem času, lahko sklepamo, da gre za ergodičen proces.

Note

Ato-korelacijska funkcija \(R_{xx}(\tau)\) je definirana z:

\[R_{xx}(\tau) = E[ x_k(t) \, x_k(t+\tau)]\]

in je enaka avto-kovariančni funkciji v primeru procesa z ničelno srednjo vrednostjo.

Note

Avto-spektralna gostota moči (PSD, Power Spectral Density) je Fourierova transformacija avtokorelacijske funkcije signala:

\[S_{xx}(f) = \mathcal{F}\big(R_{xx}(\tau)\big)\]

Za ergodičen proces jo lahko izračunamo kot:

\[S_{xx}(f) = \frac{1}{T}\, X^*(f) \, X(f)\qquad\text{PSD}.\]

Primer 4

Opazujemo naključni proces \(a(t)\), definiran z:

\[a(t) = \sin(10 \, \pi \, t) + n(t),\]

kjer je \(n(t)\) naključni šum z normalno porazdelitvijo (\(\mu=0\), \(\sigma=0.5\)).

Opazujmo avtokorelacijsko funkcijo \(R_{aa}(t)\) in gostoto spektralne moči \(S_{aa}(f)\) signala \(a(t)\).

f = 5
T = 1
fs = 100
t = np.arange(fs*T) / fs
a = np.sin(2*np.pi*f*t) + np.random.randn(len(t))*0.5

R_aa = np.correlate(a-np.mean(a), a-np.mean(a), mode='full') / np.std(a)**2 / len(a)
Hide code cell source 
plt.plot(np.hstack([-t[::-1], t[1:]]), R_aa)
plt.xlabel('t [s]')
plt.ylabel('$R_{aa}$');
../_images/c38bd61808799ba671adb7d37f923af0a6acc09025d62ee46fe7b41f90961923.png 
R_aa = R_aa[-len(a):]
Hide code cell source 
plt.plot(t, R_aa)
plt.xlabel('t [s]')
plt.ylabel('$R_{aa}$');
../_images/43589b0ae5c41e56109d613bba4b20079ab98dad626a0a68e8f38492a402e3f5.png 
G_aa_ = np.fft.rfft(R_aa) / len(R_aa) * 2 * T

A = np.fft.rfft(a) / len(a)
freq = np.fft.rfftfreq(len(a), 1/fs)
G_aa = np.conj(A) * A  * 2 * T

Izračun PSD v Scipy:

from scipy.signal import welch
f_scipy, G_aa_scipy = welch(a, fs=fs, scaling='density', window=np.ones_like(a))
Hide code cell source 
plt.plot(freq, np.abs(G_aa_), label='$\\mathcal{F}(R_{aa})$')
plt.plot(freq, np.abs(G_aa), label='$\\mathcal{F}(A)* \\cdot \\mathcal{F}(A) * T$')
plt.plot(freq, G_aa_scipy, '--', label='PSD v Scipy')
plt.legend()
plt.xlabel('f [Hz]')
plt.ylabel('PSD');
../_images/87607583bc06f243c6cf4c8fe67d066aeed6e1139a7068b4c9b406f1c26ace5e.png

Parsevalov teorem:

a.var(), np.trapz(np.abs(G_aa), freq)

(0.839728121091446, 0.8398123217642351)
Contents