Kompjutaciona fizika - seminarski rad¶

Dopplerov efekat je pojava s kojom su upoznati skoro svi, bilo svjesno ili nesvjesno. Moderna prevozna sredstva s kojim se susrećemo svakodnevno omogućavaju stalnu demonstraciju ove pojave. Širok je spektar aplikacija Dopplerovog efekta, a u ovom seminarskom radu će biti pokazano kako se Dopplerov efekat uz korištenje Furijeove analize može primjeniti za određivanje brzine automobila.
Algoritam za rješenje ovog porblema je napisan u Pythonu, koji posjeduje bogatu bibilioteku algoritama za Furijeove transformacije.

Dopplerov efekat¶

Doplerov efekat je pojava promjene frekvencije zvuka kojeg registruje slušalac usljed relativnog kretanja u odnosu na izvor zvuka. Jednostavna matematička zavisnost frekvencije izvora $f_0$, brzine izvora $v_0$, brzine prijemnika $v$ i frekvencije koju prijemnik registruje $f$ data je u obliku:
$$f=\dfrac{c \pm v}{c \mp v_0}f_0. $$
U ovom primjeru razmatrat ćemo jednostavan slučaj kretanja izvora zvuka po pravoj liniji stalnom brzinom (automobila ili voza) dok je posmatrač nepokretan, tako da prethodan formula ima oblik:
$$f_1=\dfrac{c}{c - v_0}f_0 $$ kada se izvor približava nepokretnom posmatraču, odnosno:
$$f_2=\dfrac{c}{c + v_0}f_0 $$ kada se udaljava.
Mjerenjem frekvencija $f_1$ i $f_2$ moguće je odrediti brzinu kretanja izvora zvuka:
$$ v = \frac{f_1 −f_2}{f_1+f_2}c. $$

Mjerenje frekvencije¶

Eksperimentalni podaci potrebni da se odredi brzina kretanja izvora su dobiveni snimanjem zvuka vozila koje prolazi pored nepokretnog uređaja za snimanje. Osim toga, mjerena je i temperatura vazduha kako bi se što preciznije odredila brzina zvuka u vazduhu $c$ prema relaciji:
$$ c=c_0 \sqrt{1+\frac{t}{273}}$$ gdje je $t$ mjereno u stepenima Celzijusa, a $c_0=331$ m/s brzina zvuka u vazduhu na temperaturi 0°.
Snimljeni zvučni zapis je za kompjutersku memoriju ništa drugo do niz elemenata od kojih svaki element predstavlja amplitudu zvuka u datom trenutku. Kako bi iz takvog vektora "izvukli" frekvenciju primjenit ćemo moćan matematički alat poznat kao Furijeova transformacija.

Furijeova transformacija¶

Kako je već naglašeno, zvučni zapis je u memoriji računara predstavljen u obliku vekotora $\bf{x}$ sa $N$ elemenata, gdje je svaki element vrijednost amplitude zvučnog talasa u trenutku $\Delta t$.
Vektor $\bf {x}$ se može predstaviti kao superpozicija eksponencijalnih funkcija koje odgovaraju različitim frekvencijama. U praksi, svaki zvuk je sastavljen od velikog broja talasa različitih frekvencija. Svaki element vektora $\bf{x}$ se prema Diskretnoj Furijeovoj transformaciji može pisati kao:
$$x_n=\dfrac{1}{N} \sum_{k=0}^{N-1} X_k e^{i2\pi f_k n \Delta t} $$ gdje su $f_k=\frac{k}{N\Delta t}$ frekvencije, a $X_k$ amplitude. Na taj način, modul amplitude $|X_k|^2$ se može smatrati "težinom" frekvencije $f_k$ u vektoru $\bf{x}$.
Da bismo iz našeg zvučnog zapisa, iz niza elemenata $x_n$ dobili vrijednosti odgovarajućih frekvencija, primjenit ćemo inverznu Furijeovu transformaciju:
$$ X_k=\sum_{n=0}^{N-1}x_n e^{-i2\pi f_k n \Delta t} .$$
Vrijedi ponoviti da je vektor $\bf{x}$ funkcija vremena $\bf{x}(t)$, a vektor $\bf{X}$ funkcija frekvencije $\bf{X}(f)$, tako da će rezultat inverzne Furijeove transformacije biti niz elemenata (amplituda) po različitim frekvencijama.
Računanje inverzne Furijeove transformacije zahtjeva znatnu količinu vremena i računarske memorije, zbog čega je algoritam poznat kao brze Furijeove transformacije efikasan način izračunavanja iste. Prelaskom sa Disretne (DFT) na Brzu (FFT) Furijeovu transformaciju ušteda u vremenu je značajna. Naime za DFT je potrebno izvršiti $N^2$ operacija, dok je u slučaju FFT taj broj reduciran na $N lnN$.
U Pythonu postoji paket numpy sa gotovim funkcijama fft i ifft koje ćemo koristiti za računanje FFT.

Rješavanje problema u Pythonu¶

