Código Matlab para acondicionar y filtrar señales de electrocardiogramas.


Sergey Chernenko nos muestra, explica y regala el código Matlab para detectar una característica importante de las señales provenientes de un electrocardiógrafo: el pico R de un electrocardiograma.

Para estudiar este excelente ejemplo de procesamiento de señales en ciencias medicas, lo mejor es bajar el archivo zip, con el código original.

Hay que estudiarlo para comentarlo en la clase.

Como nos comenta el mismo Sergey en su sitio Web, la tarea básica en el procesamiento de electro-cardiogramas (ECG) es la detección de los picos R. Sin embargo, debido a la respiración del paciente, la detección es complicada; pues se presenta en los ECGs artefactos que complican la detección de los picos R. Por ejemplo, las distancias irregulares entre los picos, la forma irregular en los picos, la presencia de componentes de baja frecuencia en los ECGs.

Es por ello que se deben contar con varias etapas de procesamiento que reduzcan la influencia de estos factores. Sergey nos muestra un caso donde una señal ECG original se le eliminan las frecuencias bajas por medio de FFT (Fast-Fourier-Transform). En la segunda etapa, se usa un filtro pasa banda, el cual localiza los máximos significativos en la señal. En la etapa final, se eliminan los datos más alejados del promedio. Este proceso permite la detección de los picos R.

Código para obtener la primera derivada por coeficientes de una función

Como se puede consultar en el libro: transformada de Fourier para peatones, podemos definir la derivada Yk como la “primera diferencia central” de la función F:


La cual tiene una función de transferencia, que no presenta problemas de corrimientos de fase:


Sin embargo, presenta un efecto de filtro muy pronunciado a frecuencias altas.

El siguiente código matlab aplicado a la función sin^2(x) con un paso de 0.1 funciona muy bien. En las figuras de salida del script se observan la función original, la derivada teórica y la derivada obtenida. Observe que presenta un comportamiento muy similar las dos derivadas. Efectivamente, cuando comparamos la derivada obtenida sobre la deriva teórica obtenemos en la gráfica una recta identidad. Por todo lo anterior, afirmamos que es un buen algoritmo para operar la derivada.

%%Script que calcúla numéricamente la primera derivada de la función:
%%sin(x)^2.
clear all
close all
clc
%Definiendo intervalos
x=0:0.01:1;
%función original:
funcion=(sin(x)).^2;
%preparando parámetros:
s=length(funcion); %longitud de la funcion
Y=zeros(1,s); % espacio de memoria para la derivada numerica
% Algoritmo de primera derivada
for k=2:s-1
Y(k)=(funcion(k+1)-funcion(k-1))/(2*0.01);
end
% Graficación de funciones teoricas y su compraración con la estimación
% numerica de la primer derivada
figure
subplot(2,1,1)
hold on
plot(x,funcion,'r') % función original
plot(x,2*sin(x).*cos(x), '-b') % derivada teorica de la funcion original
plot(x,Y,'ok') % derivada numerica obtenida
axis([0,1.04,0,1.1]);

xlabel('Vector <>')

ylabel('Funciones utilizadas')
text(0.02,0.8,'Funciones Sin^2(x) (rojo), derivada teorica (azul), derivada numerica (negro)')
box on
hold off 
subplot(2,1,2)
plot(2*sin(x).*cos(x),Y,'-b')% comparacion entre funciones;
% derivada teorica contra la derivada numerica
xlabel('Función 2*sin(x).*cos(x), derivada teorica')
ylabel('Derivada numerica obtenida')
axis([0,1.04,0,1.1]);
grid
%%% fin del script 


Ejercicios:

Matlab tiene su comando/función para hacer la derrivada: diff, su algoritmo de trabajo es:


1) ¿Cómo se compara este algoritmo con el propuesto en el post?

2) Muestre ejemplos del uso de este comando

3) ¿Cuales son las limitaciones de esta forma de definir la derivada?

Enlaces relacionados:

Capítulo 2, Derivación numérica del curso de métodos numéricos de Agustín A. Rosas.

