Diseño de filtros digitales FIR

Método de Ventaneo

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...)
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
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
34.4s
Julia
→

El problema del filtro ideal

Charles Ventaneo necesitaba armar un filtro pasabajos FIR de frecuencia de corte inline_formula not implemented. Quería conseguir un vector de números h con el que pudiera filtrar cualquier señal x haciendo conv(h, x). 

H_ideal(Ω) = u(Ω + π / 4) - u(Ω - π / 4)
plot(H_ideal, -π, π)
3.2s
Julia
→

La respuesta al impulso sabemos que es inline_formula not implemented

h_ideal(n) = 1/4 * sinc(n / 4)
let ns = -50:50
  stem(ns, h_ideal.(ns))
end
0.5s
Julia
→

Charles era buen alumno de SyS, y notó rápido que este filtro ideal tiene dos problemones:

  • No es causal. La función conv sólo funciona para sistemas causales. Asume siempre que las señales a convolucionar tienen su soporte en los tiempos positivos.

  • No es de soporte acotado. No es FIR. La función conv recibe un vector finito, que representa a una señal de soporte acotado.

Se le ocurre una idea: "¿Y si desplazo la respuesta al impulso a la derecha, tratando de que sólo queden ceros en los tiempos negativos?

  • El módulo del espectro no me cambia al desplazar así que sigue siendo un pasabajos ideal

  • Ok, cierto, se defasa el espectro, pero de forma lineal así que no es tan grave. La señal de salida va a estar demorada un poquito, pero más allá de eso va a estar filtrada con el filtro ideal. Puedo vivir con eso."

Pero -- siguió pensando -- esto no me arregla todos los problemas. El soporte de la sinc es infinito: no importa cuánto desplace, nunca voy a conseguir que quede un sistema causal. Además está el asunto del soporte no acotado.

Aún acostado en la bañera con agua calentita, seguía con las cejas fruncidas. Esa sinc. Cómo se achican las bochas a medida que te alejás de n=0. No puede ser que bochas tan chiquitas me compliquen tanto la vida. ¡Si sólo fueran cero en vez de casi cero!

EUREKA!

"¿Y si las hago cero a la fuerza? Ok, no será el filtro ideal, pero tan diferente no puede ser. Dejemos las bochas del medio que son las que tienen casi toda la energía y hagamos cero las de un costado. Luego, a esa señal con soporte acotado la puedo desplazar para lograr que el sistema sea causal".

Y así nació el Método de Ventaneo.

Método de Ventaneo

La propuesta de Charlie es la siguiente:

  1. Armar el filtro pasabajos ideal, analíticamente.

  2. Multiplicarlo por una ventana centrada. Acá hay parámetros

    1. El tipo de ventana

    2. El ancho de la ventana. Si va a tener un soporte de N puntos, irá desde inline_formula not implemented a inline_formula not implemented. Note que N debe ser impar pues debe ser una ventana centrada). A esta N se la conoce como el orden del filtro.

  3. Desplazar a derecha para convertirlo en causal; i.e., desplazar inline_formula not implemented puntos.

Análisis

Con los conocimientos de la materia, podés analizar fácilmente cual es la respuesta en frecuencia real de este filtro no ideal.

El paso 2 modifica el espectro del filtro ideal, convolucionándolo con la transformada de la ventana. Esto introduce los artefactos que ya nos son familiares:

  • Los flancos ascendentes y descendentes instantáneos del filtro ideal dejan de serlo. El tiempo que le lleva pasar de un valor a otro está relacionado con el ancho del lóbulo principal de la transformada de la ventana, y esto depende fuertemente del ancho N de la ventana, y en menor medida del tipo de ventana.

  • Los lóbulos secundarios de la transformada de la ventana se traducen en sobrepicos y bajopicos en la respuesta del filtro ideal (efecto Gibbs). Las amplitudes dependen exclusivamente del tipo de ventana; la velocidad a la que se atenúan estas amplitudes sucesivas, del tipo de ventana y del orden del filtro.

#= Gráfico horrible y mal etiquetado sólo para mostrar cualitativamente cómo difieren las ventanas en ancho de lóbulo principal, ritmo de decrecimiento, y 
amplitud de primer lóbulo secundario =#
let
  normalize(x) = x ./ sum(x)
   plot()
   plotwin(win) = plot!(log.(abs.(fft(padright(normalize(win(41)), 501)))); 
      xlims=(0, 250), ylims=(-15, 2), ylabel="dB")
  
  foreach(plotwin, [ones, triang, hamming, hanning, blackman])
  plot!()
end
1.0s
Julia
→

