header ads

Primjer Dopplerovog efekta i Furijeove transformacije

Kvant svibnja 13, 2018


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

f1 = 524.422500, f2 = 484.732500
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)

Brzina izvora je 48.18 km/h.
Tags:

Objavi komentar

0 Primjedbe