Comparación del desempeño de un filtro convolución y un filtro pasa bajas de promedios ponderados en una misma función/señal

Ahora mostramos dos ejemplos de filtros. La funcion a trabajar es:




La cual tiene un pico extendido en t=5 y dos pequeños picos de estructura en t=2. Pues bien la entrada a un instrumento  cuenta con  filtro descrito por:




Efectivamente, la señal filtrada por  convolución será la transforma de Fourier inversa de la multiplicación de las transformadas de Fourier de la señal y del filtro; es decir:



A continuación aplicaremos el filtro en código Matlab.

Filtro 1: Convolución.

Para comparar la señal y la salida filtrada, a ambas las normalizamos. Observamos que hay un corrimiento en el tiempo. Este efecto en los indices y en la amplitud son muy comunes cuando se usa la transformada de Fourier.



Adicionalmente, una vista cercana a la estructura de los dos picos pequeños muestra que la salida del filtro cuenta con estos picos más pequeños y casi sin desernir. Con todo, esta forma de presentar los datos muestra a detalle que el filtro es bueno sin ser excelente para eliminar la estructura de la curva principal .


A continuación, mostramos el código Matlab utilizado para hacer estas curvas.


 %%% Inicio código 1: Filtro convolución
 clc
clear
close all

x = 0:0.0001:10;

y = 4*((sin(x-2)).^2).*exp(-50*((x-2)).^2) + 0.5*( ((x-5).^2 +0.25).^(-.25) );

r = ones(size(x));
for k= 1:length(x)
    if x(k) <= 0.15
        r(k)= .50;
    elseif x(k) <= 0.30
        r(k)= 0.25;
    elseif x(k) <= 0.60
        r(k)= 0.10;
    else
        r(k)= 0;
    end
end
plot(x,r,  'LineWidth',3)
axis([-.2,3,-0.2,0.6])
Y= ifft( fft(y).*fft(r));
figure
hold on
plot(x, y/max(y), '--b','LineWidth',3)
plot(x,Y/max(Y), 'k', 'LineWidth',3)
axis([0,10,0.28,1.05])
hold off
%%% Fin código 1: Filtro convolución

Filtro 2: Pasa bajas con promedio ponderado


Ahora bien, utilizaremos el filtro pasa bajas (PB), que ya mostramos en otro post, en este caso usamos 6,0000 iteraciones. Dado que la estructura en la función puede representar una oscilación rápida aplicaremos el filtro PB.   Observamos en la salida que tanto el cambio en la amplitud como el corrimiento en los indices de los vectores de las curvas es casi minimo. Además de que se ve más disminuida la estructura en la función.



Este filtro PB presenta mejores resultados que el de convolución; esto es por el tipo de filtro r(t) utilizado. Sin embargo, el filtro PB  tarda mucho más en desarrollarse. Puedes  usar los comandos de Matlab TIC yTOC antes y despues de que inicien las rutinas para observar la velocidad de ejecución de cada función o escript.


Podemos asegurar que dependiendo que se busquemos filtrar, hay un filtro que nos funcionara mejor que otros, pero hay que saberlo encontrar :)

A continuacion el segundo código Matlab:

%%% Inicio código 2: filtro pasa bajas con promedio ponderado
clc
clear
close all

tic % inicia el contador de tiempo
t = 0:0.001:5*pi;
f = 4*((sin(t-2)).^2).*exp((-50)*(t-2).^2) +0.5*((t-5).^2+0.25).^(-0.25); %señal de trabajo

hold on
plot(t, f, 'r') % para visualizar la señal con ruido

for corrida = 1:60000 %iteracion del filtro
for k = 2:(length(t)) -1 ;
f(k) = 0.25*(f(k-1) + 2*f(k) + f(k+1)); %filtro
end
end

plot(t, f, '--b') %para visualizar la señal filtrada
hold off

toc %finaliza el contador del tiempo

%%% Fin código 2: filtro pasa bajas con promedio ponderado

Ejercicios: 

