Tarea 27-07-2020

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...)
# 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
47.2s
Julia
→
specplot (generic function with 3 methods)

Ej 2)

# Secuencia arbitraria de longitud 10
n = 10
x = rand(10) .* 2 .- 1
stem(x)
47.3s
Julia
→
# Armo una función para submuestrear
function downsample(x, n)
  out = zeros(eltype(x), length(x) * n)
  out[1:n:end] .= x
  return out
end
0.7s
Julia
→
downsample (generic function with 1 method)
downsample(x, 3)
1.4s
Julia
→
30-element Array{Float64,1}: 0.301791 0.0 0.0 -0.676953 0.0 0.0 -0.601105 0.0 0.0 -0.248725 ⋮ -0.336294 0.0 0.0 0.698397 0.0 0.0 -0.728044 0.0 0.0
stem(abs.(fft(x)))
0.5s
Julia
→
stem(abs.(fft(downsample(x, 2))))
0.5s
Julia
→
stem(abs.(fft(downsample(x, 3))))
0.5s
Julia
→

Se ve que queda igual. Sólo estamos graficando más períodos.

Teóricamente

formula not implemented

Tomamos inline_formula not implemented y sumamos para los m enteros nomás, ya que para los otros términos la señal vale cero

formula not implemented

11) a) b) y d)

a)

Es una ventana rectangular de longitud 20 que arranca en 0. Da una sinc periódica.

Tomando la DFT de 20 puntos o más, consigo muestras de la transformada entre 0 y 2pi.

xa = padright(ones(20), 200)
xaf = fft(xa)
plot(range(-pi; stop=pi, length=length(xaf)), fftshift(abs.(xaf)))
0.7s
Julia
→

b)

xb = padright([1., -1.], 200)
xbf = fft(xb)
plot(range(-pi; stop=pi, length=length(xbf)), fftshift(abs.(xbf)))
0.6s
Julia
→

d)

Están entrando armónicas puras, que son autofunciones de sistemas LTI. Sale exactamente lo mismo multiplicado por su autovalor, que es la transferencia -- transformada de h(n) -- en la frecuencia. Necesito entonces saber la transferencia en 2pi*3/100 y 2pi*1.5/100.

Si muestreo la transformada en una grilla de 200 puntos entre 0 y 2pi, el 7mo cae en 2pi*6/200 = 2pi*3/100, y el 4to en 2pi*3/200 = 2pi*1.5/100

xaf = fft(padright(ones(20), 200))
avalA1 = xaf[7]
avalA2 = xaf[4]
xbf = fft(padright([1., -1.], 200))
avalB1 = xbf[7]
avalB2 = xbf[4]
;
0.2s
Julia
→

13) a) c)

La segunda porque tiene cambios más abruptos

Voy a obtener 8 muestras de las TF en ambos casos. En módulo, van a coincidir, ya que son iguales a menos de un desplazamiento circular que sólo cambia la fase. Estas 8 muestras no alcanzan para ver que efectivamente la segunda tiene más componentes de alta frecuencia, pero tomando más, esperaría que se viera.

xl = [0, 1, 2, 3, 4, 3, 2, 1]
xr = [4, 3, 2, 1, 0, 1, 2, 3]
stem(abs.(fft(xl)))
0.7s
Julia
→
stem(abs.(fft(xr)))
0.6s
Julia
→
# Igual, como esperábamos
# Ahora, si tomo más muestras
plot(abs.(fft(padright(xl, 100))))
0.5s
Julia
→
plot(abs.(fft(padright(xr, 100))))
0.6s
Julia
→

Recordemos que las altas frecuencias son las que están en el medio acá (cerca de frecuencia pi).

Si graficáramos este último como corresponde...

plot(
  range(-pi; stop=pi, length=100),
  fftshift(abs.(fft(padright(xr, 100)))))
0.5s
Julia
→

Runtimes (1)