Espectrograma
Código auxiliar
# Cargo paquetitos
using DSP, FFTW
# Y armo un par de funciones auxiliares
stem(args...; kwargs...) = sticks(args...;
marker=:circle,
leg=false,
kwargs...)
# Delta
d(n) = n == 0 ? 1. : 0.
# Escalón
u(n) = n >= 0 ? 1. : 0.
using Plots # esto está demorando como un minuto :S
32.5s
# Pad vector with zeros on the right until its length is `n`
padright(x, n) = copyto!(zeros(eltype(x), n), x)
"""
Función módulo pero con offset (opcional)
Manda a `t` al intervalo [from, from+length)
sumándole o restándole algún múltiplo de `len`
"""
cshift(t, len, from=0) = mod(t - from, len) + from
1.1s
cshift
Armo función custom
using IterTools
function stft(x; overlap, window, nfft, rest...)
nwin = length(window)
overlap < nwin
res = [ fft(padright(xseg .* window, nfft))
for xseg in partition(x, nwin, nwin - overlap)]
return [ res[i][j] for j in 1:nfft, i in eachindex(res)]
end
function specgram(x;
overlap=0.5,
window=hamming(div(length(x), 16)),
nfft=length(window),
rest...)
window isa Integer && (window = rect(window))
overlap isa AbstractFloat && (overlap = round(Int, length(window) * overlap))
return stft(x; overlap=overlap, window=window, nfft=nfft)
end
function specplot(x::AbstractVector;
fs=1,
onesided=false,
xaxis="Tiempo (s)",
yaxis="Frecuencia (Hz)",
kws...)
mat = specgram(x; kws...)
fmax = fs
if onesided
mat = mat[1:div(size(mat, 1) + 2, 2), :]
fmax = fs/2
end
times = range(0; length=size(mat, 2), stop=length(x)/fs) # aprox
freqs = range(0; length=size(mat, 1), stop=fmax)
# Reubico las frecuencias negativas arriba de todo
if !onesided
freqs = cshift.(freqs, fs, -fs/2)
ord = sortperm(freqs)
mat = mat[ord, :]
freqs = freqs[ord]
end
return heatmap(times, freqs, log.(abs.(mat) .+ eps());
xaxis=xaxis, yaxis=yaxis,
seriescolor=:bluesreds, legend=true, kws...)
return times, freqs, mat
end
function specplot(x :: AbstractVector{<:AbstractFloat}; kws...)
return specplot(convert.(Complex, x); onesided=true, kws...)
end
0.2s
specplot (generic function with 2 methods)
let
x = cos.(2π .* range(0; length=1024, step=0.1))
x += 1e-2 .* rand(size(x)...)
specplot(x;
overlap=0.9,
window=rect(512), fs=4)
end
1.4s
let
x = cos.(2π .* range(0; length=1024, step=0.1))
x += 1e-2 .* rand(size(x)...)
specplot(x;
overlap=0.9,
window=hamming(256), onesided=false, fs=4)
end
2.4s
let
fs = 10e3
N = 1e5
amp = 2 * sqrt(2)
noise_power = 0.01 * fs / 2
time = (0:N-1) ./ fs
mod = 1000 .* cos.(2pi * 0.25 .* time) # moduladora
carrier = amp .* sin.(2pi * 3e3 .* time .+ mod)
noise = randn(size(time)) .* sqrt(noise_power)
noise .*= exp.(-time ./ 5)
x = carrier .+ noise
specplot(x;
overlap=0.9,
window=tukey(256, 0.5),
fs=fs)
end
4.1s
Audio
using WAV
x, sr = wavread(fafafa.wav)
specplot(x; fs=sr, window=hamming(1024), ylims=(0, 10e3))
6.3s
specplot(x; fs=sr, window=hamming(4086), overlap=0.95, ylims=(0, 1e3))
20.0s
La del paquete DSP
s = let
x = cos.(2π .* range(0; length=1024, step=0.1))
x += 1e-2 .* rand(size(x)...)
spectrogram(x, 256, 220; fs=4, window=hamming)
end
;
1.9s
function Plots.plot(spec::DSP.Periodograms.Spectrogram; kwargs...)
time, freq, pow = spec.time, spec.freq, spec.power
# Reubico las frecuencias negativas al comienzo
nneg = findfirst(x -> x<0, spec.freq)
if nneg !== nothing
npos = nneg - 1
freq = [ spec.freq[nneg:end]; spec.freq[1:npos] ]
pow = [ spec.power[nneg:end, :];
spec.power[1:npos , :] ]
end
return heatmap(time, freq, pow .+ eps() .|> log; seriescolor=:bluesreds, kwargs...)
end
0.2s
plot(s)
1.9s
let
x = cos.(2π .* range(0; length=1024, step=0.1))
x += 1e-2 .* rand(size(x)...)
plot(spectrogram(x, 256, 220; fs=4, window=hamming, onesided=false))
end
1.2s