Proyecto Especial

Detección de latidos cardíacos a partir de lagrabación de video de un smartphone

# tmp

2 + 10

Carga de datos

Cargar el archivo `intensidad_RGB.mat` y visualizar 1la señal RGB en superposición (3 canales). Elija alguno de los 3 canales que, según su criterio, posea la señal más útil a los efectos la dinámica de los latidos y utilícelo para realizar los puntos siguientes.

using DSP, SySTools

@load 
intensidad_RGB.jld2 sr x

println("La frecuencia de muestreo es de $sr.")

println("Los datos son una matriz de $(size(x)) donde cada columna es un canal RGB")

let

  cols = [:red :green :blue]

  ts = range(0, length = size(x, 1), step = 1/sr)

	plot(ts, x, color=cols, label=cols, xaxis="t (s)")

end

Tiene pinta que el verde tiene la data, lalala. Tiro a la goma el primer cachín

x_gr = x[100:end, 2]

;

Estimación manual de latidos a partir de la señal temporal

Del gráfico del punto anterior, estimar aproximadamente los latidos por minuto (LPM). Identifique el momento a partir del cual la frecuencia cardiaca comienza a incrementarse, y en cuánto se incrementa.

Estimación manual de latidos a partir del espectro

Rehacer el punto anterior, pero utilizando DFT. ¿Puede en este caso identificar cuándo se da el cambio de LPM?

let

  len = div(length(x_gr), 2)

  xf = abs.(fft(x_gr)[1:len])

  fs = range(0, stop=sr/2, length=len)

  plot(fs, xf, xaxis="f (Hz)", leg=false)

end

No pude :(

Diseño de filtro

Diseñar un filtro pasa-banda tipo Butterworth con banda de paso entre 0.5 Hz y 10 Hz.

bp = digitalfilter(Bandpass(0.5, 10, fs=sr), Butterworth(20))

;

Graficar respuesta en frecuencia (módulo y fase), diagrama de polos y ceros y respuesta al impulso

let

  fs = range(0, stop=sr/2, length=200)

  ys = abs.(freqz(bp, fs, sr))

  plot(fs, ys, xaxis="f (Hz)", leg=false, title="Módulo del filtro")

end

zplane(bp)

¿Qué papel juega el orden del filtro seleccionado en su diseño?

Filtrado y análisis

Filtrar la señal FPG utilizando el filtro diseñado en el punto anterior.

x_f1 = filt(bp, x_gr)

;

Grafique en superposición la señal original con la filtrada y comente acerca de:

let 

  ts = range(0, length = size(x_f1, 1), step = 1/sr)

	plot(ts, [x_gr x_f1], xaxis="t (s)", label=["Original" "Filtrada"])

end

Remoción de derivas

Blablabla derivas se fueron blah.

Cambios en la forma de la señal

Blablabla forma cambia un toque blah

Retardo de la señal filtrada respecto de la original

Blah bloh retardo un toque orden blah. Mirá

plot([x_gr x_f1][300:360,:])

Retardo temporal

A partir de la respuesta en fase del filtro, calcule su retardo temporal y compare con lo observado en el punto 5c.

let

  fs = range(0, stop=sr/2, length=200)

  ys = phase(freqz(bp, fs, sr))

  plot(fs, ys, xaxis="f (Hz)", yaxis="Rads",

      title="Fase del filtro", leg=false)

end

let

  fs = range(0, stop=sr/2, length=200)

  ys = phase(freqz(bp, fs, sr))

  plot(fs[50:90], ys[50:90], xaxis="f (Hz)", yaxis="Rads", 

    title="Zoom a la fase")

end

Pendiente de aprox balbal

Filtrado sin defasaje

Justificar teóricamente el funcionamiento de este tipo de filtrado y cuál resulta su ventaja.

Va y viene, par, fase lineal, blabla.

Implementar un filtrado IIR ida y vuelta para anular la fase del filtro. Filtrar nuevamente la señal FPG y comparar el resultado con lo obtenido en el punto anterior, particularmente en la forma de la señal y su retardo.

x_f2 = filtfilt(bp, x_gr)

plot([x_gr x_f1 x_f2][300:360,:], 

  label = ["Original" "Filtrado anterior" "Filtrado ida y vuelta"])

Realización de un espectrograma

Realizar un espectrograma de la señal antes y después de filtrar.

let

  wlen  = round(Int, sr*2)

  wtype = hamming

	specplot(x_gr, wlen, fs=sr, window=wtype, nfft=2^9, title="Sin filtrar")

end

let

  wlen  = round(Int, sr*2)

  wtype = hamming

	specplot(x_f2, wlen, fs=sr, window=wtype, nfft=2^9, title="filtrado")

end

Justificar la longitud de ventana elegida y comente acerca del resultado obtenido, relacionándolo con los puntos 2 y 3.

Calcule la resolución en frecuencia de la ventana mediante DFT.

¿Cómo haría para obtener mejor resolución en frecuencia y qué se pierde con esto?

Análisis del espectrograma

Identificar en el espectrograma la zona donde el pulso se acelera.

Observar con detenimiento los componentes de frecuencia que posee la señal y justificar el origen de cada uno (para esto último, necesitará hacer uso de la señal audio_det.mat para explicar todos los componentes observados).

Detector automático de latidos

Realizar un detector automático de latidos. El mismo debe tomar como entrada la señal FPG y producir como salida un vector de tiempos, donde cada tiempo corresponde a la detección de un latido en la señal. Para esto, se sugiere implementar los siguientes pasos:

Filtrado pasa-banda

filtrado(x) = filtfilt(bp, x)

;

Filtro de derivada

deriv(x) = circshift(filt([-2, -1, 0, 1, 2], x), -2)

;

Normalización con energía instantánea

Primero calcular la energía instantánea de la señal mediante un filtro MA1 de la señal del punto 10a elevada al cuadrado.

Lluego dividir la señal del punto b por el vector obtenido. Esto tiene como objeto reducir el impacto de la presión sanguínea sobre el nivel de señal.

mavg(n, x)    = filtfilt(ones(n)./n, x)

energia(x)    = mavg(10, x.^2)

normalizado(x) = deriv(x)./energia(x)

;

Remuestreo

Sobre-muestreo en un factor 4 para obtener mayor resolución temporal: implemente el sobre-muestreo utilizando la función upsample y diseñe un filtro interpolador FIR.

function resample(n, x; order=100) 

  h    = lowpass(1/n, order=order)

  x_up = upsample(n, x)

  return filtfilt(h, x_up)

end

const factor = 4

remuestreo(x) = resample(factor, x)

;

Grafique: respuesta en frecuencia del filtro en módulo y fase, y señal original y sobre-muestreada en superposición.

let

  ord = 100

  len = div(ord, 2)

  fs = range(0, stop=sr/2, length=len)

  ys = fft(lowpass(1/factor, order=ord))[1:len]

  plot(fs, [abs.(ys) phase(ys)], xaxis="f (Hz)", layout=(2, 1), leg=false)

end

Detector de picos

Detector de picos mediante umbral (puede definir como umbral un valor arbitrario).

# a ver cómo viene la señal hasta ahora

x_so_far = (remuestreo ∘ normalizado ∘ filtrado)(x_gr)

plot(x_so_far[1000:2000])

const umbral = 3e-7

function peaks(x) 

  pos_up_down = findall(v -> abs(v) == 1, diff(x .> umbral))

  pos = map(mean, Iterators.partition(pos_up_down, 2))

  return pos

end

function detect_beats(x, sr)

  f = peaks ∘ remuestreo ∘ normalizado ∘ filtrado

  raw = f(x)

  return raw ./ (factor*sr)

end

;

Resultados

Gráfico en superposición de la señal con las marcas de los picos detectados.

let

  ts = range(0, step = 1/sr, length=length(x_gr))

  ts_beats = detect_beats(x_gr, sr)

  vline(ts_beats)

	plot!(ts, x_gr, xlims=(40, 70), leg=false)

end

Sugerencia para el punto e: en general, un pico en la señal producirá la detección de múltiples muestras por encima del umbral. Para reducirlas a sólo una, puede utilizar la función diff para evaluar la primera derivada sobre el vector de muestras detectadas (la cual contiene sólo unos y ceros), quedándose solo con las muestras positivas que resulten y buscando, para cada una de estas, la posición del máximo inmediato de la señal original dentro de un intervalo pre-establecido (este deberá contener mínimamente a la duración aproximada de los picos).

Tiempo entre latidos

En base a los resultados del punto anterior, calcule y grafique el intervalo temporal instantáneo entre latidos (IBI: inter-beat interval) y los LPM instantáneos.

function bpm(x, sr, db=detect_beats)

  ts_beats = db(x, sr)

  bpms = 60 ./ diff(ts_beats)

  ts_avg = (ts_beats[1:end-1] .+ ts_beats[2:end]) ./ 2

  return ts_avg, bpms

end

plot(bpm(x_gr, sr)..., ylims=(50, 120), 

  xaxis="Tiempo (s)", yaxis="Latidos (BPM)", leg=false)

Opcional: mejorar el detector de latidos

Opcional. Mejorar el detector de latidos aplicando las reglas de [5]:

a. Establecer como regla que, si dos latidos se detectaron con una separación temporal menor a 200ms, sobrevive sólo aquel que corresponda al pico de la señal mayor entre ambos.

const min_beat_ts = 200e-3

function detect_beats_v2(x, sr)

	ts_beats = detect_beats(x, sr)

  bad_pos = findall(diff(ts_beats) .< min_beat_ts)

  goods = trues(length(ts_beats))

  for bad_cand in bad_pos

    bad_cand_ns = round.(Int, ts_beats[bad_cand:bad_cand+1] .* sr)

    second_is_worse  = x[bad_cand_ns[1]] > x[bad_cand_ns[2]]

  	goods[bad_cand + second_is_worse ] = false

  end

  return ts_beats[goods]

end

plot(bpm(x_gr, sr, detect_beats_v2)..., ylims=(50, 120), 

  xaxis="Tiempo (s)", yaxis="Latidos (BPM)", leg=false)

b. Establecer como regla que si el IBI instantáneo aumenta repentinamente en al menos 1.5 veces entre muestras consecutivas, puede haberse perdido la detección de un latido. Cuando este sea el caso, realice una nueva detección de picos dentro del intervalo correspondiente, utilizando un umbral la mitad del nominal. Si de esta manera se halla un nuevo pico, distanciado al menos 360ms de la detección precedente, entonces clasificarlo como latido.

Opcional: datos propios

Obtener su propia señal FPG mediante la cámara de su celular y aplicar los análisis y algoritmos desarrollados a esta señal. Para esto, realice una grabación tal como se explicó en la sección 1 y utilice los scripts auxiliares de Matlab para extraer los archivos de intensidades y audio del video. Nota: para un mejor resultado en la señal registrada, utilizar el dedo índice aplicando una leve presión sobre la lente, evitando resultar excesivamente suave o fuerte.