c
c ----------------------------------------------------------------------
C INT_QT.FOR BUSCA EL INICIO DE LA ONDA Q Y EL FIN DE LA T
C PARA MEDIR EL QT
C AUTOR: PABLO LAGUNA
C FECHA: 17-9-87
C ----------------------------------------------------------------------
C -----------------------------------------------------------------------
C Copyright (C) 2002 Pablo Laguna
C
C This program is free software; you can redistribute it and/or modify it
C under the terms of the GNU General Public License as published by the
C Free Software Foundation; either version 2 of the License, or (at your
C option) any later version.
C
C This program is distributed in the hope that it will be useful, but
C WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
C General Public License for more details.
C
C You should have received a copy of the GNU General Public License along
C with this program; if not, write to the Free Software Foundation, Inc.,
C 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
C
C You may contact the author by e-mail (laguna@posta.unizar.es) or postal
C mail (Dpt. Ingenieria Electrónica y Comunicaciones, Grupo de Tecnología
C de las Comunicaciones (GTC), Centro Politécnico Superior. Universidad
C de Zaragoza, María de Luna, 3 (Pol. Actur), Edificio A, Despacho 3.16,
C 50015 Zaragoza. Spain). For updates to this software, please visit
C PhysioNet (http://www.physionet.org/).
C _______________________________________________________________________
C Se llama desde ecgmain
C ----------------------------------------------------------------------
C BUSCA EL INTERVALO QT
C ----------------------------------------------------------------------
subroutine int_qt(sen,der1,iqrs,iondaq,iondat,ifm,nqrs,
& iqt,iqtc,iask)
dimension sen(100000),der1(100000),iqrs(8000),iondaq(8000)
dimension iqt(8000),iqtc(8000),iondat(8000)
if (iask.eq.2) then
write(6,*)
write(6,*) 'que fraccion quieres de la derivada?'
5 write(6,10)
10 format('$ der_limite = der_maxima / x x= ')
read(6,*,err=5) frac
end if
call frec_cardiaca(ifm,iqrs,nqrs,irrint)
J=0
do i=1,nqrs
call busca_q(sen,der1,iqrs(i),iondaq(i),ifm)
if (iask.eq.1) then
call busca_t1(sen,der1,iqrs(i),iondat(i),ifm,irrint)
else
call busca_t2(sen,der1,iqrs(i),iondat(i),ifm,irrint,frac)
end if
iqt(i)=(iondat(i)-iondaq(i))*1000/ifm
rr=(iqrs(i)-iqrs(i-1))/float(ifm)
if (i.eq.1) rr=(iqrs(2)-iqrs(1))/float(ifm)
call qt_corregido(rr,irrint,ifm,iqt(i),iqtc(i))
if (i.gt.2) call ajustrr(i,iqrs,irrint,no_uso)
end do
do i=nqrs+1,8000
iqt(i)=0
iondaq(i)=0
iondat(i)=0
iqtc(i)=0
end do
return
end
C ----------------------------------------------------------------------
C BUSCA EL INCIO DE LA ONDA Q
C ----------------------------------------------------------------------
subroutine busca_q(sen,der1,iqrs,iqb,ifm)
dimension sen(100000),der1(100000)
dimension der2(8000),der1f(8000)
c sen ES LA SEQAL EN LA QUE MIRO LA ONA Q
c der1 ES SU DERIVADA PRIMERA
c iqrs SON LAS POSICIONES LOS QRS DETECTADOS EN LA SEQAL sen
c iondaq SON LAS POSICIONES DEL INICIO DEL QRS
c der2 ES LA DERIVADA SEGUNDA SE sen EN ZONAS QUE INTERESA
c ifm ES LA FRECUENCIA A LA QUE ESTA MUESTREADA LA SEQAL
C nqrs ES EL NUMERO DE QRS QUE HAY SELECCIONADOS EN LA SEQAL
in=nint(ifm*0.04)
n=iqrs
if (sen(n).gt.0) then
call r_q(ifm,der1,iqrs,iq)
ir=iqrs
else
call posimaxmin(der1,imaxd,imind,iqrs,iqrs+in)
call posimaxmin(der1,imaxi,imini,iqrs-in,iqrs)
if (der1(imaxd).gt.(-der1(imini))) then
call detectar_cero(der1,n,'d')
ir=n
iq=iqrs
else
call detectar_cero(der1,n,'i')
ir=n
call r_q(ifm,der1,n,iq)
end if
end if
c SI HAY Q BUSCO DONDE EMPIEZA: PRIMERO VEO DONDE ES MAS PENDIENTE
C EL FLANCO DE BAJADA DE LA ONDA Q (j) Y LUEGO MIRO PARA ATRAS DONDE
C CAMBIA MAS RAPIDAMENTE LA DERIVADA-> MAXIMO DE LA DERIVADA 2
in=nint(ifm*0.04)
if (iq.gt.0) then
j=iq
do while (der1(j-1).le.der1(j))
j=j-1
end do
call derivada2(j-in,j,der1,der2,ifm)
l=j
call detectar_cero(der1,l,'i')
k=in
call detectar_cero(der2,k,'i')
if (l.gt.j-(in-k)) then
iqb=l
else
iqb=j-(in-k)
end if
if (iqb.lt.iq-ifm*0.06) iqb=iq
else
c PARTO HACIA ATRAS DE LA ONDA R POSITIVA
call fpbc(ir-ifm/10,ifm/10,ifm,30,der1,der1f,ret_pb)
j=ifm/10
do while (der1f(j-1).ge.der1f(j))
j=j-1
end do
idermax=ir-ifm/10+j
call derivada2(j-in,j,der1f,der2,ifm)
l=idermax
call detectar_cero(der1,l,'i')
k=in-1
call detectar_cero(der2,k,'i')
if (l.gt.idermax-(in-k)) then
iqb=l
else
iqb=idermax-(in-k)
end if
end if
return
end
c ------------------------------------------------------------------
C VA DE LA POSICON DE LA ONDA R A LA DE LA ONDA Q SI LA HAY
C ------------------------------------------------------------------
subroutine r_q(ifm,der1,ir,iq)
dimension der1(100000)
c BUSCO LA POSICION DELA ONDA Q 'IQ' 40 MS.
C ANTES DE LA ONDA R SI NO ESTA AQUI IQ=0
n=ir
call detectar_cero(der1,n,'i')
in=nint(ifm*0.035)
if (n.le.ir-in) then
iq=0
else
iq=n
end if
return
end
C -------------------------------------------------------------------
C CALCULA LA DERIVADA DEL TROZO DE SEQAL QUE MANDO
C -------------------------------------------------------------------
subroutine derivada2(ni,nf,x,y,ifm)
dimension x(100000),y(8000)
c ni,nf SON LOS PTOS. INICIAL Y FINAL DE LOS QUE CALCULO LA DERIVADA
C EN X PONGO LA SEQAL DE ENTRADA Y EN Y LA DE SALIDA
k=ifm/8
y(1)=(ifm/4)*(x(ni+2)+2*x(ni+1)-3*x(ni))
y(2)=(ifm/7)*(x(ni+3)+2*x(ni+2)-3*x(ni))
do i=3,nf-ni
y(i)=k*(x(ni+i+1)+2*x(ni+i)-2*x(ni+i-2)-x(ni+i-3))
end do
do i=nf-ni+1,8000
y(i)=0
end do
return
end
C --------------------------------------------------------------------
C DETECTA EL CERO SIGUIENTE DE LA SEQAL DE ENTRADA SEN A LA DERECHA
C O LA IZQUIERDA 'I' 'D' Y SACA LAPOSICION EN N
C --------------------------------------------------------------------
subroutine detectar_cero(sen,n,sentido)
dimension sen(100000)
character*1 sentido
if (sentido.eq.'i') then
i=n-1
do while (sen(i)*sen(i-1).gt.0)
i=i-1
end do
else
i=n+1
do while (sen(i)*sen(i-1).gt.0)
i=i+1
end do
end if
n=i
return
end
C --------------------------------------------------------------------
C SACA LA POSICION DEL MAXIMO Y DEL MINIMO DE LA SEQAL
C --------------------------------------------------------------------
subroutine posimaxmin(sen,imax,imin,ni,nf)
dimension sen(100000)
imax=ni
imin=ni
do i=ni,nf
if (sen(i).gt.sen(imax)) imax=i
if (sen(i).lt.sen(imin)) imin=i
end do
return
end
C ----------------------------------------------------------------------
C BUSCA EL FIN DE LA ONDA T SEGUN EL METODO DE LA DERIVADA SEGUNDA
C ----------------------------------------------------------------------
subroutine busca_t1(sen,der1,iqrs,ite,ifm,irrint)
dimension sen(100000),der1(100000),der1f(8000)
dimension der2(300)
c sen eS LA SEQAL EN LA QUE MIRO LA ONA Q
c der1f ES LA DERIVADA PRIMERA FILTRADA
c der1 eS SU DERIVADA PRIMERA
c iqrs SON LAS POSICIONES LOS QRS DETECTADOS EN LA SEQAL sen
c iondat SON LAS POSICIONES DEL INICIO DEL QRS
c der2 ES LA DERIVADA SEGUNDA DE sen EN ZONAS QUE INTERESA
c ifm ES LA FRECUENCIA A LA QUE ESTA MUESTREADA LA SEQAL
C nqrs ES EL NUMERO DE QRS QUE HAY SELECCIONADOS EN LA SEQAL
c IT ES LA POSICION DE LA ONDA t
it=0
call busca_it(ifm,irrint,sen,iqrs,der1,it)
c BUSCO EL PUNTO SE MAYOR DERIVADA NEGATIVA ENTRE IT Y IT+100MS.
call fpbc(it,ifm/3,ifm,30,der1,der1f,ret_pb)
call posimaxmin(der1,imax,j1,it,it+ifm/10)
call posimaxmin(der1f,imax,j2,1,ifm/10)
if (j1.gt.it+j2) then
m=j1
else
m=it+j2
end if
c BUSCO EL FINAL DE LA ONDA T
in=nint(ifm*0.2)
call derivada2(m-it,m-it+in,der1f,der2,ifm)
l=m
call detectar_cero(der1,l,'d')
k=2
call detectar_cero(der2,k,'d')
if (l.lt.m+k) then
ite=l
else
ite=m+k
end if
return
end
C ----------------------------------------------------------------------
C BUSCA EL FIN DE LA ONDA T SEGUN EL METODO DEL % DE DERIVADA
C ----------------------------------------------------------------------
subroutine busca_t2(sen,der1,iqrs,ite,ifm,irrint,frac)
dimension sen(100000),der1(100000)
dimension der1F(8000)
c sen eS LA SEQAL EN LA QUE MIRO LA ONA Q
c der1 eS SU DERIVADA PRIMERA
c iqrs SON LAS POSICIONES LOS QRS DETECTADOS EN LA SEQAL sen
c iondat SON LAS POSICIONES DEL INICIO DEL QRS
c der1f ES LA DERIVADA PRIMERA FILTRADA PASO BAJO
c ifm ES LA FRECUENCIA A LA QUE ESTA MUESTREADA LA SEQAL
C nqrs ES EL NUMERO DE QRS QUE HAY SELECCIONADOS EN LA SEQAL
c IT ES LA POSICION DE LA ONDA t
it=0
call busca_it(ifm,irrint,sen,iqrs,der1,it)
c BUSCO EL PUNTO SE MAYOR DERIVADA NEGATIVA ENTRE IT Y IT+100MS.
call fpbc(it,ifm/5,ifm,30,der1,der1f,ret_pb)
call posimaxmin(der1f,imax,j2,1,ifm/10)
derlim=der1f(j2)/frac
c BUSCO EL FINAL DE LA ONDA T
j=j2
do while (der1f(j).lt.derlim)
j=j+1
end do
ite=j+it
if (i.gt.2) call ajustrr(i,iqrs,irrint,no_uso)
return
end
C ------------------------------------------------------------
C BUSCA LA POSICION DE LA ONDA T
C ------------------------------------------------------------
subroutine busca_it(ifm,irrint,sen,iqr,der1,it)
dimension sen(100000),der1(100000)
c BUSCO LA ONDA t ENTRE 100 Y 400 MS. DESPUES DEL QRS
ilim1=nint(ifm*0.1)
if (irrint.lt.ifm*3/5) then
ilim2=nint(irrint/2.)
else
ilim2=nint(ifm*0.4)
end if
rmax=-abs(sen(iqr))
j=iqr+ilim1
do while (j.lt.iqr+ilim2)
call detectar_cero(der1,j,'d')
if (sen(j).gt.rmax) then
rmax=sen(j)
if (j.lt.iqr+ilim2) it=j
end if
end do
if (it.lt.iqr) it=j
return
end
C ------------------------------------------------------------
C CALCULA EL QT CORREGIDO
C ------------------------------------------------------------
subroutine qt_corregido(rr,irrint,ifm,iqt,iqtc)
c CALCULO EL QT CORREGIDO SI EL INTERVALO RR ESTA ACOTADO
if (rr.lt.irrint*1.1/ifm.and.rr.gt.irrint*0.9/ifm) then
iqtc=iqt/sqrt(rr)
end if
if (rr.lt.irrint*0.6/ifm) then
iqtc=0
endif
return
end
C -------------------------------------------------------------
C BUSCA LA FRECUENCIA CARDIACA MEDIA
C ------------------------------------------------------------
subroutine frec_cardiaca(ifm,iqrs,nqrs,irrint)
dimension iqrs(8000)
C BUSCO LA FRECUENCIA DEL CORAZON MAS COMUN
isum=0
max=0
min=iqrs(2)-iqrs(1)
nqrsini=10
if (nqrs.le.10) nqrsini=nqrs-1
icon=0
do i=1,nqrsini
irrint=iqrs(i+1)-iqrs(i)
c VEO SI LA FRECUENCIA CARDIACA ESTA EN (30,200)PULSACIONES/MIN.
if (irrint.gt.ifm*3/10.and.irrint.lt.ifm*2) then
if (irrint.gt.max) max=irrint
if (irrint.lt.min) min=irrint
isum=irrint+isum
icon=icon+1
end if
end do
if (icon.gt.2) irrint=(isum-max-min)/(icon-2)
return
end
C ----------------------------------------------------------------
C CREA UN FICHERO CON LOS QT
C ---------------------------------------------------------------
subroutine escribe_qt(tit,ifm,iqt,nqrs)
character*16 tit
dimension iqt(8000)
c CREA UN FICHERO NONDE GUARDA EL QT
open(unit=1,file=tit,status='new')
do i=1,nqrs
write(1,*) iqt(i)*1000/ifm
end do
close(unit=1)
return
end