Ecuaciones en diferencias

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

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(zpg; kwargs...)

	zs = getzeros(zpg)

	ps = getpoles(zpg)

	isempty(zs) && isempty(ps) && (append!(ps, 0))

	return zplane(zs, ps; 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

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

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

      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

#=

Si se quiere tener la salida inicialmente en reposo para este sistema

a[1] * y(n) + a[2] * y(n-1) +... = b[1] * x(n) + b[2] * x(n-1) + ...

frente a una entrada definida pro el vector x (x(n) = 0 para n<0, x(0)=x[1], etc)

uno usa

filt(b, a, x)

La salida tendrá tantos elementos como x

=#

# Caso FIR

filt([0, 1], [1], 1:10)

# Caso IIR

a = [1, 3.4] # Sin pérdida de generalidad, a(1) puede ser siempre 1

b = [3.2, 1.3, 6.7]

filt(b, a, 1:3)

# Si querés más valores, meté una entrada más larga. Si terminó la entrada, agregar ceros

filt(b, a, padright(1:3, 5))

zs = [-1.34, 0]

ps = 0.8 .* exp.(im .* range(0; step=2pi/10, length=10))

zplane(zs, ps)

# Podemos armar un objeto que representa un filtro IIR a partir de los polos, ceros y ganancia

gain = 1

zpg = zeropolegain(zs, ps, gain);

# y recuperar los zeros

getzeros(zpg)

# o polos

getpoles(zpg)

# o ganancia

getgain(zpg)

# o tb podemos armar un objeto en términos de los vectores de coeficientes a y b definidos como para filt

pr = polynomialratio(b, a);

# y recuperar numerador

getnumcoefs(pr)

# o denominador

getdencoefs(pr)

# Cualquiera de estas se pueden usar con zplane

zplane(zpg)

zplane(pr)

# y cualquiera permite extraer polos, ceros

getzeros(pr), getpoles(pr)

# o coeficientes del polinomio numerador y denominador

getnumcoefs(zpg), getdencoefs(zpg)

# filt acepta estos objetos también

filt(pr, 1:3), filt(b, a, 1:3), filt(zpg, 1:3)

# También podés usar estos objetos en la función freqz, que

# devuelve la respuesta en frecuencia del sistema estable

# evaluada en cierto omegón (ojo que el estable puede no ser el

# inicialmente en reposo)

w = pi/4

freqz(pr, w)

# No hace más que la cuentita de evaluar H(z=exp(im*w))

cuentita(w) = sum(b .* exp(-im*w) .^ (0 : length(b) - 1)) / 

  sum(a .* exp(-im*w) .^ (0 : length(a) - 1))

cuentita(w)

fs = range(0; stop=pi, length=300)

plot(fs, abs.(freqz(pr, fs)))