El paso 3, introduce un delay en el filtro, igual al desplazamiento de inline_formula not implemented. El resultado final es un filtro de fase lineal (la función de fase es una función lineal); estos filtros desplazan temporalmente a todas las componentes de frecuencia en el mismo valor y entonces no se dispersa temporalmente la energía de la señal de entrada.

Ejemplo

order = 231
ns    = range(-(order - 1) / 2; length=order)
winds = hamming(order) 
hs    = h_ideal.(ns) .* winds
stem(0:order-1, hs)
0.6s
Julia
→
let
  samples = 500
  hfs = fft(padright(hs, samples))
  fs = range(0; step=2pi/samples, length=samples)
  
  plot(fs, abs.(hfs); xlims=(0, π))
  plot!(H_ideal, 0, π)
end
1.1s
Julia
→

Pasaaltos

Armemos uno con frecuencia de corte pi/8 y orden 231.

Se puede seguir el mismo método que antes, o convertir un pasabajos en pasaaltos.

Mismo método que antes

h_hp_ideal(n) = 7/8 * sinc(n * 7/8) * (-1)^n
order = 231
ns    = range(-(order - 1) / 2; length=order)
winds = hamming(order) 
hs_hp    = h_hp_ideal.(ns) .* winds
stem(0:order-1, hs_hp)
0.8s
Julia
→
plot(range(0; stop=2pi, length=500), 
  abs.(fft(padright(hs_hp, 500))); xlims=(0, pi))
0.5s
Julia
→

Convirtiendo un pasabajos en pasaaltos

Si le hago un desplazamiento en frecuencia al pasabajos anterior, que dejaba pasar entre -pi/4 y pi/4, y lo desplazo en pi, va a dejar pasar entre 3/4 pi y 5/4 pi

hs_hp2 = hs .* (-1) .^ range(0; length=length(hs))
plot(range(0; stop=2pi, length=500), 
  abs.(fft(padright(hs_hp2, 500))); xlims=(0, pi))
0.5s
Julia
→

Pasabandas

Nuevamente, se puede hacer directo, o a partir de otros filtros.

Combinando pasabajos y pasaaltos

Los pongo en cascada y sus transferencias se multiplican. La cascada va a eliminar lo que cualquiera de los dos elimina, y sólo va a dejar pasar lo que ambos dejan pasar.

let fs = range(0; stop=2pi, length=500)
  plot(fs, abs.(fft(padright(
          hs_hp, 
          500))))
  
  plot!(fs, abs.(fft(padright(
          hs, 
          500))))
  
  plot!(fs, abs.(fft(padright(
          conv(hs, hs_hp), 
          500))); xlim=(0, pi))
end
0.6s
Julia
→

Mismo método que antes

Armo la ideal, ventaneo y muestreo, entre pi*2/5 y pi*4/5

h_bp_ideal(n) = 2/5 * sinc(n * 1/5) * cos(pi*3/5*n)
order = 231
ns    = range(-(order - 1) / 2; length=order)
winds = hamming(order) 
hs_bp    = h_bp_ideal.(ns) .* winds
stem(0:order-1, hs_bp)
0.7s
Julia
→
plot(range(0; stop=2pi, length=500), 
  abs.(fft(padright(hs_bp, 500))); xlims=(0, pi))
1.5s
Julia
→

¿Restando pasabajos?

Ojo con la fase

h_lp_ideal(n; Ωc) = Ωc/π * sinc(n * Ωc/π)
order = 231
ns    = range(-(order - 1) / 2; length=order)
winds = hamming(order) 
# los hago del mismo orden
h_lp1 = h_lp_ideal.(ns; Ωc=2/5*pi) .* winds
h_lp2 = h_lp_ideal.(ns; Ωc=4/5*pi) .* winds
h_bp = h_lp2 .- h_lp1
plot(
  plot(range(0; stop=2pi, length=500),
     abs.(fft(padright(h_bp, 500))); xlims=(0, pi), title="|H|"),
  plot(range(0; stop=2pi, length=500),
     unwrap(angle.(fft(padright(h_bp, 500)))); xlims=(0, pi), title="fase");
  layout=(2, 1)
  )
0.7s
Julia
→

Btw, acá la pendiente (negada) de la fase te da el delay del filtro. Debería ser (231-1)/2=115 en este caso. A ojo tiene pinta de que da bien.

# Pero si uso órdenes diferentes para cada uno...
order = 181
ns    = range(-(order - 1) / 2; length=order)
winds = hamming(order) 
h_lp2 = h_lp_ideal.(ns; Ωc=4/5*pi) .* winds
h_bp = padright(h_lp2, length(h_lp1)) .- h_lp1
plot(range(0; stop=2pi, length=500),
     abs.(fft(padright(h_bp, 500))); xlims=(0, pi), title="|H|")
0.5s
Julia
→

Runtimes (1)