Show 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
LabView program za zajem signalov lahko prenesete v obliki zip arhiva
.
Pripravite kratko poročilo (obseg od 3 do 10 celic s kodo v okolju Jupyter Notebook), iz katerega naj bodo razvidni podatki naloge (iz tabele), ter da ste vse parametre in zahteve naloge naloge tudi upoštevali. Vsebuje naj prikaz in kratek komentar naslednjih lastnosti zajetega signala:
zajeto frekvenčno območje in frekvenčna ločljivost,
ocena gostote porazdelitve verjetnosti vrednosti
zajetega signala ,enostranska ter dvostranska avtospektralna gostota moči,
in ,prvi štirje centralni statistični momenti (
, , , ) 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,
Note
Diskretna Fourierova transformacija:
Velja
Ker je DFT periodična z
Modul numpy.fft
izračuna
Primer 1
Poglejmo diskretizacijo frekvenčnega vektorja pri DFT signala, vzorčenega pri časih:
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
Lastnosti takih signalov je pogosto smiselno obravnavati v frekvenčni domeni.
Primer 2
Normalno porazdeljena naključna spremenljivka
Ocenimo prva dva centralna statistična momenta funkcije gostote verjetnosti (
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
5.186961971048731e-16
s_2 = np.mean((x_center)**2)
s_2
3.9035971550047157
Bolj splošen generator naključnih vrednosti v Numpy:
rng = np.random.default_rng()
y = rng.normal(loc=mu, scale=sigma, size=N)
Show 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()

Porazdelitev verjetnosti ni nujno normalna (Gaussova):
z = rng.weibull(a=2.5, size=N)
Show 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$ [/]');

Show 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()

Note
Naključni proces
Naključni proces
Primer 3
Opazujemo naključni proces
Opazujmo razvoj srednje vrednosti in variance signala v času,
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)
Show code cell source
for p in ansambel:
plt.plot(t, p, '.')
plt.xlabel('t')
plt.ylabel('x(t) (N realizacij)');

mu_t = np.mean(ansambel, axis=0)
var_t = np.mean((ansambel-np.mean(ansambel, axis=0))**2, axis=0)
Show 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)$');

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)
Show 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)$');

Ker se srednja vrednost in varianca naključnega procesa
Note
Ato-korelacijska funkcija
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:
Za ergodičen proces jo lahko izračunamo kot:
Primer 4
Opazujemo naključni proces
kjer je
Opazujmo avtokorelacijsko funkcijo
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)
Show code cell source
plt.plot(np.hstack([-t[::-1], t[1:]]), R_aa)
plt.xlabel('t [s]')
plt.ylabel('$R_{aa}$');

R_aa = R_aa[-len(a):]
Show code cell source
plt.plot(t, R_aa)
plt.xlabel('t [s]')
plt.ylabel('$R_{aa}$');

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))
Show 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');

Parsevalov teorem:
a.var(), np.trapz(np.abs(G_aa), freq)
(0.78897363152495, 0.7893752336562335)