Pretpostavimo da smo uzorak približavanja vozila posmatraču memorisali u vektor uzorak1, a uzorak udaljavanja vozila u uzorak2. Ove dvije varijable su vektori koji u sebi sadrže zavisnost amplitude zvuka od vremena. Prvi korak je izračunati Furijeove transformacije ovih vektora, tj. primjeniti funkciju fft u Pythonu na vektore uzorak1 i uzorak2.
U drugom koraku se izračunava kvadrat modula amplituda, kvadriranjem svakog elementa iz dobivenih transformacija. Konačno, vrijednosti frekvencija se izračunavaju prema formuli $f=\frac{k}{N\Delta t}$.

Vrijednosti kvadrata amplituda su normalizirane na 1, tako da crtanjem grafika P(f), za oba uzorka, dobijamo određene pikove na frekvencijama. Ono što se primjeti da su ovi pikovi pomjereni za drugi uzorak ka manjim vrijednostima. (navesti sliku). To je u suglasnosti sa onim što posmatrač čuje kada pored njega prolazi vozilo sa uključenom sirenom. Brzinu izvora ćemo odrediti biranjem najvećih vrijednosti frekvencija u oba uzorka, tj. $f_1=max(P1)$ i $f_2=max(P2)$. Nakon toga, jednostavnim računom se dolazi do vrijednosti brzine izvora zvuka.

UBACITI SHEMATSKI PRIKAZ ALGORITMA¶

Python kod¶

# IMPORTOVANJE POTREBNIH BIBILIOTEKA
%matplotlib inline
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fft, ifft
import scipy.io.wavfile as wav
from IPython.core.display import HTML, display

# Set common figure parameters:
newparams = {'axes.labelsize': 11, 'axes.linewidth': 1, 'savefig.dpi': 300, 
             'lines.linewidth': 1.0, 'figure.figsize': (8, 3),
             'ytick.labelsize': 10, 'xtick.labelsize': 10,
             'ytick.major.pad': 5, 'xtick.major.pad': 5,
             'legend.fontsize': 10, 'legend.frameon': True, 
             'legend.handlelength': 1.5}
plt.rcParams.update(newparams)

# DEFINISANJE FUNKCIJE KOJA UČITAVA ZVUČNI ZAPIS i VRAĆA
# FREKVNECIJE I ODGOVARAJUĆE KVADRATE AMPLITUDA

def fftwrapper():
    """
    Output
       P1      :   1D vector
       P2      :   1D vector
       f       :   1D vector
    
    Example usage
       (P1,P2,f) = fftwrapper()
    """
    
    N=120000       # Broj uzoraka (tacaka) u svakom snimku
    shift1=220000  # Prvi index za uzorak 1
    shift2=352000  # Prvi index za uzorak 2
    fcutoff=550    # Maksimalna frekvencija u spektru koji se dobije
    
    # Učitati zvučni zapis i iz stereo prebaciti u mono kanalni
    Fs, ystereo = wav.read('Doppler_car2.wav', 'r')
    ymono = (ystereo[:,0] + ystereo[:,1])/2
    ymono = ymono/max(abs(ymono))
    deltat = 1/Fs
    
    sample1 = ymono[shift1:N+shift1-1]
    sample2 = ymono[shift2:N+shift2-1]
    
    # Izvrsiti FFTs
    p1 = fft(sample1)
    p2 = fft(sample2)
    P1 = np.absolute(p1)**2
    P2 = np.absolute(p2)**2
    f = np.linspace(0,N-1,N)/(N*deltat)
    
    # Isjeci vektore na veličine koje su nam od interesa:
    ifcutoff= np.nonzero(abs(f-fcutoff)==min(abs(f-fcutoff)))[0]-1
    ifcutoff=int(ifcutoff)
    f = f[0:ifcutoff]
    P1 = P1[0:ifcutoff]
    P2 = P2[0:ifcutoff]
    
    return (P1, P2, f)

# pozivanje PRETHODNE FUNKCIJE
(P1, P2, f) = fftwrapper()

# cRTANJE GRAFA
plt.plot(f, P1/max(P1), f, P2/max(P2))
plt.xlabel(r"$f$ (Hz)");
plt.ylabel(r"$P/P_{max}$")
plt.legend(['Uzorak1','Uzorak2','Location','NorthWest'], loc=2);
plt.show()

# ODABIR MAKSIMALNIH VRIJEDNOSTI FREKVENCIJA U OBA UZORKA
f1 = f[P1==max(P1)]
f2 = f[P2==max(P2)]

print("f1 = %f, f2 = %f" % (f1, f2))

f1 = 524.422500, f2 = 484.732500

# iZRAČUNAVANJE BRZINE IZVORA ZVUKA
c = 340.29              # Speed of sound [m/s]
v = (f1-f2)/(f1+f2)*c   # Speed of train [m/s]
v = 3.6*v               # Speed of train [km/h]

print("Brzina izvora je %0.2f km/h." % v)

Brzina izvora je 48.18 km/h.

Kompjutaciona Fizika

Primjer Dopplerovog efekta i Furijeove transformacije