Ecuaciones en diferencias
Código auxiliar
# 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)
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...) = "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
Algunas funciones para trabajar con sistemas definidos por ecuaciones en diferencias
Los sistemas LTI definidos por ecuaciones en diferencias tiene una transferencia H(z) racional. Quedan totalmente definidos por los coeficientes del polinomino del numerador y del denominador. Alternativamente, por la lista de polos, ceros, y un valor de ganancia.
Los sistemas FIR (respuesta al impulso finita) expresan la salida como una convolución que ya es una ecuación en diferencias, de orden 0 (sólo aparece y(n)). El denominador de H(z) es 1. No tiene polos (salvo quizás en 0 o infinito).
Si queremos ver la salida de un sistema FIR, usamos la función `conv`.
Los sistemas IIR no se pueden implementar con la función `conv` ya que h(n) es infinita. Pero, los que vienen de ecuaciones en diferencias, sí. La PC implementa la versión inicialmente en reposo, iterando hasta cierto n.
La función que hace esto se llama `filt` (en MATLAB y otros, `filter`).
#=
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))
Se puede hacer un diagramita de polos y ceros a partir del vector de polos y de ceros con la función zplane
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)))