Serie de Fourier

# Cargo paquetitos

using DSP, FFTW, Statistics, WAV

function wavread_mono(file)

  x, sr = wavread(file)

  return mean(x; dims=2)[:], sr

end

# Y armo un par de funciones auxiliares

stem(args...; kwargs...) = sticks(args...; 

  																marker=:circle, 

  																leg=false, 

  																kwargs...)

stem!(args...; kwargs...) = sticks!(args...; 

  																marker=:circle, 

  																leg=false, 

  																kwargs...)

zeropolegain(pr) = DSP.ZeroPoleGain(pr)

zeropolegain(z, p, g) = DSP.ZeroPoleGain(z, p, g)

polynomialratio(zpg) = DSP.PolynomialRatio(zpg)

function polynomialratio(b, a)

  n = max(length(a), length(b))

  return DSP.PolynomialRatio(padright(b, n), padright(a, n))

end

getpoles(zpg) = DSP.ZeroPoleGain(zpg).p

getzeros(zpg) = DSP.ZeroPoleGain(zpg).z

getgain(zpg) = DSP.ZeroPoleGain(zpg).k

getnumcoefs(pr) = trimlastzeros!(reverse(DSP.PolynomialRatio(pr).b.coeffs))

getdencoefs(pr) = trimlastzeros!(reverse(DSP.PolynomialRatio(pr).a.coeffs))

function trimlastzeros!(a)

  !iszero(a[end]) && return a

  pop!(a)

  return trimlastzeros!(a)

end

DSP.filt(zpg::DSP.ZeroPoleGain, r...; kwargs...) = filt(polynomialratio(zpg), r...; kwargs...)

function zplane(zs, ps; kwargs...)

	scatter(real.(zs), imag.(zs);

		  marker = (:black, :circle), label="Cero", kwargs...)

	scatter!( real.(ps), imag.(ps);

	  	marker = (:red, :xcross), label="Polo", kwargs...)  

  ts = range(0,stop=2pi;length=100)

  plot!(cos.(ts), sin.(ts); aspect_ratio = 1, kwargs...)

end

function zplane(zs, ps; kwargs...)

	scatter(real.(zs), imag.(zs);

		  marker = (:black, :circle), label="Cero", kwargs...)

	scatter!( real.(ps), imag.(ps);

	  	marker = (:red, :xcross), label="Polo", kwargs...)  

  ts = range(0,stop=2pi;length=100)

  plot!(cos.(ts), sin.(ts); aspect_ratio = 1, kwargs...)

end

zplane(pr::DSP.PolynomialRatio; kwargs...) = zplane(DSP.ZeroPoleGain(pr); kwargs...)

# Delta

d(n) = n == 0 ? 1. : 0. 

# Escalón

u(n) = n >= 0 ? 1. : 0. 

using Plots

Plots.default(:legend, false)

# 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

# Espectrograma

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

specplot(x::AbstractMatrix; kwargs...) = @error "You are entering a Matrix (2D Array). I need a Vector (1D Array)."

function specplot(x::AbstractVector; 

      fs=1, 

      onesided=false, 

      xaxis="Tiempo (s)", 

      yaxis="Frecuencia (Hz)",

      window=hamming(div(length(x), 16)),

      overlap=0.5,

      nfft=length(window),

      kws...)

    window isa Integer && (window = rect(window))

    overlap isa AbstractFloat && (overlap = round(Int, length(window) * overlap))

    mat = stft(x; overlap=overlap, window=window, nfft=nfft)

    fmax = fs

    if onesided

      mat = mat[1:div(size(mat, 1) + 2, 2), :]

      fmax = fs/2

    end

  toffset = length(window) / 2sr

  times = range(toffset; length=size(mat, 2), stop=length(x)/fs - toffset)

  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

[ x(0), x(1), x(2), ..., x(N - 1) ]

[ a(0), a(1), a(2), ..., a(N - 1) ] 

sum(2i for i in 1:3)

sum( 2 .* 1:3 )

function coefs_serie(x)

   N = length(x)

   a = zeros(N)

   for k in 0:N-1

       a[k + 1] = 1/N * sum( x[n + 1] * exp(-im* 2π / N * k * n) for n in 0:N-1)

   end

   return a

end

function serie_synth(a)

   N = length(a)

   x = zeros(ComplexF64, N)

   for n in 0:N-1

       x[n + 1] = sum( a[k + 1] * exp(im* 2π / N * k * n) for k in 0:N-1)

   end

   return x

end

coefs_serie([ 1, 0, 0, 0])

serie_synth([0.25, 0.25, 0.25, 0.25]) |> real

x = [ 1, 0, 0, 0, 0, 0, 0, 0 ]

a = [ 1, 1, 1, 1, 1, 1, 1, 1 ] ./ length(x)