Espectrograma

# 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

# 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

using IterTools

function stft(x; overlap, window, nfft, rest...)

  nwin = length(window)

  @assert 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

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

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

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

using WAV

x, sr = wavread(
fafafa.wav)

specplot(x; fs=sr, window=hamming(1024), ylims=(0, 10e3))

specplot(x; fs=sr, window=hamming(4086), overlap=0.95, ylims=(0, 1e3))

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

;

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

plot(s)

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