1) Muestre el grado de atenuación de las protuberancias después de pasar por el filtro.

2) Diseñe otro filtro de convolución que atenue mejor a los dos picos que forman protuberancias en la curva original

3) Calcule la salida del filtro de convolución con r(w)*I(w) en lugar de  I(w)*r(w). Describa sus observaciones.

Ejemplos de código para formar señales discretas en Matlab

El objeto más básico en Matlab es una matriz numérica con la posibilidad de almacenar números complejos. Por supuesto, los datos obtenidos en el estudio de señales y sistemas son muy bien representados en forma de matrices. En este post usaremos Matlab para la generación de señales elementales: cuadrada, triangular, entre otras.

El ToolBox de procesamiento de señales de Matlab posee una gran variedad de funciones para la generación de señales, estas señales requieren de una representación vectorial de la variable tiempo, de manera continua o discreta. Para realizar una simulación de un intervalo continuo, se usa un vector de valores discretos con un intervalo de muestreo muy pequeño.

Como vimos en post anteriores, el siguiente comando genera un vector llamado t de valores que representan la variable tiempo, con un intervalo de muestreo de 1 ms entre 0 y 1 segundo.

t = 0:0.001:1;

Después de creado el vector que representa la variable tiempo, es posible iniciar el desarrollo de alguna señal de interés.

En Matlab una señal discreta en el tiempo se representa exactamente, porque los valores de la señal son representados como los elementos de un vector. Sin embargo las señales de tiempo continuo en Matlab son tan solo aproximaciones. La aproximación consiste de un vector cuyos elementos son muestras de la verdadera señal de tiempo continuo. Cuando se usa esta técnica para la representación de señales continuas es importante escoger el intervalo de muestreo lo suficientemente pequeño para asegurar que las muestras capturan todos los detalles de la señal.

EJEMPLOS DE SEÑALES EN MATLAB 

La generación de señales periódicas tales como ondas cuadradas y triangulares es una actividad muy fácil de realizar en MATLAB.

1) SEÑAL CUADRADA
Consideremos primero la generación de una onda cuadrada de amplitud A, frecuencia fundamental w (medida en radianes por segundo) y ciclo útil rho. Recordemos que el ciclo útil es la fracción de cada periodo en donde la señal es positiva.


Para generar dicha señal se pueden escribir lo siguiente en la linea de comandos:


>> A = 1;
>> w = 10 * pi;
>> rho = 0.5;
>> t = 0:0.001:1;
>> sq = A*square(w*t+rho);
>> plot(t,sq);

En la segunda línea de este ejemplo, pi es una función interna de Matlab que calcula el número más cercano a la constante PI en formato de coma flotante. El último comando es usado para vizualizar la señal generada. El comando plot dibuja líneas conectando los valores sucesivos de la señal y así da la apariencia de una señal en tiempo continuo.

2) SEÑAL TRINGULAR 

Consideremos ahora la generación de una onda triangular de amplitud A, frecuencia fundamental w y ancho Wdt . El periodo de la onda triangular será T con el máximo valor de la señal ocurriendo en t = WT . El comando básico para generar esta señal es:

A * sawtooth(w * t + Wdt)

El resultado se puede observar en la gráfica a la izquierda


3) SEÑAL ESCALÓN

En Matlab, el comando ones(M, N) genera una matriz de unos de tamaño MxN, y el comando zeros(M, N) es una matriz de ceros del mismo tamaño. Se puede hacer uso de estas dos matrices para generar dos señales comúnmente usadas: la señal escalón y la señal impulso.

Una señal paso de amplitud uno, puede ser generada con el siguiente comando.


U = [zeros(1, 10), ones(1, 11)];


Para la versión continua creamos un vector que represente el tiempo el cual tenga muestras de un intervalo separados por valores muy pequeños. Los comandos y los resultados se muestran a continuación:

