Abychom mohli pracovat se zvukovými daty, musíme je nejprve importovat. K tomu účelu využijeme balíček WAV.jl. Pokud nemáte tento balíček nainstalován, stačí v interaktivní příkazové řádce Julia (REPL) vydat příkaz Pkg.add("WAV").
In [1]:
# balíček pro práci s WAV soubory
using WAV
Začněme s jednoduchou ukázkou. Příkaz wavread načte ze zadaného WAV souboru vzorky, vzorkovací frekvenci a počet bitů použitých na kódování vzorků. Ve výchozím nastavení jsou hodnoty vzorků v rozmezí od $-1.0$ do $1.0$.
In [2]:
# načtení dat
y, fs, nbits = wavread("violin.wav")
# info
println("Frekvence: $fs")
println("Bitů: $nbits")
println("Prvních deset vzorků:")
println(y[1:10])
Záznam si můžeme přehrát, stačí zavolat wavplay s vzorky a frekvencí.
In [21]:
wavplay(y,fs)
Pojďme se nejprve podrobněji podívat na data v souboru, který jsme v předchozí části načetli. K tomu použijeme velmi jednoduchý balíček na generování grafů Winston.
In [3]:
# balíček pro generování grafů
using Winston
Vykresleme celý záznam.
In [15]:
# vzorkovací frekvence je fs, časový krok mezi vzorky
# proto je 1/fs sekund.
p = FramedPlot(
title="violin.wav",
xlabel="Čas [s]",
ylabel="Amplituda"
)
setattr(p.frame, draw_grid=true)
add(p, Curve(Float64[k*1/fs for k = 1:length(y)], y, color="blue", linewidth=1))
savefig(p,"violin.pdf",width=800,height=600) # uložíme obrázek pro další zpracování
p
Out[15]:
Nenechte se zmást. Předchozí graf sice vypadá jako skvrna, ale skutečně jde o záznam jistých oscilací. Podívejme se na kratší časový interval.
In [20]:
indmin = 10000
indmax = 11000
p = FramedPlot(
title="violin.wav",
xlabel="Čas [s]",
ylabel="Amplituda"
)
setattr(p.frame, draw_grid=true)
add(p, Curve(Float64[k*1/fs for k = indmin:indmax], y[indmin:indmax], color="blue", linewidth=4))
Out[20]:
Abych prozkoumali jaké frekvence jsou v signálu obsaženy, použijeme diskrétní Fourieorovu transformaci (DFT). Pro připomenutí uvádíme definiční vztah z přednášky (bez normalizačního faktoru):
$$ \hat x_k = \sum_{j=0}^{N-1} x_k e^{-2i\pi kj/N}\,, \quad k = 0,1,\ldots,N-1\,. $$Vstupnímu vektoru $(x_0,x_1,\ldots,x_{N-1}) \in \mathbb{C}^N$ délky $N$ je tímto předpisem přiřazen vektor $(\hat x_0, \hat x_1,\ldots,\hat x_{N-1})\in\mathbb{C}^N$. Poznamenejme, že jsme vynechali normalizační člen $\frac{1}{\sqrt{N}}$, na výsledky v tomto notebooku nemá žádný vliv. Naivní implementací DFT by tedy bylo prosté přepsání tohoto vztahu:
In [21]:
# naivní implementace DFT
function naiveDFT(data::Vector{Complex64})
# délka vektoru
n = length(data)
# alokace paměti pro uložení výsledku
ft = Array(Complex64, (n,1))
# napočtení složek výsledného vektoru
for k in 0:n-1
ft[k+1] = sum(j->data[j+1]*exp(-2im*pi*k*j/n), 0:n-1)
end
return ft
end
# metoda pro reálný vektor
naiveDFT(data::Vector{Float64}) = naiveDFT(convert(Vector{Complex64}, data))
Out[21]:
Na přednášce jsme si ukázali Cooley-Tukey algoritmus, který počítá DFT na vektorech, jejichž délka je mocninou dvou. Protože jeho složitost (vzhledem k délce vstupu) je $\mathcal{O}(N\log N)$, na rozdíl od složitosti $\mathcal{O}(N^2)$ naivní implementace, mluvíme o něm jako o rychlé Fourierově transformaci (FFT).
In [22]:
# FFT, rekurzivní Cooley-Tukey algoritmus
function myFFT(data::Vector{Complex64})
# délka vektoru a poloviční délka vektoru
n, halfn = length(data), div(length(data), 2)
if n == 2
# DFT dvousložkového vektoru
return Complex64[data[1] + data[2], data[1] + exp(-pi*1im)*data[2]]
else
# "decimace", N.B. v původním vzorci indexujeme
# od 0, Julia indexuje od 1.
even = myFFT(data[1:2:n-1])
odd = myFFT(data[2:2:n])
return Complex64[ even[mod(k,halfn)+1] + exp(-2im*pi*k/n)*odd[mod(k,halfn)+1] for k = 0:n-1]
end
end
# metoda pro reálný vektor
myFFT(data::Vector{Float64}) = myFFT(convert(Vector{Complex64}, data))
Out[22]:
Julia má samozřejmě FFT k dispozici i samostatně pomocí funkce fft. Dle dostupných informací (zdrojový kód julia/base/fftw.jl) Julia využívá knihovnu FFTW. Tato knihovna nabízí i další zobecnění Cooley-Tukey algoritmu i na vektory libovolné délky.
Porovnejme rychlosti (čas) a výsledky těchto tří funkcí (naiveDFT, myFFT a fft) na náhodných datech. Upozorňujeme, že je vždy potřeba funkce invokovat aspoň dvakrát, při prvním volání se kompilují (Julia využívá JIT).
In [24]:
# náhodný vektor dély 1024
data = rand(1024);
println("naiveDFT:")
@time naivedata = naiveDFT(data);
println("myFFT:")
@time mydata = myFFT(data);
println("fft:")
@time jdata = fft(data);
Zkontrolujme, že v rámci chyb při aritmetických operacích v konečné přesnosti jsme získali stejné výsledky.
In [11]:
println("Maximální absolutní rozdíly ve složkách výstupu:")
println("Naivní vs. naše FFT: ", maximum(abs(naivedata - mydata)))
println("Naše FFT vs. FFTW ", maximum(abs(mydata - jdata)))
In [33]:
kstart = 3
klen = 10
telapsed = Array(Float64,(klen,4))
for k in kstart:kstart+klen-1
data = rand(2^k)
telapsed[k-kstart+1,:] = [2^k,
mean([@elapsed(naiveDFT(data)) for j in 1:10]),
mean([@elapsed(myFFT(data)) for j in 1:10]),
mean([@elapsed(fft(data)) for j in 1:10])
]
end
println("Done!")
Následující tabulka porovnává průměrné doby (v sekundách) běhu tří výše zmíněných funkcí na vektorech délky n.
In [34]:
using DataFrames
DataFrame(n=telapsed[:,1], dft=telapsed[:,2], myfft=telapsed[:,3], fftw=telapsed[:,4])
Out[34]:
Vraťme se k významu samotné diskrétní Fourierovy transformace. Stručně lze říci, že hodnota její složky, $\hat x_n$ udává, jakou mírou do vstupního signálu $(x_k)_{k=0}^{N-1}$ přispívá harmonická oscilace s frekvencí $n/T$, kde $T$ je celková délka vstupního signálu $(x_k)_{k=0}^{N-1}$ (v sekundách). Tedy $T = N/f_s$, kde $f_s$ je vzorkovací frekvence. Frekvence je pak dána v Herzech. Například pro zpracování záznamu řeši se často volí délka okna 25 ms a používá se překryv 10 ms.
V této sekci sestrojíme tak zvaný spektrogram zvukového záznamu. Jedná se o poměrně přímočarý proces. Vstupní signál rozdělíme na kratší úseky, tzv. "okna", která se překrývají. Dále vypočteme DFT každého okna. Do grafu pak na vodorovnou osu vynášíme čas a na svislou osu frekvenci. Bod s danou souřadnicí pak vybarvíme podle toho, jaká je hodnota Fourierovy transformace (typicky dle kvadrátu absolutní hodnoty, intenzity).
Před provedením DFT se záznam v okně ještě často upravuje přenásobením tzv. okénkovou funkcí. My budeme nejčastěji používat Hannovu funkci. Vliv této úpravy signálu na výsledný spektrogram si ukážeme o několik řádek níže.
In [25]:
# Hannova okénková funkce
function hanning(n::Integer)
[0.5*(1 - cos(2*pi*k/(n-1))) for k=0:(n-1)]
end
# Pravoúhla okénková funkce
function rectangular(n::Integer)
ones(Float64,n)
end;
Graf pro lepší představu.
In [27]:
p = FramedPlot(
title="Okénkové funkce",
xlabel="Délka záznamu",
ylabel="Amplituda"
)
setattr(p.frame, draw_grid=true)
fighann = Curve(1:1024, hanning(1024), color="red", linewidth=7)
setattr(fighann, label="Hann")
figrect = Curve(1:1024, rectangular(1024), color="blue", linewidth=7)
setattr(figrect, label="Pravoúhlá")
leg = Legend(.8, .9, {fighann, figrect})
add(p, figrect, fighann, leg)
savefig(p,"window.pdf",width=800,height=600)
p
Out[27]:
Vykreslování provedeme pomocí následující funkce.
In [37]:
function spec(data::Vector{Float64}, # vstupní signál
freq::Real; # vzorkovací frekvence
wlen::Float64=0.025, overlap::Float64=0.01, # nepovinná délka okna, překryv (v sekundách)
maxfreq=pi*freq, logscale=true, wfunc=hanning, # ořez frekvence, měřítko intenzity, okénková funkce
mytitle="Spektrogram", # titulek grafu
filename="" # ukládání do souboru
)
# délka záznamu (počet vzorků
n = length(data)
# délka okénka (počet vzorků)
wn = round(Int64, wlen*freq)
# základní frekvence
bfreq = freq/wn
# délka překryvu (počet vzorků)
overlapn = round(Int64, overlap*freq)
# předpočítání Hanningovy funkce
wdat = wfunc(wn)
# pole kam budeme zapisovat výsledky
spectrum = Array(Complex64, (wn,0))
# průběžný index
k = 0
while k + wn < n
# aplikace Hanningovy funkce
spectrum = hcat(spectrum, fft(data[k+1:k+wn] .* wdat))
k += wn - overlapn
end
if logscale
imgdata = map(x -> log(1+abs(x)^2), spectrum)
else
imgdata = map(x -> abs(x)^2, spectrum)
end
p = imagesc((0,n/freq), (maxfreq,0), imgdata[1:floor(Int64,maxfreq/bfreq),:])
setattr(p,title=mytitle)
setattr(p,xlabel="Čas [s]")
setattr(p,ylabel="Frekvence [Hz]")
if !isempty(filename)
savefig(filename,width=800,height=600)
end
p
end
Out[37]:
Spektrogram záznamu houslí, okénková funkce Hannova, okénko má délku 25ms a překryv 10ms.
In [39]:
spec(y[:,1], fs, maxfreq=5000, mytitle="Spektrogram: violin.wav",filename="spec_violin.pdf")
Out[39]:
Spektrogram záznamu houslí, okénková funkce pravoúhlá, okénko má délku 25ms a překryv 10ms.
In [20]:
spec(y[:,1], fs, maxfreq=5000, mytitle="Spektrogram: violin.wav", wfunc=rectangular)
Out[20]:
Pozorujeme, že okénková funkce má vliv na rozlišení, rozmazání, spektrogramu.
Nyní se na problém podívejme z druhé strany. Naopak vygenerujme signál, přehrajme ho a podívejme se jak dopadne jeho spektrogram. Zadefinujme si nejprve frekvence standardního ladění šestistrunné kytary.
In [47]:
gtuning = {
:e4 => 329.63,
:b3 => 246.94,
:g3 => 196.00,
:d3 => 146.83,
:a2 => 110.00,
:e2 => 82.41
}
Out[47]:
Následující funkce vygeneruje harmonickou oscilaci s předepsanou frekvencí freq, vzorkovací frekvencí sfreq a délkou signálu duration. Frekvence je opět udávána v Herzech a čas v sekundách.
In [48]:
function sinwave(freq::Float64, sfreq::Real, duration::Float64)
n = floor(Int64,duration*sfreq)
return Float64[ sin(2*pi*freq * k/sfreq) for k in 1:n]
end
Out[48]:
Například:
In [49]:
fs = 44100
y = sinwave(gtuning[:e4], fs, 3.0);
Graf tohoto záznamu dává skutečně očekávanou sinusovku.
In [50]:
p = FramedPlot(
title="E4",
xlabel="Čas [s]",
ylabel="Amplituda"
)
setattr(p.frame, draw_grid=true)
add(p, Curve(Float64[k*1/fs for k = 1:500], y[1:500], color="blue", linewidth=7))
Out[50]:
Můžeme si tento signál i přehrát!
In [51]:
wavplay(y,fs)
Ve spektrogramu očekáváme jednu frekvenci.
In [53]:
spec(y, fs, maxfreq=400, wlen=0.3, overlap=0.01, filename="spec_sin329.pdf")
Out[53]:
Všechny tóny od nejvyššího k nejnižšímu.
In [32]:
wave = vcat(
sinwave(gtuning[:e4], fs, 2.0),
sinwave(gtuning[:b3], fs, 2.0),
sinwave(gtuning[:g3], fs, 2.0),
sinwave(gtuning[:d3], fs, 2.0),
sinwave(gtuning[:a2], fs, 2.0),
sinwave(gtuning[:e2], fs, 2.0)
);
wavplay(wave, fs)
Spektrogram tohoto záznamu:
In [63]:
spec(wave, fs, maxfreq=500, wlen=0.045, overlap=0.01)
Out[63]:
Následující pomocná funkce pouze sestaví superpozici více signálu a výsledek normalizuje na jedničku (ve smyslu absolutní hodnoty vzorku).
In [33]:
function normalizesum(args...)
s = +(args...)
m = maximum(abs(s))
return 1/m .* s
end
Out[33]:
Sestavme aproximaci 'zubového' signálu, jehož fourierovu řadu jsme počítali na přednášce. Za tím účelem použijeme šest harmonických oscilací.
In [35]:
saw = normalizesum(
2*sinwave(gtuning[:e4], fs, 3.0),
-1*sinwave(2*gtuning[:e4], fs, 3.0),
2/3 * sinwave(3*gtuning[:e4], fs, 3.0),
-2/4 * sinwave(4*gtuning[:e4], fs, 3.0),
+2/5 * sinwave(5*gtuning[:e4], fs, 3.0),
-2/6 * sinwave(6*gtuning[:e4], fs, 3.0),
);
Graf (z letáku).
In [36]:
p = FramedPlot(
title="Pilový signál",
xlabel="Čas [s]",
ylabel="Amplituda"
)
setattr(p.frame, draw_grid=true)
add(p, Curve(Float64[k*1/fs for k = 1:500], saw[1:500], color="blue", linewidth=7))
Out[36]:
Jak tento signál zní?
In [37]:
wavplay(saw, fs)
Spektrogram:
In [38]:
spec(saw, fs, maxfreq=2200)
Out[38]:
Pro zajímavost se podívejme na několik zvukových záznamů elektrické kytary. Nejprve prázdné struny.
In [40]:
y, fs, nbits = wavread("guitar.wav");
println("Frekvence: $fs")
println("Bitů: $nbits")
Přehrání záznamu.
In [41]:
wavplay(y,fs)
Spektrogram.
In [44]:
spec(y[:,1], fs, maxfreq=2000, mytitle="Kytara, prázdné struny",wlen=0.06, overlap=0.025,filename="spec_guitar.pdf")
Out[44]:
Podívejme se na menší výřez.
In [72]:
spec(y[:,1], fs, maxfreq=1000, mytitle="Kytara, prázdné struny", wlen=0.05, overlap=0.025)
Out[72]:
Nyní akord G.
In [45]:
y, fs, nbits = wavread("guitar_g.wav");
wavplay(y,fs)
spec(y[:,1], fs, maxfreq=1500, mytitle="Kytara, akord G", wlen=0.1, overlap=0.025)
Out[45]:
Akord G a A s kvákadlem.
In [46]:
y, fs, nbits = wavread("guitar_gaw.wav");
wavplay(y,fs)
spec(y[:,1], fs, maxfreq=1500, mytitle="Kytara, akord G, A a kvákadlo", wlen=0.06, overlap=0.03)
Out[46]:
Riff z 'For whom the bell tolls' od Metallica.
In [47]:
y, fs, nbits = wavread("guitar_m.wav");
wavplay(y,fs)
spec(y[:,1], fs, maxfreq=1000, mytitle="Rozostřena kytara", wlen=0.06, overlap=0.03)
Out[47]: