PRÁCTICA 1.
INTRODUCCIÓN AL TRATAMIENTO DE
SEÑALES ECG
El objetivo de esta práctica es trabajar el tratamiento de señales
electrocardiográficas (ECG) con el objetivo de familiarizarnos con el
procesamiento digital de señales mediante el uso de MATLAB. En esta práctica,
se emplean registros ECG de 12 segundos obtenidos de un paciente con
síndrome de Brugada, a una frecuencia de muestreo de 1000 Hz, utilizando un
sistema de grabación Holter.
Debemos cargar, filtrar y analizar estas señales, incluyendo operaciones como la
obtención de la transformada discreta de Fourier y la densidad espectral de
potencia, entre otros análisis.
Para ello, trabajaremos los siguientes apartados;
a) Cargar las señales. Puede usar el comando load para cargar las
variables ecg y fs que están el fichero ecg_12leads_12sec.mat. → >>
load ecg_12leads_12sec
Cargamos las variables usando el siguiente comando en el script y ejecutamos.
En primer lugar, cargamos el fichero con los datos, luego la frecuencia de
muestreo. Creamos un vector de tiempo que vaya desde la muestra 1, pero
antes usamos la función size para las 12 derivaciones.
% apartado a)
% load("ecg_12leads_12sec.mat")
load ecg_12leads_12sec.mat
fs = 1000;
[N, leads] = size(ecg);
t = (1:N)'/fs;
figure(1);
plot(t,ecg)
xlabel('tiempo(s),'), ylabel('mV')
Ejecutamos y nos dará la siguiente grafica;
Figura 1. El tiempo en segundos vs amplitud en mv de todo el
electrocardiograma (los 12 canales)
A continuación, creamos otra figura usando el detrend;
ecg = ecg - mean (ecg);
ecg = detrend(ecg);
figure(2)
plot(t,ecg)
xlabel('tiempo(s),'), ylabel('mV')
title('ECG Ideal (derivación) 2')
grid minor
Con la función detrend quitamos la tendencia del polinomio número 1 en el caso
de que vaya subiendo. Aparecen las 12 derivaciones solapadas.
Le ponemos una rejilla usando el grid minor para que se vea mejor la gráfica
(con la cuadricula detrás).
Figura 2. Muestra las 12 derivaciones en función del tiempo
% extraemos lead II, v4
sII = ecg(:,2);
sV4 = ecg(:,10);
Se extraen las señales de todas las filas de la columna 2 y se extraen de la
columna 10, ya que son las posiciones en las cuales se encuentran la señal 2 y la
v4 respectivamente.
b) Obtenga y grafique la señal ECG filtrada paso alto a 0.6 Hz y paso
bajo de 40 Hz de la derivación II. Puede usar un filtro Butterworth de
orden 3 (ayuda: ver Nota 3).
Para implementar un filtro en Matlab usamos el filtro butterworth que tiene 2
parámetros de entrada. Usamos los parámetros Ba, Bb, Ab y Aa para ejecutar la
señal paso bajo (40Hz) y luego el paso alto (0,6Hz).
La función filfilt hace un filtrado hacia delante y hacia detrás.
%% apartado b filtrado
fN = fs/2; % fN es la frecuencia de Nyquist la máxima que podemos recuperar
[Bb,Ab]=butter(3,40/fN); % paso bajo
[Ba,Aa]=butter(3,0.6/fN,'high'); %paso alto
sIIfa = filtfilt(Ba,Aa,sII);
sIIf = filtfilt(Bb,Ab,sIIfa);
figure(3)
plot(t,sII)
xlabel('tiempo(s)')
ylabel('Amplitud(mV)')
legend('ECG lead derivación II')
grid, grid minor
Figura 3. Señal filtrada de la derivación II.
figure(4)
plot(t,sII,t,sIIf)
xlabel('tiempo(s),'), ylabel('Amplitud(mV)');
legend('ECG lead II original', 'ECG lead IIfiltrado')
title('ECG lead II original vs lead filtrada');
grid, grid minor
Figura 4. Derivación II y V4 original y filtrada del 1 al 12
c) Obtenga el valor máximo y mínimo de amplitud del ECG de las
derivaciones II y V4 del primer latido, calculando también en qué
muestras de tiempo ocurren.
%% apartado c max y min de leads II y V4
sIIbeat = sIIf(430:1320);
sV4 = ecg(:,10);
sV4beat = sV4(430:1320);
figure(5)
plot(sIIf)
[x, ~] = ginput(2); %marcar de cruz hecha por mouse 2 puntos inicio y final
x = round(x);
sIIbeat = sIIf(430:1320);
tb = (1:length(sIIbeat))'/fs;
plot(tb, sIIbeat)
[MII,posMII] = max(sIIbeat);
posMIIseg = posMII/fs;
Cogemos el latido de la sII filtrada desde el valor máximo de 1320 y mínimo a
430 (a ojo).
Pero con la función getinput podremos coger los dos puntos con el ratón, y los
cogemos del eje x ya quw es el que nos interesa.
Figura 5. Marcamos con el ratón aproximadamente el punto de inicio y final.
A continuación, representamos el valor máximo y mínimo del latido en la zona
marcada anteriormente;
figure (6)
sV4beat = sV4(x(1):x(2));
sIIbeat = sIIf(x(1):x(2));
tbeat = (1:length(sIIbeat))'/fs;
plot(tbeat, sIIbeat, tbeat, sV4beat)
xlabel('Tiempo(s)')
ylabel('Amplitud(mV)')
grid
grid minor
Figura 6. Máximo y mínimo del latido del síndrome de Brugada que se ve en
color rojo, da una forma de “aleta de tiburón”
d) Obtenga y presente gráficamente la magnitud de la transformada
discreta de Fourier de la derivación II (ayuda: ver Nota 1).
%%d) Transformada de fourier en el lead II
XII=fft(sII);
M = abs(XII);
figure(7)
f=(0:N-1*fs/N);
plot(f,M)
xlabel('frecuencia (Hz)')
ylabel('(mV)')
En la señal filtrada, le vamos a aplicar la transformada de Fourier. EN primer
lugar calculamos la transformada en el XII y luego sacamos el valor absoluto M.
Figura 7. Frecuencia de muestreo de la transformada de Fourier aplicada en la
derivación II.
Graficamos 1000 muestras por segundo, el espectro de frecuencias es desde 0
hasta la mitad de frecuencia de muestreo (Nquist).
Es una manera de mostrarlo del espectro desde 0 hasta 500, y la segunda parte
seria como la imagen especular.
e) Obtenga y presente gráficamente la densidad espectral de potencia
(PSD) en dB de toda la señal de derivación II y de un latido de la
derivación II (ayuda: ver Nota 2).
%%e) densidad espectral de potencia o pediodograma
% nota 2
x = fft(sII);
Px = (abs(x).^2) / (N*fs);
f=(0:N-1)'*fs/N;
%para graficar la PSD:
figure(8)
plot(f,Px)
xlabel('frecuencia(Hz)')
ylabel('PSD(mV^2/Hz)')
if mod (N,2) == 0
select = 1:N/2+1;
else
select = 1:(N+1)/2;
end
f_unlado=f(select);
Px_unlado=2*Px(select);
Figura 8. Densidad espectral de potencia.
Representamos la densidad espectral de potencia vs a la frecuencia.
figure(9) % un latido de sIIf: sIIbeat
plot(f_unlado, Px_unlado)
xlabel('frecuencia (Hz)')
ylabel('PSD(mV^2/Hz')
Nb = length(sIIbeat);
tb = (1:Nb)'/fs;
plot(tb,sIIbeat), grid
xlabel('tiempo(s)')
ylabel('amplitud(mV)')
grid
grid minor
Figura 9. Representación del latido frente al tiempo.
figure(10)
plot(f_unlado,10*log10(Px_unlado))
xlabel('frecuencia(Hz)')
ylabel('PSD(dB/Hz)')
ylim([-120 0])
grid
[Px,f1] = periodogram(sIIbeat,[],Nb,fs);
[Px2,f2] = periodogram(sIIbeat,hamming(Nb),Nb,fs);
Aplicamos la ventana de Hamming a toda la señal.
Figura 10. Representación de un solo lado del PSD hasta 500.
figure(11)
plot(f1,10*log10(Px),f2,10*log10(Px2))
grid
xlabel('frecuencia(Hz)')
ylabel('PSD (dB/Hz)')
legend('Pediodograma estándar',' Periodograma Hamming')
title('Latido de derivación II')
Para representar el periodograma tenemos que aplicar el logaritmo en base 10
de Px1 y de Px2.
Figura 11. Representación del periodograma estándar frente al de Hamming.
En conclusión, esta práctica permitió aplicar conceptos fundamentales de
procesamiento digital de señales al analizar y procesar señales ECG reales en
MATLAB. A través de procesos como el filtrado, el análisis en frecuencia y la
obtención de parámetros característicos del ECG, se logró comprender cómo las
herramientas matemáticas y de programación permiten extraer información
relevante de señales biomédicas.
El uso de filtros pasabajos y pasoaltos permite mejorar la calidad de las señales
eliminando el ruido, mientras que la transformada de Fourier y el análisis de la
densidad espectral de potencia permitieron estudiar las frecuencias dominantes
en el ECG, fundamentales para el diagnóstico y monitoreo de condiciones de
latido cardíaco.
Finalmente hay que añadir que, esta práctica nos ha ofrecido una comprensión
práctica de cómo el tratamiento de señales ECG es esencial en el ámbito de la
bioingeniería.
