Usar la PC para ver espectros de señales

Código auxiliar

# 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. 
"""
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 Plots
Plots.default(:legend, false)
30.1s
Julia
false

Señales discretas

inline_formula not implemented para inline_formula not implemented

x1(n) = cos(2*π*n/10) * (u(n) - u(n-100));
0.2s
Julia

La graficamos

let
  ns = -5:105
  xs = x1.(ns)
  stem(ns, xs)
end
38.9s
Julia

Queremos el espectro de esto, aunque sea aproximado, sin cuentas analíticas.

n_tot = 500 # la cantidad de muestras que voy a tomar. ¡Incluye todo el soporte!
ns    = range(0; length=n_tot) # un vector de índices `n`
x1s = x1.(ns) # un vector con las primeras `n_tot` muestras
# Calculo su DFT
x1s_dft = fft(x1s)
# Armo un vector de las frecuencias donde caen las muestras
fs = range(0; step=2π/n_tot, length=n_tot)
# Grafico el módulo, uniendo las muestras con una recta
plot(fs, abs.(x1s_dft))
30.1s
Julia
# Graficamos pero en el intervalo central 
let
 # Mandamos la segunda mitad al inicio (fftshift)
 # ...y reacomodamos las frecuencias en el intervalo [-pi pi)
  fs_aux = cshift.(fftshift(fs), 2π, from=-π)
  
  plot(fs_aux, fftshift(abs.(x1s_dft)) )
end
1.9s
Julia

Señal continua

Esto es lo mismo que antes, luego de un muestreo. Imaginemos una señal triangular que arranca en cero

x2(t) = tr(t/2 - 0.5)
# tr(t) es una funcióen triángulo centrado entre -0.5 y 0.5
function tr(t)
  res = 1 - 2 * abs(t)
  return -0.5 < t < 0.5 ? res : zero(res)
end
;
0.1s
Julia
plot(x2, -1, 5)
2.7s
Julia

La muestreamos con alguna frecuencia de muestreo `sr` lo suficientemente grande como para que la energía de la señal esté (casi) toda contenida en inline_formula not implemented.

sr = 100
x2d(n) = x2(n / sr)
;
0.1s
Julia

y ahora es un problema como el de antes

  • Tomamos un número de muestras que contenga todo el soporte - i.e. mayor a inline_formula not implemented

  • DFT, y el resultado son muestras del espectro de la discreta

  • Lo grafico centrado

Ahora, como muestreamos "con Nyquist", podemos interpretar a las frecuencias discretas como frecuencias continuas inline_formula not implemented

Hagamos todo directo, sin pasar por la señal discreta.

# Los tiempos donde voy a muestrear para armar el vector
ts = range(0; stop = 3, step=1/sr)
x2s = x2.(ts)
x2s_fft = fft(x2s)
plot(
  range(-sr/2; stop=sr/2, length=length(ts)),  # esto es aprox
  fftshift(abs.(x2s_fft)))
0.9s
Julia

Podríamos armar una funcioncita que grafica siguiendo estos pasos

function plotspec(x, sr, tmax; kw...)
  ts = range(0; stop = tmax, step=1/sr)
  xs = x.(ts)
  xs_fft = fft(xs)
  plot(
    range(-sr/2; stop=sr/2, length=length(ts)),  # esto es aprox
    fftshift(abs.(xs_fft)); 
    xlabel = "f (Hz)",
    kw...)  
end
0.6s
Julia
plotspec (generic function with 1 method)
plotspec(x2, sr, 3; title="|X2(f)|")
1.4s
Julia
0.1s
Julia
Runtimes (1)