>> u=[zeros(1,10),ones(1,11);
>> t=-1:0.1:1;
>> plot(t,u)


Como se menciono anteriormente, una señal generada en Matlab es inherentemente de naturaleza discreta. Para visualizar una señal en tiempo discreto se puede hacer uso del comando stem. Específicamente stem(n, x), bosqueja los datos contenidos en el vector x como una señal de tiempo discreto con los valores de tiempo definidos por el vector n. Los vectores n y x deben tener dimensiones compatibles, es decir deben tener el mismo número de elementos. Así, para este caso para obtener la representación de esta señal en tiempo discreto creamos un vector-tiempo el cual debe tener valores separados por una unidad.

>> u=[zeros(1,10), ones(1,11)];
>> n=-10:10;
>> stem(n,u)

Recuerde que para poder usar las funciones plot y stem, es requisito que los vectores (t y u) ó (n y u) tengan iguales dimensiones. Por esta razón el vector u se forma como una composición de diez ceros y 11 unos, debido a que los arreglos t y n, tienen dimensión 21 dado que incluyen un elemento central el cual es el número cero. Para probar este hecho, se puede hacer uso de la función Matlab llamada size que devuelve como resultado un vector con las dimensiones de la matriz que se le pasa como parámetro.

4) SEÑAL IMPULSO:


La versión discreta de la señal impulso se puede también generar con ayuda de las funciones zeros y ones, realizando una composición como sigue:


>> delta = [ zeros( 1 ,10 ), 1 , zeros( 1 ,10 ) ];
>> n = -10:10;
>> stem(n,delta);



5) SEÑAL RAMPA 

Para generar la señal rampa, tan solo es necesario recordar que esta función puede ser creada, como la composición de una recta Y(x) = x a partir de cero y de la recta Y(x) = 0 para valores de x menores de cero, así la versión discreta se muestra a continuación:


>> t1=0:0.1:10;
>> rampa1=t1;
>> rampa=[zeros(1,101),rampa1];
>> t2=-10:0.1:0;
>> t=[t2,t1];
>> plot(t,rampa)

Ejercicios:

1) Desarrollar un conjunto de comandos Matlab para aproximar las siguientes señales periódicas en tiempo continuo, dibujando 5 ciclos de cada una:

a) Onda Cuadrada, de amplitud 5 Volts, frecuencia fundamental 20 Hz y ciclo útil del 60%.

b) Señal diente de sierra, amplitud 5 Volts y frecuencia fundamental 20Hz

Práctica 2: Señales harmónicas moduladas y su filtrado

Ya pueden consultar su práctica 2, para lo cual les recomiendo que leen sus libros sobre transformada de Fourier.

Por cierto, esta presentación sobre función de transferencia les será útil en los siguientes semanas.









Por cierto, también deben consultar el artículo: How can Trigonometry help to understand Lock-in Amplifier operation? De Ernesto Marín, R. Ivanov. Pues la técnica lock-in también es un método de filtrado de señales, muy aplicado cuando la razón señal/ruido es muy baja.

Código y ejemplos con diferentes iteraciones de un filtro pasa bajas de promedio ponderado


Estas dos figuras se obtuvieron con pequeñas modificaciones del código que aquí se presenta. En el primer caso se realizaron 10,000 iteraciones del mismo filtro. El tiempo que tomo fue considerable y su filtradrado fue bueno. Sin embargo con 1,000 iteracciones el filtrado es malo y el tiempo no se reduce considerablemente, lo cual es lógico pues disminuyo en un factor de 10 las iteraciones, aunque siguen siendo muchas iteraciones 1,000. Estos son sólo ejemplos de cómo aplicar filtros de promedios ponderados hay filtros alternativos que son más rápidos y eficaces en la tarea asignada. El ejercicio sirve para ver las entrañas de estos filtros.


% Inicio de código
clc
clear
close all

x = 0.0:0.001:3*pi;

yy = sin(0.5*x)+sin(40*x); % la frecuencia a limpiar es la de 40, la cual es relativamente grande

subplot(3,1,1)
plot(x,yy) % primero vemos como es la señal original

% limpiar señal

% y(k) = yy(k-1) + 2*yy(k) + yy(k+1)

final = length(yy);

y = yy;
for m =1:10000 % hacemos muchas iteraciónes pues el filtro es poco eficiente

for k =2:final-1
y(k) = (0.25)*(y(k-1) + 2*y(k) + y(k+1)); % filtro pasa bajas
end
end

subplot(3,1,2)

plot(x,y) % el resultado lo vemos aqui.

subplot(3,1,3)
hold on
plot(sin(0.5*x),y) % el resultado lo vemos aqui.
plot(y,y, '--k') % el resultado lo vemos aqui.

hold off

Código para encontrar la señal principal mediante el filtro de la función max

Les dejo uno de los ejemplos que me enviaron para hacer la evaluación.


% Inicio  del código
%% ejemplo del uso de las funciones FFT y MAX para limpiar una señal

%% con ruido
%% 12 sept 2011

% limpieza de la memoria
clc
clear
close all

% Señales con ruido
x=0:1/255:1;
y= 10*sin(2*pi*29*x)+  ...
   3*sin(2*5*x)+ ...
   2.5*sin(rand*25*2*pi*x)+...
   3.5*cos(rand*155*2*pi*x)+...
   2.5*sin(300*pi*x);

% Transformada de Fourier simple
Y=fft(y);
[a,b]=max(Y);
%indiceAmplitudmax=b;
% haciendo el vector con ceros
Y1=Y*0;
% sustituyendo en el espacio del indice el numero complejo que representa
% el maximo
Y1(b)=(a);
%tomando unicamente la parte real de la parte inversa de la transformada de
%Fourier. La parte imaginaria no nos da más información en esta parte
inversa= real(ifft(Y1));

% visualizando en la misma grafica la señal con ruido y limpia
hold on
plot(x,y,'r');
plot(x,inversa, 'k')
hold off

% Fin del código

Pregunta para pensar:

¿qué implica en la señal de principal (color negro) que sea un número complejo el máximo del espectro en la Transformada de Fourier?,  ¿acaso son son dos señales?  o ¿cambios de fase?.

Códigos adicionales para gráficas y código de filtro pasabajas para mejorarlo

Les pongo el Código para la tarea, deben hacer que funcione el filtro.



%%%%%%%%%%%%%%%%%%
clc
clear
close all

%señal a trabajar

x = 0:0.0001:5*pi;
 y = sin(x) + sin(100*x);

 subplot(2,1,1)
 plot(x, y, 'c')


 for corrida = 1:6
     for k = 2:(length(x)) -1 ;
   
 y(k)  = 0.25*(y(k-1) + 2*y(k) + y(k+1));
     end
 end

 subplot(2,1,2)
  plot(x, y, 'y')
%%%%%%%%%%%%%%%%%%

Por otro lado, estos son los comandos extra para hacer otro tipo de gráficas:

X = [linspace(0,pi*3/8,20); linspace(pi*5/8, pi, 20)];
Y = sin(X); % Tanto la X como la Y son matrices
plot(X,Y) % Esto no es lo que se esperaba!!! Hay que teneren cuenta que se imprime por columnas


pie(X)
pie(X, Explode)

x = [1 1 2 3 5 8];
ex = [0 1 0 0 1 0];
pie(x, ex);

hist(X)
hist(X, nbins)
x = randn(1,10000);
hist(x);

% Generamos números aleatorios con distribución normal y
uniforme
x = [randn(10000, 1) rand(10000, 1)*6-3]
hist(x, 50);

x = 0.5:0.5:4;
y = 1./x;
barh(x, y)

Barras apiladas

x = 0.5:0.5:4;
y = 1./x;
Y = [y' fliplr(y)'];
bar(x, Y)

x = 0.5:0.5:4;
y = 1./x;
Y = [y' fliplr(y)'];
bar(x, Y, 'stacked')
x = randn(1,100);
y = randn(1,100);
scatter(x, y)

scatter(x, y, rand(1, 100)*1000, [1 0.5 0])
Related Posts Plugin for WordPress, Blogger...