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¶
In [16]:
# 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)
In [17]:
# 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)
In [18]:
# pozivanje PRETHODNE FUNKCIJE
(P1, P2, f) = fftwrapper()
In [19]:
# 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()
In [20]:
# ODABIR MAKSIMALNIH VRIJEDNOSTI FREKVENCIJA U OBA UZORKA
f1 = f[P1==max(P1)]
f2 = f[P2==max(P2)]
print("f1 = %f, f2 = %f" % (f1, f2))
In [21]:
# 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)
0 Primjedbe