C	 ======================================================================
C 	 punts.f                                                              =
C	 SUBROUTINES QUE ENS TROBEN TOTS ELS PUNTS SIGNIFICATIUS DEL          =
C      	 ECG A PARTIR DE LES POSICIONS DELS QRS I DEL SENYAL FILTRAT I DERIVAT=
c         Eudald Bogatell    (1-6-91)
c         David Vigo Anglada (  9-92)
C         Last modified: 24-2-2006 (by G. Moody, to avoid compiler warnings)
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       _______________________________________________________________________

         subroutine p_significatius(dbuf,jqrs,iqrs,if,is,ns,f,ipbeg, 
     &		ippos1, ipend, iqbeg, iqpos, ispos, isend, itbeg, 
     &		itpos, itend, irpos, ecgpb, derfi, irrpos, iqend,isbeg,
     &		rrmed, ecg, deri, itpos2, morf, idis, nep)
	 
c	 Aquesta subrutina es cridada desde ECGMAIN com una opcio, i ens
c	 troba tots els punts significatius de les ones Q,S,T, i P.

	 
	 dimension dbuf(100000),irpos(8000),ecgpb(100000),iqrs(8000)
	 dimension ipbeg(8000),ippos1(8000),ippos2(8000),ipend(8000),iqbeg(8000)
	 dimension iqpos(8000),ispos(8000),isend(8000), itbeg(8000),itpos(8000)
	 dimension itend(8000), derfi(100000), deri(100000), basel(8000)
	 dimension irrpos(8000), iqend(8000), isbeg(8000), ecg(100000)   
	 dimension itpos2(8000), morf(5,8000), idis(3,20),eauxs(100000)
	 dimension ecg_off(100000)
	 logical comp_qs, eonaq, eonas

c	 isf: mostra final a tractar
c	 dbuf: vector que conti el senyal filtrat i derivat
c	 ecgpb: senyal ECG filtrat a passa banda
c	 ecg: senyal ECG
c	 deri: senyal ECG derivat
c	 derfi: senyal dbuf filtrat
c 	 jqrs: nombre de qrs detectats i confirmats que tractarem
c 	 samp: interval de mostreig en ms
c	 iqrs: posicio dels complexes QRS detectats i confirmats
c	 irpos: vector que conti la posicio de l'ona R
c 	 ipbeg, ippos, ipend: inici, pic i final de l'ona P
c 	 iqbeg, iqpos, iqend: inici, pic i final de l'ona Q
c 	 isbeg, ispos, isend: inici, pic i final de l'ona S
c 	 itbeg, itpos, itpos2, itend: inici, pics i final de l'ona T
c 	 basel: conte les linies de base de cada QRS
c	 idis: conte els posibles episodis de VT i VF
c	 nep: nombre d'episodis de VT i VF
c	 comp_qs: es certa si es detecta un complexe QS
c	 eonaq: es certa si es detecta ona Q
c	 eonas: es certa si es detecta ona S




	 kr=5                        
c	 per definir el principi i final de l'ona R
	 isf=(is+ns)*if
	 samp=1000.0/if
c	 inicialitzem rrmed
	 rrmed=(iqrs(jqrs)-iqrs(1))/(jqrs-1)
	 n=1
         do i=1,300
		iqbeg(i)=0
		iqpos(i)=0
		iqend(i)=0
		irpos(i)=0
		irrpos(i)=0
		isbeg(i)=0
		ispos(i)=0
		isend(i)=0
		ipbeg(i)=0
		ippos1(i)=0
		ipend(i)=0
		itbeg(i)=0
		itpos(i)=0
		itpos2(i)=0
		itend(i)=0
	 end do

c	   Apliquem un filtre passa baix de 2 grau a la derivada per tal de
c 	   eliminar les components de alta frequencia i determinar millor 
c	   les ones P i T.
c	   write(6,*) n, 'antes de los filtros pb'
	   call fpb(isf, if, 40, ecgpb, eauxs, retard)
c           call fpb(isf, if, 40, ecg_off, eauxs, retard)
c	   retard=2*retard
           call normaliz_i(if,10,daux,eauxs)
           call der(isf, if, eauxs,derfi)
           call normaliz_i(if, 2, df_rmax, derfi)
c           write(6,*) if, df_rmax, derfi(1),derfi(100000)

	
	 jep=1	 
c	 QRS que estem tractant
 	  do while(n.le.jqrs.and.isf-iqrs(n).gt.500/samp)
c         afegim una proteccio en el cas de que el iqrs(n) sigui igual a 
c         cero; cas que es pot donar provenint de la subrutina multideriv, 
c         en les derivacions que s'hi ha detectat un fals negatiu.
          if (iqrs(n).ne.0) then
c	 si el qrs cau dins dun episodi de VT o VF no hi busquem cap complexe
c	 QRS
	   if (nep.gt.0.and.iqrs(n).ge.(idis(1,jep)+is)*if-if*0.3.and.
     &		jep.le.nep) then
	     if (iqrs(n).le.(idis(2,jep)+is)*if) then
	     	irpos(n)=iqrs(n)
		go to 100
	      else
		jep=jep+1
	     end if
	   end if	
		
c	   write(6,*) n, 'antes de onar'	
	   call onar(dbuf, ecgpb, iqrs, irpos, irrpos, iqbeg,iqpos, 
     &		ispos, isend, iqend, isbeg, n, samp, ecg,comp_qs,ymaxaux)
c	   write(6,*) n, 'despues de onar'
c	   busquem el pic anterior i posterior a la posicio de la ona R
	   call busca_pic2(irpos(n), dbuf, ima, 'i')
	   ymax=dbuf(ima)
	   call busca_pic2(irpos(n), dbuf, imi, 'd')
	   ymin=dbuf(imi)
c	 protegim per R amples on la derivada ti molts pics pero que la r
c	 sigui significativa
           if (.not.comp_qs) then
	   if (ymax.gt.ymaxaux/4) then
	   inicia=irpos(n)-nint(70/samp)
	   if (inicia.lt.0) inicia=0
	   call buscamaxmin(inicia,irpos(n),
     &	   	dbuf, iaux, yaux, imaa, ymaxa)
	   ilim=irpos(n)+nint(70/samp)
	  call buscamaxmin(irpos(n),ilim,dbuf, imia, ymina, iaux, yaux)
	   if (ymaxa.gt.ymax) then 
		ymax=ymaxa
		ima=imaa
	   end if
	   if (ymina.lt.ymin) then
c	 write(6,*) n, iqrs(n),irpos(n), ima,'primero'
 
		ymin=ymina
		imi=imia
	   end if
	   end if
	   if (abs(ymax).gt.abs(ymin)) then
	       	 dermax=abs(ymax)
	       else
	 	 dermax=abs(ymin)
	   end if

c	 i encara afegim unaltra proteccio si no tenim una Q o S molt marcada
c	 es a dir amb un pic de la derivada gran.
	    ilim=ima-nint(70/samp)
	    ilim2=ima-nint(30/samp)
	   if (ymax.gt.ymaxaux/4) then
	  call buscamaxmin(ilim,ilim2,dbuf, iaux, yaux, imaa, ymaxa)
	   if (ymaxa.gt.dermax/5) then
c	   write(6,*) n, iqrs(n),irpos(n), ima,'ultimo'

		ymax=ymaxa
		ima=imaa
 	   end if
            ilim=imi+nint(40/samp)
	    ilim2=imi+nint(100/samp)
 	   call buscamaxmin(ilim,ilim2,dbuf, imia, ymina, iaux, yaux)
	   if (ymina.lt.-1*dermax/5) then 
		ymin=ymina
		imi=imia
	   end if
	   end if
c	   write(6,*) n, iqrs(n),irpos(n), ima
	   else 
c	 tractem els complexes QS
	   inicia=irpos(n)-nint(150/samp)
	   if (inicia.lt.0) inicia=0
	   call buscamaxmin(inicia,irpos(n),
     &	   	dbuf,imi, ymin ,iaux, yaux)
	   ilim=irpos(n)+nint(180/samp)
	  call buscamaxmin(irpos(n),ilim,dbuf, iaux, yaux, ima, ymax)

c          write(6,2) imin,imax 
c 2 	   format('posicion QS=',i4,'isend=',i4 )
	
	   if (abs(ymax).gt.abs(ymin)) then
	       	 dermax=abs(ymax)
	       else
	 	 dermax=abs(ymin)
	   end if

c	 busquem el principi del QS amb el umbral de la derivada 
	   umbral=ymin/Kr
	   call creuar_umbral (dbuf, umbral, imi, iqbeg(n), 'i')
c          comprobem que el complexe QS no comengi amb un petit pic de pujada.
c          si fos aixi no l'hauriem detectat i l'inici trobat fora dolent
c          perque no el tenia amb compte. L'umbral ser` a Kr=3 i no 5 perque
c          is un petit pic ara.
	   ilim=iqbeg(n)-nint(35/samp)
           call buscamaxmin(ilim,   
     &                       iqbeg(n),dbuf,iaux,yaux,ima2,ymax2)
c	  write(6,*) ymax2, dermax
           if (ymax2.ge.dermax/30) then
               imi=ima2
               umbral=ymax2/2
    	       call creuar_umbral (dbuf, umbral, imi, iumb2, 'i')
	       if (iumb2.gt.iqbeg(n)-30/samp) iqbeg(n)=iumb2
           end if
	   if (-yaux.ge.dermax/30.and.iaux.lt.ima2) then
               imi=iaux
               umbral=yaux/2
    	       call creuar_umbral (dbuf, umbral, imi, iumb2, 'i')
	       if (iumb2.gt.iqbeg(n)-50/samp) iqbeg(n)=iumb2
            end if


c	 Buscamos el final del complejo QS
c	write(6,*) n,  ' QS'

	   umbral=ymax/kr
	   call creuar_umbral (dbuf, umbral, ima, isend(n), 'd')
	    ilim=irpos(n)+nint(180/samp)
 	   if (isend(n)-iqbeg(n).lt.80/samp) then
	    call buscamaxmin(irpos(n),ilim,
     &          dbuf,iaux,yaux,ima2,ymax2)
	    if (ymax2.gt.ymax) then
		umbral=ymax2/kr
		call creuar_umbral (dbuf, umbral, ima2, isend(n), 'd')
	    end if
	   end if

c          Ajuste para picos a la derecha de un QS 
           ilim=isend(n)+nint(20/samp)
           call buscamaxmin(isend(n),ilim,dbuf,imin2,ymin2,iaux,yaux)
c	   write(6,*) ymax2, dermax
           if (-ymin2.ge.dermax/20) then
               umbral=ymin2/2
    	       call creuar_umbral (dbuf, umbral, imin2, iumb2,'d')
	       if (iumb2.lt.isend(n)+30/samp) isend(n)=iumb2
           end if
	   

c           write(6,2) irpos(n),isend(n) 
c 2 	   format('posicion QS=',i4,'isend=',i4 )
	
	   end if


           eonaq=.true.
	   eonas=.true.
          
c	   write(6,*) n, 'antes de onaq'	
	   if (irrpos(n).ne.0.or.comp_qs) then
		go to 40
	   end if
 20	   call onaq (dbuf, dermax, ima, irpos, n, iqbeg, iqpos,
     &		 	 iqend, samp, ecg, deri, eonaq, itend,iqrs) 
c	   write(6,*) n, 'antes de onas'	
 30	   call onas (dbuf, dermax, imi, irpos, n, isend, ispos,
     &		      isbeg, samp, ecg, iqbeg, deri, eonas,iqrs)
c          write(6,*) n, 'antes de onap'	
 40	   call onap (n, derfi, isf, irpos, iqbeg, ipbeg, ippos1,
     &                ipend, samp, if, dermax, itend, ecgpb, ecg,iqrs)
c          write(6,*)  n, 'despues de onap'	
              
 	   call cal_base(n,ippos,ecg,basel,samp,ipend,iqbeg)     
c	 trobem el final de la Q i el principi de la S en el creuament de
c	 la linea de base.
	   
	   if (iqpos(n).ne.0.and..not.comp_qs.and.irrpos(n).eq.0) then
	      if (ecg(irpos(n)).gt.0) then
		 call creuar_umbral (ecg, basel(n), iqpos(n), iqend(n), 'd')
              else
		 iqend(n)=irpos(n)
	      end if
c	   si hem anat mes enlla de irpos no existeix la ona Q
	      if ((irpos(n)-iqend(n).lt.0).and.
     &            (iqrs(n).gt.iqpos(n)+10/samp) ) then 
		 iqbeg(n)=0
		 iqpos(n)=0
		 iqend(n)=0
		 eonaq=.false.
		 go to 20
	      end if
	   end if
	   if (ispos(n).ne.0.and..not.comp_qs.and.irrpos(n).eq.0) then
	   if (ecg(irpos(n)).gt.0) then
             call creuar_umbral (ecg, basel(n), ispos(n), isbeg(n), 'i')
	     else
              isbeg(n)=irpos(n)
	    end if

c	   si hem anat mes enlla de irpos no existeix la ona S
	   if ( (irpos(n)-isbeg(n).gt.0).and.
     &          (iqrs(n).lt.ispos(n)-10/samp) ) then 
		isbeg(n)=0
		ispos(n)=0
		isend(n)=0
		eonas=.false.
		go to 30
	   end if
	   end if
           call calcula_rmedio(rrmed, n, samp, iqrs)
	  
c            write(6,*) n, 'antes de onat'	
	   call onat (n, derfi, rrmed, isf, irpos, iqbeg, itbeg, itpos,
     &	   itpos2,itend, samp, if, isend, ecg, basel, morf,iqrs) 
c	   call onatold (n, derfi, rrmed, isf, irpos, iqbeg, itbeg, itpos,
c     &				itend, samp, if, isend,iqrs) 
	   if (irpos(n).ne.0) then 
	 	if(comp_qs) then
		   call test_pic (ecg, irpos(n), samp, 'v')
                 else
	           call test_pic (ecg, irpos(n), samp, 'p')
		end if
	   end if
          end if
 100	  n=n+1
	 end do

c	 omplim el reste dels vectors a cero
	 do i=jqrs+1,300
	    irpos(i)=0
	    ipbeg(i)=0
	    ippos1(i)=0
	    ippos2(i)=0
	    ipend(i)=0
	    iqbeg(i)=0
	    iqpos(i)=0
	    iqend(i)=0
	    isbeg(i)=0
	    ispos(i)=0
	    isend(i)=0
	    itbeg(i)=0
	    itpos(i)=0
	    itend(i)=0
	    irrpos(n)=0
	 end do

c	 a continuacio tenim unes ordres que permeten escriure en uns fitxers
c	 les mostres que es desitgin de tots els senyals de cara a treure'n

	 return
	 end


c-----------------------------------------------------------------------------
  
	 subroutine buscamaxmin (ib, ie, date, imi, ymin, ima, ymax)

c	 Ens treu per min i max els minim i maxim valors del senyal date
c	 entre les posicions ib i ie. Les posicions del min i del max ens
c	 les treu per imi i ima respectivament.

	 dimension date(100000)
	
         ima=ib
	 imi=ib
	 ymax=date(ib)
	 ymin=date(ib)
	 i=ib
	 do while (i.lt.ie)
	   if (date(i).ge.ymax) then
		ymax=date(i)
		ima=i
	   end if
	   if (date(i).le.ymin) then
		ymin=date(i)
		imi=i
	   end if
	   i=i+2
	 end do
	 if (date(imi-1).lt.ymin) then
	     ymin=date(imi-1)
	     imi=imi-1
	 end if
	 if (date(imi+1).lt.ymin) then
	     ymin=date(imi+1)
	     imi=imi+1
	 end if
	 if (date(ima-1).gt.ymax) then
	     ymax=date(ima-1)
	     ima=ima-1
	 end if
	 if (date(ima+1).gt.ymax) then
	     ymax=date(ima+1)
	     ima=ima+1
	 end if
	 return
	 end


c-------------------------------------------------------------------------


	 subroutine creuar_umbral (seny, umbral, inici, iumb, sentit)

c	 Ens treu per iumb la posicio del senyal 'seny' quan aquest
c 	 es inferior, en modul, al valor umbral, sempre buscan en un sentit
c	 dreta 'd' o esquerra 'i', a partir d'una posicio inicial inici

	 dimension seny(100000)
	 character*1 sentit

	 iumb=inici
	 if (sentit.eq.'d') then
	     if(seny(inici).gt.umbral) then
 		do while (seny(iumb).gt.umbral)
		  iumb=iumb+1
		end do    
                if (abs(seny(iumb-1)-umbral).lt.abs(seny(iumb)-umbral))
     &          iumb=iumb-1
	     else
 		do while (seny(iumb).lt.umbral)
		  iumb=iumb+1
		end do    
                if (abs(seny(iumb-1)-umbral).lt.abs(seny(iumb)-umbral))
     &          iumb=iumb-1
	     end if
	 else
	     if(seny(inici).gt.umbral) then
       		do while (seny(iumb).gt.umbral)
		  iumb=iumb-1
		end do    
                if (abs(seny(iumb+1)-umbral).lt.abs(seny(iumb)-umbral))
     &          iumb=iumb+1
	     else
       		do while (seny(iumb).lt.umbral)
		  iumb=iumb-1
		end do    
                if (abs(seny(iumb+1)-umbral).lt.abs(seny(iumb)-umbral))
     &          iumb=iumb+1
	     end if
	 end if
	 return
	 end

c------------------------------------------------------------------------------


	 subroutine busca_pic2 (inici, seny, ipic, sentit)

c	 Ens treu per ipic la posicio del primer maxim o minim relatiu 
c	 del senyal, trobat en sentit (d , i ), a partir de la posicio inici
c        Evita picos d euna sola muestra


	 dimension seny(100000)
	 character*1 sentit

	 if (sentit.eq.'d') then             
	        ipic=inici+1
 		do while ((seny(ipic).lt.seny(ipic+1).or.
     &			seny(ipic).lt.seny(ipic+2)).and.
     &			(seny(ipic).gt.seny(ipic-1).or.      
     &			seny(ipic).gt.seny(ipic-2)).or.
     &			(seny(ipic).gt.seny(ipic+1).or.  
     &			seny(ipic).gt.seny(ipic+2)).and.
     &			(seny(ipic).lt.seny(ipic-1).or.  
     &			seny(ipic).lt.seny(ipic-2)))
		  ipic=ipic+1
		end do    
	    else
		ipic=inici-1
 		do while ((seny(ipic).lt.seny(ipic+1).or.
     &			seny(ipic).lt.seny(ipic+2)).and.
     &			(seny(ipic).gt.seny(ipic-1).or.      
     &			seny(ipic).gt.seny(ipic-2)).or.
     &			(seny(ipic).gt.seny(ipic+1).or.  
     &			seny(ipic).gt.seny(ipic+2)).and.
     &			(seny(ipic).lt.seny(ipic-1).or.  
     &			seny(ipic).lt.seny(ipic-2)))
	 	ipic=ipic-1
		end do    
	 end if
	 return
	 end
c------------------------------------------------------------------------------


	 subroutine busca_pic1 (inici, seny, ipic, sentit)

c	 Ens treu per ipic la posicio del primer maxim o minim relatiu 
c	 del senyal, trobat en sentit (d , i ), a partir de la posicio inici


	 dimension seny(100000)
	 character*1 sentit

	 if (sentit.eq.'d') then             
	        ipic=inici+1
 		do while (seny(ipic).lt.seny(ipic+1).and.
     &			seny(ipic).gt.seny(ipic-1).or.      
     &			seny(ipic).gt.seny(ipic+1).and.  
     &			seny(ipic).lt.seny(ipic-1))  
		  ipic=ipic+1
		end do    
	    else
		ipic=inici-1
 		do while (seny(ipic).lt.seny(ipic+1).and.
     &			seny(ipic).gt.seny(ipic-1).or.      
     &			seny(ipic).gt.seny(ipic+1).and.  
     &			seny(ipic).lt.seny(ipic-1))
	 	ipic=ipic-1
		end do    
	 end if
	 return
	 end


c------------------------------------------------------------------------------

         subroutine onar(dbuf, ecgpb, iqrs, irpos, irrpos, iqbeg,iqpos, 
     &		ispos, isend, iqend, isbeg, n, samp, ecg, comp_qs,ymaxaux)

c	 Ens troba la posicio de la ona R en cas de que la deteccio del 
c	 complex QRS correspongui a una ona Q o S per ser aquestes anormalment
c 	 grans.
c	 Tambe ens distingeix els complexos RSR'

         dimension dbuf(100000), irpos(8000), iqrs(8000), ecgpb(100000)
	 dimension irrpos(8000), iqbeg(8000), isend(8000), ecg(100000)
         dimension iqpos(8000), ispos(8000), iqend(8000), isbeg(8000)
         logical comp_qs

c	 iqrs(n): posicio del QRS del qual buscarem la ona R
c	 irpos(n): posicio de la ona R del impuls que estem tractant 

	 comp_qs=.false.                                                   

c 	 Busquem els pics a dreta i esquerra de la posicio iqrs en dbuf
	 call busca_pic2(iqrs(n), dbuf, mpicd, 'd')
         call busca_pic2(iqrs(n), dbuf, mpici, 'i')
c 	 si el senyal en iqrs es menor que cero, o el pic de la derivada 
c	 a l'esquerra es > 0 i el pic de la derivada a la dreta es < 0,
c	 sent al mateix temps l'un molt mes gran que l'altre,
c	 llavors estem en els casos on iqrs no correspon a la R
	 ydbufi=dbuf(mpici)
	 ydbufd=dbuf(mpicd)
	 if(abs(ydbufi).gt.abs(ydbufd)) then
		ymaxaux=abs(ydbufi)
	    else
		ymaxaux=abs(ydbufd)
         end if
	 kpi=2
	 if (ecgpb(iqrs(n)).lt.0.or.ydbufi.gt.0.and.ydbufd.lt.0.and.
     &	     (kpi*ydbufi.lt.-1*ydbufd.or.kpi*(-1)*ydbufd.lt.ydbufi)) then
	 perc=0.25
	 if (ecgpb(iqrs(n)).gt.0.and.ydbufi.gt.0.and.ydbufd.lt.0
     &			.or.(1+perc)*(-1)*ydbufi.gt.ydbufd.and.
     &			(1-perc)*(-1)*ydbufi.lt.ydbufd) then
c	 estem en el cas RSR'
c	 write(6,*) n,  'RsR'

		perc=0.35
	 	kr=5
		krr=5
	 	if (ecgpb(iqrs(n)).lt.0) then
c	 iqrs correspon a la ona s, llavors la R' estara a la dreta i la R
c	 a la esquerra
		  ncero=mpicd
		  call detectar_cero(dbuf, ncero, 'd')
c	 mirem que el pendent de despres de R' sigui mes gran que un umbral
	       	  call busca_pic2(ncero, dbuf, mpda, 'd')
		  if (-1*dbuf(mpda).lt.ydbufd/2) go to 5
		  irrpos(n)=ncero
		  ispos(n)=iqrs(n)
	          ncero=mpici
		  call detectar_cero(dbuf, ncero, 'i')
c	      mirem que el pendent de antes de R sea mes gran que un umbral
	       	  call busca_pic2(ncero, dbuf, mpda, 'i')
c	 write(6,*) dbuf(mpda), mpda, ydbufi
		  if (dbuf(mpda).lt.-ydbufi/2) go to 5
 		  irpos(n)=ncero
		else if (abs(ydbufi).lt.abs(ydbufd)) then
c	 iqrs correspon a la ona R
		  ncero=mpicd
		  call detectar_cero(dbuf, ncero, 'd')
c 	 afegim una nova condicio verificadora
		  call busca_pic2(ncero, dbuf, mpic, 'd')
		  if (.not.((1+perc)*abs(ydbufd).gt.abs(dbuf(mpic)).and.
     &			(1-perc)*abs(ydbufd).lt.abs(dbuf(mpic)))) goto 10
                  ispos(n)=ncero
		  call detectar_cero(dbuf, ncero, 'd')
                  irrpos(n)=ncero
		  irpos(n)=iqrs(n)
		else if (abs(ydbufi).gt.abs(ydbufd)) then
c	 iqrs correspon a la ona R'
		  ncero=mpici
		  call detectar_cero(dbuf, ncero, 'i')
c 	 afegim una nova condicio verificadora
		  call busca_pic2(ncero, dbuf, mpic, 'i')
		  if (.not.((1+perc)*abs(ydbufi).gt.abs(dbuf(mpic)).and.
     &			(1-perc)*abs(ydbufi).lt.abs(dbuf(mpic)))) goto 10
		  ispos(n)=ncero
		  call detectar_cero(dbuf, ncero, 'i')
                  irpos(n)=ncero
		  irrpos(n)=iqrs(n)
                end if
c	 afegim una verificacio sobre les distancies R-R'
		if (abs(irpos(n)-irrpos(n)).gt.150/samp) then
		  if (ecgpb(iqrs(n)).gt.0) then
			goto 10
	            else
			goto 5
		  end if
		end if
c	 associem l'inici del RSR' a iqbeg i el final a isend, segons el
c	 criteri de la derivada
	        call busca_pic2(irrpos(n), dbuf, mpicd, 'd')
         	call busca_pic2(irpos(n), dbuf, mpici, 'i')
c	 fem un ajust dels pics del complexe QRS
	   if (irrpos(n).ne.0) then 
	   call test_pic (ecg, irrpos(n), samp, 'p')
	   end if
	   if (ispos(n).ne.0) then 
	   call test_pic (ecg, ispos(n), samp, 'v')
	   end if
	   if (irpos(n).ne.0) then 
	   call test_pic (ecg, irpos(n), samp, 'p')
	   end if
		umbral=dbuf(mpicd)/krr
		call creuar_umbral(dbuf, umbral, mpicd, isend(n), 'd')
	 	umbral=ecg(isend(n))
		call creuar_umbral(ecg, umbral, irrpos(n), isbeg(n), 'i')
		umbral=dbuf(mpici)/kr
		call creuar_umbral(dbuf, umbral, mpici, iqbeg(n), 'i')
		umbral=ecg(iqbeg(n))
		call creuar_umbral(ecg, umbral, irpos(n), iqend(n), 'd')
		iqpos(n)=0
  	    else
c 	 estem en el cas de Q o S anormalment grans, llavors
  5	 irrpos(n)=0	
c	 el pendent mes gran es el associat a la ona R, i,que tinguem
c	 presencia aprop de ona R per tant 
c	 busquem les dos posibles posicions de la ona R
		nrted=mpicd
		call detectar_cero(dbuf, nrted, 'd')
		nrtei=mpici
		call detectar_cero(dbuf, nrtei, 'i')
	 prr=1.4
c	write(6,*) n,iqrs(n),nrtei,nrted,dbuf(mpici),dbuf(mpicd)
	 if (abs(dbuf(mpicd)).gt.prr*abs(dbuf(mpici)).or.
     &		iqrs(n)-nrtei.gt.nrted-iqrs(n)) then         
c	 tenim ona Q en la posicio iqrs(n), i la ona R estara en el primer
c	 cero a la dreta
	    	 ncero=mpicd
		call detectar_cero(dbuf, ncero, 'd')
c	 mirem que el pendent de despres de R sigui mes gran que un umbral
	       	  call busca_pic2(ncero, dbuf, mpda, 'd')
c 	 si la distancia entre ones es massa gran potser no existeix ona R
c	 saltem a les sentencies de ona S per verificar-ho.
c	        write(6,*) n, iqrs(n), ncero, mpda,dbuf(mpda),ydbufd

	       if(ncero-iqrs(n).gt.150/samp.or.-1*dbuf(mpda).lt.ydbufd/10)
     &			 then
	 		go to 7
	 	   else
			irpos(n)=ncero
		end if
	   else

c	 tenim ona S en la posicio iqrs(n), i la ona R estara en el primer
c	 cero a l'esquerra
 7	  	ncero=mpici
		call detectar_cero(dbuf, ncero, 'i')
c	 mirem que el pendent de despres de R' sigui mes gran que un umbral
	         ilim=ncero-nint(60/samp)
	     call buscamaxmin(ilim,ncero,dbuf,iau,yau,imax2,ymax2)  
c 	 si la distancia entre ones es massa gran, o be el pic buscat mpdi
c	 no supera l'umbral , no existeix ona R, associem
c	 irpos a iqrs, despres no es detectara ona S, pero sabrem que ho es 
c	 per l'amplitut negativa.
c 	    write(6,*) n,ymax2, ydbufi,ecgpb(ncero),ecgpb(iqrs(n))
	  if ( (iqrs(n)-ncero.gt.140/samp.or.ymax2.lt.-1*ydbufi/10)
     &		.or.(ecgpb(ncero).lt.-abs(ecgpb(iqrs(n))*1/6)) )  then
			irpos(n)=iqrs(n)
			comp_qs =.true.
		   else
			irpos(n)=ncero
		end if
	 end if             
	 end if
	 else 
c	 La posicio iqrs correspon a la R, sent el QRS normal
 10	 irpos(n)=iqrs(n)
	 irrpos(n)=0
	 end if  
	 return
	 end
	 

c------------------------------------------------------------------------------

	 subroutine onaq(dbuf, dermax, ima, irpos, n, iqbeg, iqpos,
     &			iqend, samp, ecg , deri, eonaq, itend,iqrs)

c	 Ens treu per iqpos i iqbeg la posicio i l'inici de la ona Q, en cas
c	 de no existir ens treu la poisicio =0 i inici de la ona R.

	 dimension dbuf(100000), irpos(8000), iqbeg(8000), iqpos(8000), iqend(8000)
	 dimension ecg(100000), deri(100000), itend(100000),iqrs(8000)
	 logical eonaq

c	 eonaq: ens diu si existeix la ona Q
c	 irpos(n): posicio de la ona R del impuls que estem tractant 
c	 kq, kr: constants per definir els umbrals de les derivades per
c		 les ones Q i R respectivament
c	 dermax: derivada maxima del QRS
c	 ima: posicio del maxim anterior a la ona R en el senyal derivat

                 
c	 

	 kq=2
	 kr=5
	 inte=nint(10/samp)
c	 eonaq=.true.
c 	 busquem el cero previ a la posicio de la ona R partint del pic,
c	 que en principi correspondra a la posicio de la ona Q.
	 ncero=ima+1
	 nceau=ima
c	 busquem el cero en la derivada del ecg sense filtrar per tal de 
c	 partir de la posicio correcta per buscar el pic posterior
	 call detectar_cero (deri, ncero, 'i')
	 call detectar_cero (dbuf, nceau, 'i')
c	write(6,*) n, ncero, nceau, 'll'
c 	 si la distancia entre ones es superior a 60 ms , no existeix ona Q
	 if (irpos(n) - ncero.gt.80/samp.and.irpos(n).le.iqrs(n))
     &                eonaq=.false. 	
c	 si tenim ona Q, busquem el pendent maxim anterior

	 if (eonaq) then 
	 call busca_pic2 (nceau, dbuf, mpic, 'i')
c	 posem una proteccio per casos en que la derivada cuasi creua el 0.
	 call busca_pic2 (ima, dbuf, icep, 'i')  	 
c	 si el pic es gran, i no tenim un cuasi creuament de cero,
c	 treballarem amb dbuf
	 if(abs(dbuf(mpic)).gt.dermax/12.and..not.(icep.gt.mpic.and.
     &		abs(dbuf(icep)).lt.dermax/50)) then

c 	 probem de detectar be quan la ona P esta unida a la Q sense cap
c	 inflexio
	 if (irpos(n) - mpic.gt.90/samp.or.nceau-mpic.gt.30/samp
     &               .and.irpos(n).le.iqrs(n)) eonaq=.false. 	
c	 busquem el inici amb el umbral de la derivada o be amb el primer
c	 pic a l'esquerra
c	   umbral=dbuf(mpic)/kq
	   umbral=dbuf(mpic)/1.8
	   call creuar_umbral(dbuf, umbral, mpic, iumb, 'i')
	   call busca_pic1 (mpic, dbuf, ipic, 'i')
c	 write(6,*) n,irpos(n),nceau,mpic,iumb,ipic,dbuf(ipic),dbuf(mpic)


	   if(ipic.gt.iumb) iumb=ipic

c	   if (dbuf(ipic).gt.abs(dbuf(mpic)*0.8)) go to 10
c	   si la distancia es superior a 80 ms, no existeix ona Q
c 	 write(6,*) n, mpic, dbuf(mpic), irpos(n), dermax
 	   if (irpos(n)-iumb.gt.120/samp)  eonaq=.false.
	       
c          comprobem que el complexe Q no comengi amb un petit pic de pujada.
c          si fos aixi no l'hauriem detectat i l'inici trobat fora dolent
c          perque no el tenia amb compte. L'umbral ser` a Kr=3 i no 5 perque
c          is un petit pic ara.
	   ilimp=nint(30/samp)
           call buscamaxmin(iumb-ilimp,   
     &                       iumb,dbuf,imin2,ymin2,imax2,ymax2)
c	  write(6,*) iumb,imax2,ymax2, imin2,ymin2,dermax
          if (abs(ymin2).ge.dermax/20) then
                  umbral=ymin2/1.8
    	          call creuar_umbral (dbuf, umbral, imin2, iumb2, 'i')
	          if (iumb2.gt.iumb-40/samp) iumb=iumb2
          end if
c	  write(6,*) iumb, iumb2
	  if  ((imax2.lt.imin2.or.abs(ymin2).lt.dermax/20)
     &            .and.abs(ymax2).gt.dermax/20) then
                  umbral=ymax2/1.8
    	          call creuar_umbral (dbuf, umbral, imax2, iumb2, 'i')
	          if (iumb2.gt.iumb-40/samp) iumb=iumb2
          end if
  	 else
  10	 kq=3
c	 per tal de detectar ones Q de component frequencial molt elevat
c 	 fem servir la derivada del ecg sense filtrar.
c	 busquem el pic corresponent a deri.
	 call busca_pic2 (ncero, deri, mpic, 'i')
	 call test_pic (deri, mpic, samp, 'v')
c 	 si el pendent maxim de la ona Q es << que el pendent maxim llavors 
c	 no existeix ona Q
	 if (abs(deri(mpic)).lt.dermax/10
     &       .or.ncero-mpic.gt.30/samp) eonaq=.false.
c	 si existeix la ona Q busquem el seu principi amb el umbral de la
c	 derivada  

	 if (eonaq) then
	   umbral=deri(mpic)/kq
	   umbral=deri(mpic)/2.8
	   call creuar_umbral(deri, umbral, mpic, iumb, 'i')
c	   si la distancia es superior a 80 ms, no existeix ona Q
	   if (irpos(n)-iumb.gt.80/samp) eonaq=.false.
	 end if
	 end if
         end if
c	 si no existeix la ona Q busquem el principi de la ona R amb el
c	 umbral de la derivada 
	 if (.not.eonaq) then
c	   write(6,*) n,ima, dbuf(ima)
	   if (dbuf(ima).ge.1.5) kr=8
	   if (dbuf(ima).ge.3) kr=18
           if (dbuf(ima).ge.4.0) kr=28
	   umbral=dbuf(ima)/kr
           kr=5
	   call creuar_umbral (dbuf, umbral, ima, iumb, 'i')
c          asegurem-nos que la ona R no comenga amb un petit pic :
              ib=iumb-abs(30.0/samp)
              call buscamaxmin(ib,iumb,dbuf,iaux,yaux,ima2,ymax2)
              if (ymax2.ge.dermax/50) then
                  ima=ima2
                  umbral=ymax2/1.5
    	          call creuar_umbral (dbuf, umbral, ima, iumb2, 'i')
		  if (iumb2.gt.iumb-36/samp) iumb=iumb2
              end if
c	      write(6,*) n, ima, dbuf(ima), iumb,iumb2,ima2,ymax2
           end if
c	 omplim els vectors de sortida
	 if (eonaq) then
		iqpos(n)=ncero
            else
		iqpos(n)=0
	 end if
	 iqbeg(n)=iumb
 
c	 write(6,*) n, iqpos(n),iqbeg(n) , 'ippb'	
         if (n.gt.1.and.iqbeg(n).le.itend(n-1)) iqbeg(n)=itend(n-1)+1
	  
	   if (iqpos(n).ne.0) then 
	   call test_pic (ecg, iqpos(n), samp, 'v')
	   end if
c	 busquem el final de la ona Q
c	 if (eonaq) then
c	   umbral=ecg(iumb)
c	   call creuar_umbral(ecg, umbral, iqpos(n), iqend(n), 'd')
c	  else
c	   iqend(n)=0
c	 end if
	 return
	 end

	   
	   	 
c---------------------------------------------------------------------------   



	 subroutine onas(dbuf, dermax, imi, irpos, n, isend, ispos,
     &			isbeg, samp, ecg, iqbeg, deri, eonas,iqrs)

c	 Ens treu per ispos i isend la posicio i el final de la ona S, en cas
c	 de no existir ens treu la poisicio igual a 0 i el final de la ona R.

	 dimension dbuf(100000), irpos(8000), isend(8000), ispos(8000), isbeg(8000)
	 dimension derfi(100000), ecg(100000),iqbeg(8000), deri(100000), iqrs(8000)
	 logical eonas

c	 eonas: ens diu si existeix la ona S
c	 irpos(n): posicio de la ona R del impuls que estem tractant
c	 ks, kr: constants per definir els umbrals de les derivades per
c		 les ones S i R respectivament
c	 dermax: derivada maxima del QRS
c	 imi: posicio del maxim posterior a la ona R en el senyal derivat


c 	 busquem el cero posterior a la posicio de la ona R partint del minim,
c	 que en principi correspondra a la posicio de la ona S.

	 ks=3
	 kr=5
	 inte=nint(10/samp)
c	 eonas=.true.
	 ncero=imi
	 nceau=imi
c	 busquem el cero en la derivada del ecg sense filtrar per tal de 
c	 partir de la posicio correcta per buscar el pic posterior
	 call detectar_cero (deri, ncero, 'd')
	 call detectar_cero (dbuf, nceau, 'd')
	 
c 	 si la distancia entre ones es superior a 130 ms , no existeix ona S
	 if ((ncero-irpos(n).gt.130/samp).and.
     &          irpos(n).ge.iqrs(n))  eonas=.false. 
         if ((nceau.lt.iqrs(n)).and.(ecg(iqrs(n)).lt.0)) then
		ncero=iqrs(n)
		nceau=iqrs(n)
	 end if

c	 si tenim ona S, busquem el pendent maxim posterior
	 if (eonas.eqv..true.) then
c	 en comptes de buscar el primer pic a la dreta busquem el maxim
c	 en un interval pero amb un valor mes restringiu per si tenim 
c	 bloqueig de branca dreta
	  
	  if (ispos(n).eq.iqrs(n)) then
		ilim=nint(nceau+140./samp)
	 	    else 
	       ilim=nint(nceau+80./samp)   
          end if
	 call buscamaxmin (nceau,ilim , dbuf, iau, yau,mpic, ypic)
 
c	 write(6,*) n, nceau, mpic 
	 if (ypic.lt.dermax/10) call busca_pic2 (nceau, dbuf, mpic, 'd')  
c	 posem una proteccio per casos en que la derivada cuasi creua el 0.
	 call busca_pic2 (imi, dbuf, icep, 'd')  	 
c	 si el pic es gran, i no tenim un cuasi creuament de cero,
c	 treballarem amb dbuf
	 if( dbuf(mpic).gt.dermax/30.and.(.not.(icep.lt.mpic.and.
     &	    abs(dbuf(icep)).lt.dermax/50).or.(iqrs(n).eq.ncero)) ) then
c          posem ks en funcis de la pendent
           if (dbuf(mpic).ge.4.) ks=8
           if (dbuf(mpic).ge.4.75) ks=9
           if (dbuf(mpic).ge.6.2) ks=10 
c	   write(6,*) n, dbuf(mpic)
 	   umbral=dbuf(mpic)/ks
           ks=3
	   call creuar_umbral(dbuf, umbral, mpic, iumb, 'd')
	   call busca_pic1 (nint(mpic+10/samp), dbuf, ipic, 'd')
	   if(ipic.lt.iumb.and.dbuf(ipic).lt.dermax/15) iumb=ipic
c		write(6,1) n, ncero, iqrs(n), iumb
c 1	 	format('c=',i4,'ncero=',i4,'iqrs(n)=',i4,'iu=',i4)
 
c	   si la distancia es superior a 200 ms, no existeix ona S
	   if (iumb-irpos(n).gt.200/samp) then
             umbral=dbuf(mpic)/kr
             call creuar_umbral(dbuf, umbral, mpic, iumb, 'd')
             call busca_pic1 (nint(mpic+10/samp), dbuf, ipic, 'd')
             if(ipic.lt.iumb.and.dbuf(ipic).lt.dermax/3) iumb=ipic
           end if 
            
c           write(6,2) iumb-irpos(n) 
c 2 	   format('QRS=',i4 )
	
	 else
           call busca_pic2 (ncero, deri, mpic, 'd')
c 	 si el pendent maxim de la ona S es << que el pendent maxim llavors 
c	 no existeix ona S
c	 per tal de detectar ones S de component frequencial molt elevat
c 	 fem servir la derivada del ecg sense filtrar.
c	 busquem el pic corresponent a deri.
	   ks=3
	   call test_pic (deri, mpic, samp, 'p')
	   if (abs(deri(mpic)).lt.dermax/10.and.irpos(n).ge.iqrs(n))
     &                    eonas=.false.
c	 si existeix la ona S busquem el seu final amb el umbral de la
c	 derivada, o be, amb el primer pic a la dreta  
	   if (eonas) then
	   umbral=deri(mpic)/ks
	   call creuar_umbral(deri, umbral, mpic, iumb, 'd')
	   call busca_pic1 (nint(mpic+10/samp), dbuf, ipic, 'd')
	   if (ipic.lt.iumb) iumb=ipic
c	   si la distancia es superior a 80 ms, no existeix ona S
	   if (iumb-irpos(n).gt.200/samp)  then
             umbral=dbuf(mpic)/kr
             call creuar_umbral(dbuf, umbral, mpic, iumb, 'd')
             call busca_pic1 (nint(mpic+10/samp), dbuf, ipic, 'd')
             if(ipic.lt.iumb.and.dbuf(ipic).lt.dermax/3) iumb=ipic
           end if 
  	 end if
	 end if
         if ((iumb-irpos(n).gt.200/samp).and.
     &            irpos(n).ge.iqrs(n))  eonas=.false.
    
	 end if
c	 si no existeix la ona S busquem el final de la ona R amb el
c	 umbral de la derivada  o be amb el primer pic a la dreta
	 if (.not.eonas) then
	   umbral=dbuf(imi)/kr
	   call creuar_umbral (dbuf, umbral, imi, iumb, 'd')
	   call busca_pic2(nint(imi+10/samp), dbuf, ipic, 'd')
	   if(ipic.lt.iumb) iumb=ipic
	 end if
c	 omplim els vectors de sortida
	 if (eonas) then
		ispos(n)=ncero
            else
		ispos(n)=0
	 end if
	 isend(n)=iumb
	 if (isend(n).lt.iqrs(n)) isend(n)=iqrs(n)+1
	   if (ispos(n).ne.0) then 
	   call test_pic (ecg, ispos(n), samp, 'v')
	   end if
c	 busquem l'inici de la ona S
c	 if (eonas) then
c	   umbral=ecg(iumb) 
c	   call creuar_umbral(ecg, umbral, ispos(n), isbeg(n), 'i')
c	  else
c	   isbeg(n)=0
c	 end if
	 return
	 end
                                         
c--------------------------------------------------------------------------


	 subroutine onat (n, dbuf, rrmed, isf, irpos, iqbeg, itbeg, itpos,
     &   itpos2,itend, samp, if, isend, ecg, basel, morf,iqrs)

c	 ens localitza els punts inicial, pics i final de la ona T. Tambe ens
C	 treu la morfologia de l'ona T pel vector morf segons el codi:
c
c		morf(4,j) -ona T	0: existeix, normal
c					1:    "    , invertida
c					2:    "    , nomes pujada
c					3:    "    , nomes baixada
c					4: bifasica -+
c					5:    "     +-
c					6: no existeix


	 dimension dbuf(100000), irpos(8000), iqbeg(8000), itbeg(8000),itpos(8000)
	 dimension itend(8000), isend(8000), ecg(100000),back(6)
 	 dimension morf(5,8000), basel(8000), itpos2(8000), iqrs(8000)
	 logical flag1, inici_t, fivec
         

c 	 flag1: variable de control de bucle
c	 inici_t: ens diu si podem buscar el inici de la ona (T normal o
c	          invertda)
c	 bwind, ewind: amplada de la finestra
c	 ibw: posicio a partir de la que comengarem a buscar
c	 iew: posicio fins la que buscarem
c 	 ktb, kte: constants per definir els umbrals de la derivada per
c 	 	   localitzar l'inici i final de la ona T respectivament
c	 pco=valor de comparacio entre els pics
	 pco=8.0
	 ktb=2
	 kdis=nint(50/samp)
	 
c	Limite hasta donde se permitira ir al fin de la T 
	 if (rrmed.gt.900/samp) then
	         itqlim=nint(280/samp)
	         ewind=800
 	      else if (rrmed.gt.800/samp) then
	         itqlim=nint(250/samp)
 	         ewind=600
 	      else if (rrmed.gt.600/samp) then
		 itqlim=nint(200/samp)
		 ewind=450
 	      else 
		itqlim=nint(170/samp)
	        ewind=450
          end if
c	  write(6,*) itqlim	
c         factor de aumento de los umbrales si la primera vez va muy lejos
 
         do  i=1,6
          back(i)=1.0
	 end do
	 
         bwind=100
	          ibw=irpos(n)+nint(bwind/samp)
	 iew=irpos(n)+nint(ewind/samp)
	 flag1=.true.
	 inici_t=.true.

c	 inicialitzem finestres
 	
	 if (n.gt.1) then
          call inifinestres(rrmed, isf, irpos(n), bwind, ewind,  
     &				ibw, iew, samp, fivec)
         end if

c	 posem un reajust per QRS molt amples i per intervals rr anormalment
c	 petits.
	 if(ibw.le.isend(n)+kdis) then
		iew=iew+isend(n)-ibw+kdis
		ibw=isend(n)+kdis
	 end if   
 
	 if (iqrs(n+1).gt.0.and.iew.gt.iqrs(n+1)-210/samp) 
     &		iew=iqrs(n+1)-nint(210/samp)
	 if (iqrs(n+1).eq.0) iew=irpos(n)+nint(400/samp)
         if (fivec) return
c          write(6,5) n, (iqrs(n+1)-iew)*samp 
c 5 	   format('beat=',i4,' distancia iwe QRS', f6)
 	
c	 posem una altre condicio per detectar be T de pujada o baixada
	 call buscamaxmin (ibw, iew, dbuf, imi, ymin, ima, ymax) 
	 if (ymin.gt.0.or.ymax.lt.0) then
	  if (iew.eq.iqrs(n+1)-210/samp) then
	   if (ymin.gt.0) ymin=0
	   if (ymax.lt.0) ymax=0
	  else
	   do while (ymin.gt.0.or.ymax.lt.0.and.iqrs(n+1).gt.0.and.
     &		iew.lt.iqrs(n+1)-250/samp)
	     iew=iew+nint(25/samp)
	     call buscamaxmin (ibw, iew, dbuf, imi, ymin, ima, ymax) 
           end do
	  end if
	 end if
c	 inicialitzem morf al valor de no existencia per si no es pot clasificar
c	 dintre de cap dels patrons
	 morf(4,n)=6

	 do while (flag1.and.iew.gt.ibw)
c	 si el maxim i el minim tenen valors semblants busquem els pics
c	 anteriors i posteriors segons cada cas
	 kint1=nint(250/samp)
	 kint2=nint(300/samp)
c	 kend= valor en ms de marge minim a partir de isend on buscarem pics
	 kend=50
c	 ampmi= valor minim de les amplituts de les ones T bifasiques
         ampmi=0.075
c    	 kk es el nivel de comparacion para determinar ondas bifasicas
	 kk=3
	 com=0.3
c		write(6,*) n,imi, ymin, ima, ymax, ibw, iew, iqrs(n+1)
          if (-com*ymin.lt.ymax.and.-ymin.gt.com*ymax) then
c         write(6,2) n 
c 2 	   format('beat=',i4,' with max and min comparables')
	   if (imi.lt.ima) then
	     call buscamaxmin(ibw, imi, dbuf, iau, yau, imaa, ymaxa)
	     call buscamaxmin(ima, iew, dbuf, imip, yminp, iau, yau)
c	 si algun dels dos pics novament trobats es petit, no el considerem 
c	 com part de la ona T, o be si els dos a la vegada son grans
	     if (ymaxa.lt.ymax/kk.and.-yminp.lt.-ymin/kk.or.
     &		ymaxa.ge.ymax/kk.and.-yminp.ge.-ymin/kk) then
		morf(4,n)=1
c	        write(6,*) 'invertida primera opcion'
c		write(6,*) n, imi,ima,ymin,ymax

c	   write(6,*) n, ymaxa, ymin,ymax,yminp
	     else if(ymaxa.ge.ymax/kk) then
		go to 200
c	 mirem si estem en el cas de ona bifasica
c	 	ncea=imaa
c		call detectar_cero(dbuf, ncea, 'd')
c		if (ecg(ncea)-basel(n).lt.ampmi) then
c		   morf(4,n)=1
c		 else 
c	 tenim ona t bifasica +-
c		   imap=ima
c		   ymaxp=ymax
c		   morf(4,n)=5
c		end if
	     else if(-yminp.ge.-ymin/kk) then
		go to 200
c	 mirem si estem en el cas de ona bifasica
c	 	ncep=imip
c		call detectar_cero(dbuf, ncep, 'i')
c		if (ecg(ncep)-basel(n).lt.ampmi) then
c		morf(4,n)=1
c		 else 
c	 tenim ona t bifasica -+
c		   imia=imi
c		   ymina=ymin
c		morf(4,n)=4
c		end if
	     end if	
	    else
c	 estem en el cas de minim posterior al maxim
c	     if (ima-kint1.gt.isend(n)+kend/samp) then
c		inici=ima-kint1
c	      else
c		inici=isend(n)+nint(kend/samp)
c	     end if
c	     if (inici.gt.ima) inici=ima
c             if (imi+kint2.gt.iew) kint2=iew-imi
         call buscamaxmin(imi, iew, dbuf, iau, yau, imap, ymaxp)
	     call buscamaxmin(ibw, ima, dbuf, imia, ymina, iau, yau)
c	 si algun dels dos pics novament trobats es petit, no el considerem 
c	 com part de la ona T, o be si els dos a la vegada son grans
c	     write(6,*) n,n,imia,ymina,ima,ymax,imi,ymin,imap,ymaxp
	     if (ymaxp.lt.ymax/kk.and.-ymina.lt.-ymin/kk.or.
     &		ymaxp.ge.ymax/kk.and.-ymina.ge.-ymin/kk) then
		morf(4,n)=0
	     else if(ymaxp.ge.ymax/kk) then
c	 mirem si estem en el cas de ona bifasica
	 go to 200
c	 	ncep=imap
c		call detectar_cero(dbuf, ncep, 'i')
c		if (basel(n)-ecg(ncep).lt.ampmi) then
c		morf(4,n)=0
c		 else 
c	 tenim ona t bifasica +-
c		   imaa=ima
c		   ymaxa=ymax
c		   morf(4,n)=5
c	   	end if
	     else if(-ymina.ge.-ymin/kk) then
c	 mirem si estem en el cas de ona bifasica
	 go to 200
c	 	ncea=imia
c		call detectar_cero(dbuf, ncea, 'd')
c		if (basel(n)-ecg(ncea).lt.ampmi) then
c		morf(4,n)=0
c		 else 
c	 tenim ona t bifasica -+
c		   imip=imi
c		   yminp=ymin
c		morf(4,n)=4
c		end if
	     end if	
	   end if
 	 else 
c	 estem en el cas de pics diferents
c	 sino, si el maxim es mes gran que el minim busquem els dos minims
c	 anteriors i posteriors
 200	 continue
 	   if (ymax.gt.-1*ymin) then 
c	     if (ima-kint1.gt.isend(n)+kend/samp) then
c		inici=ima-kint1
c	      else
c		inici=isend(n)+nint(kend/samp)
c	     end if
c	     if (inici.gt.ima) inici=ima
	     call buscamaxmin(ibw, ima, dbuf, imia, ymina, iau, yau)
c	     if (ima+kint2.gt.iew) kint2=iew-ima
             call buscamaxmin(ima,iew,dbuf,imip,yminp, iau, yau)
             ncea=imia
	     if (ymina.lt.0) call detectar_cero(dbuf, ncea, 'd')
	     ampa=basel(n)-ecg(ncea)
             ncep=imip
	     if (yminp.lt.0) call detectar_cero(dbuf, ncep, 'i')
	     ampp=ecg(ncep)-basel(n)     
c	 si la primera amplitud es gran podem tenir T de pujada, bifasica -+ o
c	 invertida
c	     write(6,*) n, ampa, ampp, imia,ymina,imip,yminp
	     if (ampa+ampp.gt.ampmi)  then
	     if (-ymina.lt.ymax/pco.and.-yminp.lt.ymax/pco)  then
c	 tenim ona T de nomes pujada
		morf(4,n)=2      
	     else if (-ymina.ge.ymax/pco.and.-yminp.ge.ymax/pco) then
c	 tenim ona T bifasica -+
		morf(4,n)=4
	     else if (-ymina.ge.ymax/pco.and.-yminp.lt.ymax/pco) then
c                    tenim ona T invertida
c 	  write(6,*) 'invertida segunda opcion'
 		     morf(4,n)=1
		     ymin=ymina
	 	     imi=imia
c	             write(6,*) imi,ima,ymin,ymax
            else if (-ymina.lt.ymax/pco.and.-yminp.ge.ymax/pco) then
c	             tenim ona T normal
		     morf(4,n)=0
		     ymin=yminp
	 	     imi=imip
	    end if
	   end if      		

	 else if (ymax.lt.-1*ymin) then
c	     if (imi-kint1.gt.isend(n)+kend/samp) then
c		inici=imi-kint1
c	      else
c		inici=isend(n)+nint(kend/samp)
c	     end if
c	     if (inici.gt.imi) inici=imi
	     call buscamaxmin(ibw, imi, dbuf, iau, yau, imaa, ymaxa)
c	     if (imi+kint2.gt.iew) kint2=iew-imi
             call buscamaxmin(imi, iew,dbuf,iau,yau,imap,ymaxp)
             ncea=imaa
	     if (ymaxa.gt.0) call detectar_cero(dbuf, ncea, 'd')
	     ampa=ecg(ncea)-basel(n)
             ncep=imap
	     if (ymaxp.gt.0) call detectar_cero(dbuf, ncep, 'i')
	     ampp=basel(n)-ecg(ncep)
c	 si la primera amplitut es gran podem tenir ona T de baixada, normal
c	 o be bifasica +-
c	 write(6,*) n,imaa,ymaxa,imi,ymin,imap,ymaxp,ampa,ampp
	     if (ampa+ampp.gt.ampmi) then
	     if (ymaxa.lt.-ymin/pco.and.ymaxp.lt.-ymin/pco) then
c	 tenim ona T de nomes baixada
c	        write(6,*) n, 'bajada'
		morf(4,n)=3      
	     else if (ymaxa.ge.-ymin/pco.and.ymaxp.ge.-ymin/pco) then
c	 tenim ona T bifasica +-
c              write(6,*) n, 'bifasica'
 		morf(4,n)=5
	     else if (ymaxa.ge.-ymin/pco.and.ymaxp.lt.-ymin/pco) then
c	 tenim ona T normal
c		write(6,*) n, 'normal'
 		morf(4,n)=0
		ymax=ymaxa
		ima=imaa
	    else if (ymaxa.lt.-ymin/pco.and.ymaxp.ge.-ymin/pco) then 
c	 tenim ona T invertida
		morf(4,n)=1
	        ymax=ymaxp
		ima=imap
c  		write(6,*) imi,ima,ymin,ymax
	     end if 
	    end if
	   end if
	 end if

c	 procedim a buscar pics, inicis i finals segons cada cas
 300	 go to (10, 20, 30, 40, 50, 60, 70), morf(4,n)+1
c 	 tenim ona T normal 
 10 	 umba=ymax/ktb	
	 call creuar_umbral (dbuf, umba, ima, iumba, 'i')
c 	 si hem anat mes enrera de isend, trobem el principi de la T
c	 amb un criteri del primer pic a l'esquerra
	 if (iumba.lt.isend(n)) call busca_pic2(ima, dbuf, iumba, 'e')
c 	 si encara hem anat massa enrera l'associem a isend
	 if (iumba.le.isend(n)) iumba=isend(n)+2 
	 itbeg(n)=iumba

C        kte valdr` 4,5,6,7 en funcis del valor de la derivada max.
         if (abs(ymin).ge.0.41) then
             kte=7
         else if (abs(ymin).ge.0.35) then
             kte=6
         else if (abs(ymin).ge.0.25) then
             kte=5
         else if (abs(ymin).ge.0.10) then   
             kte=4
	 else if (abs(ymin).lt.0.10) then   
             kte=3.5
         end if

         if (kte/back(1).ge.1.1) then
	      umbp=ymin*back(1)/kte
	 else 
	   umbp=ymin/1.1
	 end if
	 call creuar_umbral (dbuf, umbp, imi, iumbp, 'd')
	 itend(n)=iumbp
 
	 icero1=imi
	 call detectar_cero (dbuf, icero1, 'i')
	 icero2=ima
	 call detectar_cero (dbuf, icero2, 'd')
	 if (icero1.ne.icero2) then
		icero=(icero1+icero2)/2
	 else 
            icero=icero1
         end if
	 if (icero.ge.itend(n).or.icero.le.itbeg(n))
     &           icero=itbeg(n)+(itend(n)-itbeg(n))/2

	 itpos(n)=icero
	 itpos2(n)=icero
	 go to 100

c	 tenim ona T invertida
 20 	 umba=ymin/ktb	
	 call creuar_umbral (dbuf, umba, imi, iumba, 'i')
c 	 si hem anat mes enrera de isend, trobem el principi de la T
c	 amb un criteri del primer pic a l'esquerra
	 if (iumba.lt.isend(n)) call busca_pic2(imi, dbuf, iumba, 'e')
c 	 si encara hem anat massa enrera l'associem a isend
	 if (iumba.le.isend(n)) iumba=isend(n)+2 
	 itbeg(n)=iumba

C        kte valdr` 4,5,6,7 en fuincis del valor de la derivada max.
         if (abs(ymax).ge.0.41) then
             kte=7
         else if (abs(ymax).ge.0.35) then
             kte=6
         else if (abs(ymax).ge.0.25) then
             kte=5
         else if (abs(ymax).ge.0.10) then   
             kte=4
	 else if (abs(ymax).lt.0.10) then   
             kte=3.5
         end if

 
         if (kte/back(2).gt.1.1) then
	      umbp=ymax*back(2)/kte
	 else 
	   umbp=ymax/1.1
	 end if
 	 call creuar_umbral (dbuf, umbp, ima, iumbp, 'd')
	 itend(n)=iumbp
c        write(6,*) n, imi, ima, iumbp, itbeg(n), itend(n),ymax,umbp,kte   
	 icero1=ima
	 call detectar_cero (dbuf, icero1, 'i')
	 icero2=imi
	 call detectar_cero (dbuf, icero2, 'd')
	 if (icero1.ne.icero2) then 
              icero=(icero1+icero2)/2
         else 
            icero=icero1
         end if
	  if (icero.ge.itend(n).or.icero.le.itbeg(n))
     &           icero=itbeg(n)+(itend(n)-itbeg(n))/2

c 	 write(6,1) icero1, icero2, icero
c 1	 format('icero1=',i4, ' icero1=',i4,' icero1=',i4)
      

	 itpos(n)=icero
	 itpos2(n)=icero
	 go to 100

c	 tenim ona T de nomes pujada

 30      continue
C        kte valdr` 4,5,6,7 en fuincis del valor de la derivada max.
         if (ymax.ge.0.41) then
             kte=7
         else if (ymax.ge.0.30) then
             kte=6
         else if (ymax.ge.0.20) then
             kte=5
         else if (ymax.gt.0.10) then   
             kte=4
	 else if (ymax.le.0.10) then   
             kte=3.5
         end if
   
         if (kte/back(3).ge.1.1) then
	      umbp=ymax*back(3)/kte
	 else 
	   umbp=ymax/1.1
	 end if
     

         call creuar_umbral (dbuf, umbp, ima, iumbp, 'd')
	 itend(n)=iumbp
	 icero=ima
	 call detectar_cero (dbuf, icero, 'i')
	 call busca_pic2(ima-kdis, dbuf, ipic, 'i')
	 if (ipic.gt.icero) then
	    itpos(n)=ipic
	  else 
	    itpos(n)=icero
	 end if
	 if (itpos(n).le.isend(n)) itpos(n)=isend(n)+1
	 itpos2(n)=0
	 itbeg(n)=0
	 go to 100

c	 tenim ona T de nomes baixada   

 40      continue
C        kte valdr` 4,5,6,7 en funcis del valor de la derivada max.
          if (abs(ymin).ge.0.41) then
             kte=7
         else if (abs(ymin).ge.0.30) then
             kte=6
         else if (abs(ymin).ge.0.20) then
             kte=5
         else if (abs(ymin).ge.0.10) then   
             kte=4
	 else if (abs(ymin).lt.0.10) then   
             kte=3.5
         end if

         if (kte/back(4).ge.1.1) then
	      umbp=ymin*back(4)/kte
	 else 
	   umbp=ymin/1.1
	 end if
     
	 call creuar_umbral (dbuf, umbp, imi, iumbp, 'd')
	 itend(n)=iumbp
	 icero=imi
	 call detectar_cero (dbuf, icero, 'i')
	 call busca_pic2(imi-kdis, dbuf, ipic, 'i')
	 if (ipic.gt.icero) then
	    itpos2(n)=ipic
	  else 
	    itpos2(n)=icero
	 end if
	 if (itpos2(n).le.isend(n)) itpos2(n)=isend(n)+1

	 itpos(n)=0
	 itbeg(n)=0
         go to 100

c	 tenim ona T bifasica -+

 50	 umba=ymina/ktb
	 call creuar_umbral(dbuf, umba, imia, iumba, 'i')
	 if (iumba.lt.isend(n)) call busca_pic2(imia, dbuf, iumba, 'e')
c 	 si encara hem anat massa enrera l'associem a isend
	 if (iumba.le.isend(n)) iumba=isend(n)+2 
	 itbeg(n)=iumba

C        kte valdr` 4,5,6,7 en funcis del valor de la derivada max.
          if (abs(yminp).ge.0.41) then
             kte=7
         else if (abs(yminp).ge.0.30) then
             kte=6
         else if (abs(yminp).ge.0.20) then
             kte=5
         else if (abs(yminp).ge.0.10) then   
             kte=4
	 else if (abs(yminp).lt.0.10) then   
             kte=3.5
         end if
        
         if (kte/back(5).ge.1.1) then
	      umbp=yminp*back(5)/kte
	 else 
	   umbp=yminp/1.1
	 end if

	  
	 call creuar_umbral(dbuf, umbp, imip, iumbp, 'd')
	 itend(n)=iumbp
	 icero=imia
	 call detectar_cero (dbuf, icero, 'd')
	 itpos2(n)=icero
	 icero=imip
	 call detectar_cero (dbuf, icero, 'i')
	 itpos(n)=icero
         if (itpos(n).lt.itpos2(n)) itpos(n)=itpos2(n)
	 go to 100

c	 tenim ona T bifasica +-
 60	 umba=ymaxa/ktb
	 call creuar_umbral(dbuf, umba, imaa, iumba, 'i')
	 if (iumba.lt.isend(n)) call busca_pic2(imaa, dbuf, iumba, 'e')
c 	 si encara hem anat massa enrera l'associem a isend
	 if (iumba.le.isend(n)) iumba=isend(n)+2 
	 itbeg(n)=iumba

C        kte valdr` 4,5,6,7 en funcis del valor de la derivada max.
         if (abs(ymaxp).ge.0.41) then
             kte=7
         else if (abs(ymaxp).ge.0.30) then
             kte=6
         else if (abs(ymaxp).ge.0.20) then
             kte=5
         else if (abs(ymaxp).ge.0.10) then   
             kte=4
	 else if (abs(ymaxp).lt.0.10) then   
             kte=3.5
         end if
         
         if (kte/back(6).ge.1.1) then
	      umbp=ymaxp*back(6)/kte
	 else 
	   umbp=ymaxp/1.1
	 end if
 	
         call creuar_umbral(dbuf, umbp, imap, iumbp, 'd')
	 itend(n)=iumbp
	 icero=imaa
	 call detectar_cero (dbuf, icero, 'd')
	 itpos2(n)=icero
	 icero=imap
	 call detectar_cero (dbuf, icero, 'i')
	 itpos(n)=icero
	 if (itpos(n).lt.itpos2(n)) itpos(n)=itpos2(n)
         go to 100

c	 no tenim ona t
 70	 itbeg(n)=0
	 itpos(n)=0
	 itpos2(n)=0
	 itend(n)=0
	 go to 100

c 	 apliquem una regla de validesa respecte el interval QT
         
c 100 	    write(6,3) n,  itend(n)-iqbeg(n), 750/samp
c 3 	    format('beat=',i2,' QT=', i4, '  QT limite=', f4)
 	  
 100        if (itend(n).ge.isf) go to 70
c	    write(6,*)n, 'isf=', isf, itend(n)
	    if (itend(n)-iqbeg(n).lt.950/samp.and.
     &       (iqrs(n+1)-itend(n).gt.itqlim.or.iqrs(n+1).eq.0.or.
     &        itend(n)-iqbeg(n).lt.400/samp) ) then
		flag1=.false.
	      else
 		if (iew.gt.ibw+100/samp) then
	            iew=iew-nint(50/samp)
		  else
	            iew=iew-nint(25/samp)
	        end if
                call buscamaxmin (ibw, iew, dbuf, imi, ymin, ima, ymax)
	        back(morf(4,n)+1)=back(morf(4,n)+1)*1.8
 	      end if
c	  write(6,*) n,itend(n), iqrs(n+1),itend(n)-iqbeg(n)
  	 end do
        
         if (itend(n).gt.iqrs(n+1)-100/samp.and.iqrs(n+1).ne.0) then
           itend(n)=0
	   itpos(n)=0
	   itbeg(n)=0
	   itpos2(n)=0
	 end if
     
	 return
	 end



c------------------------------------------------------------------------------
	 subroutine onatold (n, dbuf, rrmed, isf, irpos, iqbeg, itbeg, itpos,
     &				itend, samp, if, isend,iqrs)          

c	 subrutina vella on no esta implementat la posiblitat de ones T
c	 bifasiques.
c	 ens localitza els punts inicial, pic i final de la ona T

	 dimension dbuf(100000), irpos(8000), iqbeg(8000), itbeg(8000),itpos(8000)
	 dimension itend(8000), isend(8000),iqrs(8000)
	 logical flag1, inici_t, fivec

c 	 flag1: variable de control de bucle
c	 inici_t: ens diu si podem buscar el inici de la ona (T normal o
c	          invertda)
c	 bwind, ewind: amplada de la finestra
c	 ibw: posicio a partir de la que comengarem a buscar
c	 iew: posicio fins la que buscarem
c 	 ktb, kte: constants per definir els umbrals de la derivada per
c 	 	   localitzar l'inici i final de la ona T respectivament
	 ktb=3
	 kte=2
	 kdis=nint(20/samp)
	 bwind=100
	 ewind=450
	 ibw=irpos(n)+nint(bwind/samp)
	 iew=irpos(n)+nint(ewind/samp)
c	 filtrem passa-baix la derivada del senyal compresa entre les finestres
c	 call filtre_pb (ibw, iew-ibw, if, 30, dbuf,retard)
	 flag1=.true.
	 inici_t=.true.
c	 inicialitzem finestres
	 call inifinestres(rrmed, isf, irpos(n), bwind, ewind, ibw, iew, 
     &				samp, fivec)
c	 posem un reajust per QRS molt amples i per intervals rr anormalment
c	 petits.
	 if(ibw.le.isend(n)+kdis) then
		iew=iew+isend(n)-ibw+kdis
		ibw=isend(n)+kdis
	 end if
	 if (iqrs(n+1).gt.0.and.iew.gt.iqrs(n+1)-100/samp) 
     &		iew=iqrs(n+1)-100/samp
         if (fivec) return
c	 posem una altre condicio per detectar be T de pujada o baixada
	 call buscamaxmin (ibw, iew, dbuf, imi, ymin, ima, ymax) 
	 if (ymin.gt.0.or.ymax.lt.0) then
	   do while (ymin.gt.0.or.ymax.lt.0)
	     iew=iew+25/samp
	     call buscamaxmin (ibw, iew, dbuf, imi, ymin, ima, ymax) 
           end do
	   iew=iew+100/samp
	 end if
	 if (iqrs(n+1).gt.0.and.iew.gt.iqrs(n+1)-100/samp)
     &		iew=iqrs(n+1)-100/samp
	 do while (flag1.and.iew.gt.ibw)
	   ieaux=nint(iew-50/samp)
	   call buscamaxmin (ibw, ieaux, dbuf, imi, ymin, ima, ymax) 
	   imig=imi
	   yming=ymin
	   if (imi.lt.ima) 
     &		call buscamaxmin (ima, iew, dbuf, imi, ymin, iaux1, aux2)
c	 regla de decisio per clasificar la ona T com a normal, invertida, 
c 	 depujada o de baixada
	   if (imi.ne.imig.and.ymin.gt.0.or.
     &		abs(ymin).lt.ymax/3.or.
     &		abs(yming)/4.gt.ymax.and.abs(yming)/4.gt.abs(ymin)) then
	     if (abs(yming).gt.3*ymax.or.ima.gt.iew-50/samp) then
c 	 tenim ona T de nomes baixada
	         imi=imig
	         ymin=yming
		 inici_t=.false.
	       else
c	 tenim ona T invertida o be de pujada
	   if (imi.eq.imig.and.ymax.gt.3*abs(ymin)) then
c	 tenim ona t de nomes pujada
		ymin=ymax
	 	imi=ima
		inici_t=.false.
	     else
c	 tenim ona T invertida
		 imi=ima
 	         ymin=ymax
		 ima=imig
	 	 ymax=yming
	   end if
	     end if
	   end if
c 	 altrament tenim ona T normal
	   icero=imi
c	 busquem el cero anterior al que en cada cas correspongui a la posicio
c	 imi, i que en principi sera el pic de la ona T
	   call detectar_cero (dbuf, icero, 'i')
c	 si hem anat mes enrera de isend associem el pic a la primera vall 
c	 a la esquerra
   	 if (icero.lt.isend(n)) call busca_pic2(imi, dbuf,icero, 'e')
c	 buquem el final de la ona T definint un umbral de la derivada
	   umbral=ymin/kte
	   call creuar_umbral (dbuf, umbral, imi, iumb, 'd')
c 	 apliquem una regla de validesa respecte el interval QT
	   if (iumb-iqbeg(n).lt.620/samp) then
		flag1=.false.
		itpos(n)=icero
		itend(n)=iumb
	     else
		iew=iew-50/samp
	   end if
	 end do
	 itpos(n)=icero
c	 busquem el inici de la ona T mitjangant un umbral
	 if(inici_t) then
		umbral=ymax/ktb
		call creuar_umbral(dbuf, umbral, ima, iumb, 'i')
c 	 si hem anat mes enrera de isend, trobem el principi de la T
c	 amb un criteri del primer pic a l'esquerra
		if (iumb.lt.isend(n)) call busca_pic2(ima, dbuf, iumb, 'e')
c 	 si encara hem anat massa enrera l'associem a isend
		if (iumb.lt.isend(n)) iumb=isend(n)
		itbeg(n)=iumb
	    else
		itbeg(n)=0
	 endif
	 return
	 end



c------------------------------------------------------------------------------


	 subroutine inifinestres (rrmed, isf, irp, bwind, ewind, ibw,
     &					 iew,samp, fivec)

c	 inicialitza les finestres a on es buscara l'ona T.

	 logical fivec
         fivec=.false.
	 if (irp+ewind/samp.lt.isf)then
		if (rrmed.lt.750/samp) then
c			ibw=nint(100/samp)+irp
			ibw=nint(125/samp)+irp
			iew=nint(rrmed*0.6)+irp
		else
			ibw=nint(bwind/samp)+irp
			iew=nint(ewind/samp)+irp
		end if
	    else
		if (irp+bwind/samp.lt.isf) then
			ibw=nint(bwind/samp)+irp
			iew=isf
		else
			fivec=.true.
		endif
	 end if
	 return
	 end				 
c ---------------------------------------------------------------------------

         subroutine filtre_pb (ni,n,if,ifc,x,retard)
         dimension x(100000),y(400)

c        APLICA A LA SEQAL QUE ENTRA EN EL VECTOR x UN FILTRO PASO BAJO. lA 
C        SEQAL QUE ENTRA LA SUPONGO MUESTREADA A if hZ.
C         EL FILTRO TIENE UN POLO DOBLE EN (0,1) Y ifc CEROS DOBLES EN
C        EL CIRCULO UNIDAD .LA GANANCIA DEL
C        FILTRO SERIA DEl NUMERO de CEROS AL CUADRADO. 
C         EL RETARDO DEL FILTRO ES DE NUMERO de CEROS MENOS UNO
C         LAS CONDICIONES INICIALES SUPONGO QUE LA SEQAL ES CONSTANTE ANTES DEL
C        INSTANTE INICIAL CON VALOR 0
C         lA ECUACION SERIA Y(nT)=2Y(nT-T)-Y(nT-2T)+X(nT)-2X(nT-cerosT)
c                                 +X(nT-2*cerosT).
C         lA SEQAL FILTRADA LA DEVUELVE EN EL VECTOR y


C        RM ES EL VALOR MEDIO DE LOS 10 PRIMEROS PUNTOS
c         sum=x(ni)
c         do i=ni+2,ni+10
c          sum=sum+x(i)
c         end do
         rm=x(ni)
       
c        l ES EL NUMERO DE CEROS
C        ifc ES LA FRECUENCIA A LA QUE EL FILTRO CAE A CERO
C        if ES LA FRECUENCIA A LA QUE ESTA MUESTREADA LA SEQAL 
C        ni ES LA POSICION INICIAL A PARTIR DE LA CUAL QUIERO FILTRAR
c        n ES EL NUMERO DE PUNTOS A FILTRAR
C        SE DEBE CUMPLIR (n<1000-l) SI NO SE SALDRA DE LOS LIMITES DEL VECTOR
         l=nint(1.*if/ifc)
         y(1)=x(ni+1)+(l**2-1)*rm
         y(2)=(2*y(1)+x(ni+2)-(l**2+1)*rm)
         do i=3,l
             y(i)=(2*y(i-1)-y(i-2)+x(ni+i)-rm)
         end do
         do i=l+1,2*l
             y(i)=(2*y(i-1)-y(i-2)+x(ni+i)-2*x(ni+i-l)+rm)
         end do
         do i=2*l+1,n+l
             y(i)=(2*y(i-1)-y(i-2)+x(ni+i)-2*x(ni+i-l)+x(ni+i-2*l))
         end do
 
C        QUITO EL RETARDO DE LA SEQAL
         retard=l
         if (n.le.(400-l+1)) then
                              do i=1,n
                                x(ni+i)=y(i+l)
                              end do
                             else
                              do i=1,400-l+1
                                x(ni+i)=y(i+l)
                              end do
c                             do i=400-l+1,n
c                               y(i)=0
c                             end do
         end if 

C        PONGO EL RESTO A CERO
         do i=n+1,400
            y(i)=0
         end do
 
         return
         end

C        ---------------------------------------------------------------


	 subroutine onap (n, dbuf, isf, irpos, iqbeg, ipbeg, ippos1,
     &	   ipend, samp, if, dermax, itend, ecgpb,ecg,iqrs)
                         
c	 Aquesta subroutina ens busca inici i final de la ona P. Considera
c	 tant ones P normals com invertides.

	 dimension dbuf(100000), irpos(8000), iqbeg(8000), ipbeg(8000), ipend(8000)
	 dimension ippos1(8000), ippos2(8000), itend(8000) ,ecgpb(100000)
         dimension ecg(100000), iqrs(8000)
	 logical nofi

c	 bwind, ewind: amplada de la finestra
c	 ibw: posicio a partir de la que comengarem a buscar
c	 iew: posicio fins la que buscarem
c 	 kpb, kpe: constants per definir els umbrals de la derivada per
c 	 	   localitzar l'inici i final de la ona P respectivament
	 counter=0
	 nofi=.true.
	 kpb=2
	 kpe=2
	 bwind=200
c	 bwind=380
	 ewind=30
	 ippos1(n)=0
	 iew=iqbeg(n)-nint(ewind/samp)
	 ibw=iqbeg(n)-nint(bwind/samp)
	 if (ibw.lt.0) ibw=0
	 if (iew.lt.0) iew=0

         do while (nofi.and.ippos1(n).eq.0.and.iqbeg(n)-ibw.lt.300/samp)
c 	 no busquem mes enlla del final de T anterior, i a una distancia minima
c	 de qbeg (Vigo : he canviat 330 per 400)
	    if (n.eq.1.or.itend(n-1).eq.0) then
	       nofi=.false.
	    else
	       if (ibw.lt.itend(n-1)) then
		  ibw=itend(n-1)
		  nofi=.false.
	       else
		  nofi=.true.
	       end if
	    end if
	    kdis=nint(20/samp)
c              Estaba cuando los limites se amrcaban desde el QRS
c	 if(iqbeg(n)-iew.lt.kdis) then
c 	   iew=iqbeg(n)-kdis
c	 end if
	    if (iew.lt.0) iew=0

c 	 busquem el maxim i el minim dins la finestra
 10	    call buscamaxmin (ibw, iew, dbuf, imi, ymin, ima, ymax)
C---------------------------------------------------------------------------
C        Aixs is el criteri de l'Audal:
C 	 si les pendents son molt inferiors a les del QRS o be 
c	 si la posicio del minim esmenor que la del maxim, no tenim ona P
c	 if (imi.le.ima.or.ymax.lt.dermax/50.or.abs(ymin).lt.dermax/50) then

c	 tambe detectem ones P invertides, pertant nomes si els pendents
c	 son molt inferiors als del QRS, sent mes exigents en les invertides 
c	 per tal de no detectar-ne de falses.
C---------------------------------------------------------------------------
c   Vigo: canvio 70 per 30 per detectar millor ones p molt prrximes a la q
	    if (imi.le.ima.and.iqbeg(n)-ima.lt.30/samp) then
	       iew=iqbeg(n)-nint(30/samp)
	       ibw=ibw-nint(30/samp)
	       call buscamaxmin (ibw, iew, dbuf, imi, ymin, ima, ymax)
	    end if
C---------------------------------------------------------------------------
C     Aixs is el criteri de l'Audal:
C	 if (ima.le.imi.and.(ymax.lt.dermax/50.or.abs(ymin).lt.dermax/50)
c     &.or.ima.ge.imi.and.(ymax.lt.dermax/25.or.abs(ymin).lt.dermax/25))
c     & then
C-------------------------------------------------------------------------
C     Ara se'n adopta un altre:
C     si la diferencia d'amplituds entre la posicio de R i la del punt de 
C     maxima amplitud a la finestra on busquem P, es superior a un cert valor
C     considerarem que no existeix ona P. Les amplituds ser`n buscades en 
C     la senyal ecgpb. Primer calculem el punt de m`xima amplitud absoluta en
C     la finestra de bzsqueda de P, i despres podrem comparar-la. De fet con-
C     siderarem bones aquelles ones P l'amplitud de les quals sigui superior a
C     un 7% de la de la ona R, o que tinguin una pendent superior al 1.4% la
C     del complexe QRS. 
C
	    base=0
	    itimes=0
	    do ii=1, nint(15/samp)
	       base=base+(ecgpb(iqbeg(n)-ii))
	       itimes=itimes+1
	    end do
	    base=base/itimes
	    ecgpbmax=0.
	    do ii=ibw,iew
	       if (abs(ecgpb(ii)-base).ge.ecgpbmax) then
		  ecgpbmax=abs(ecgpb(ii)-base)
		  posmax=ii
	       end if
	    end do
c	write(6,*) n, ecgpbmax, abs(ecgpb(iqrs(n))), ymax, ymin,dermax
	    if (ecgpbmax.le.(abs(ecgpb(iqrs(n))-base)/30).or.      
     &          ((ymax.lt.dermax/100.and.abs(ymin).lt.dermax/100).and.
     &          (ymax.lt.abs(ymin)/1.5.or.ymax.gt.abs(ymin)*1.5)) )then
	       ippos1(n)=0
	       ippos2(n)=0
	       ipbeg(n)=0
	       ipend(n)=0
	    else
c	 definim si es normal o invertida
	       if (imi.le.ima) then
c	 tenim ona P invertida
		  iaux=imi
		  yaux=ymin
		  imi=ima
		  ymin=ymax
		  ima=iaux
		  ymax=yaux
	       end if	 

c 	 busquem els dos ceros a dreta i esquerra de la posicio del maxim i
c	 minim respectivament, que correspondran als posibles pics de la ona P
	       icero1=ima
	       call detectar_cero (dbuf, icero1, 'd')
	       icero2=imi
	       call detectar_cero (dbuf, icero2, 'i')
c	 si els dos pics disten mes de un cert valor
c	 if (icero2-icero1.gt.20/samp) then
c	 considerem que la ona p te dos pics 
c	   	ippos1(n)=icero1
c	  	ippos2(n)=icero2
c	   else
c	 considerem que el pic esta en el cente dels dos pics
	       ippos1(n)=nint((icero1+icero2)/2.)
	       ippos2(n)=0
c	 end if
c	 busquem l'inici i final de la ona P segons el criteri del umbral de
c	 la derivada
c        Vigo: fem Kpb=1.65 i ara 1.5 i ara 1.35
	       umbral=ymax/1.35
c        umbral=ymax/kpb
	       ima2=ima
 20	       call creuar_umbral (dbuf, umbral, ima2, iumb, 'i')
	       if ((iqbeg(n)-iumb).ge.240.or.
     &         (iumb.le.itend(n-1).and.n.gt.1.and.itend(n-1).ne.0)) then 
		  ibw=ibw+20
		  if (ibw.gt.iew-20/samp) go to 100
		  call buscamaxmin (ibw,iew,dbuf,imi2,ymin2,ima2,ymax2)
		  goto 20
	       end if

	       call busca_soroll(iumb-40,iumb-5,ecg,soroll)
c        D. Vigo: confirmem que l'inici trobat es bo; fem un test comparant
c              l'algada del'ona p trobada amb l'algada de soroll, i comprobant
c              que l'inici de p no estgui massa aprop del pic de p.
	       if (abs(ecgpb(iumb)-ecgpb(ippos1(n))).lt.(1.5*soroll).and.
     &             (ippos1(n)-iumb).lt.40/samp) then
		  iew2=ima2-nint(15./samp)
		  if (iew2.gt.itend(n-1).and.iew2.ge.ibw) then
		     call buscamaxmin (ibw,iew2,dbuf,imi2,ymin2,ima2,ymax2)
		     if (ymax2.ge.ymax/2) then
			go to 20
		     end if
		  end if
	       end if
	       ipbeg(n)=iumb
	       umbral=ymin/kpe
	 
	       call creuar_umbral (dbuf, umbral, imi, iumb, 'd')
c	 si el umbral es rebassa mes enlla de iqbeg definim el final amb un
c	 criteri de la derivada minima
	       if (iumb.ge.iqbeg(n)) 
     &             call buscamaxmin(imi,iqbeg(n),dbuf,iumb,ymin,iau,yau)
	       ipend(n)=iumb
c	 si encara hem anat massa enll`, fem ipend=iqbeg
	       if (ipend(n).ge.iqbeg(n)) then 
		  if (counter.eq.0) then
		     iew=iew-nint(10/samp) 
		     counter=1
c         write(6,1) ipend(n), iqbeg(n)
c 1 	 format('Pend=',i4,' Qbeg=',i4)
		     go to 10
		  else
		     ipend(n)=iqbeg(n)-1
		  end if
	       end if 

C        comprobem que la diferencia d'algades entre pic i final de P sigui
C        significativa respecte el soroll calculat
	       call busca_soroll(ibw,iew,ecg,soroll)
	       if(abs(ecgpb(ippos1(n))-ecgpb(ipend(n))).le.(1.5*soroll))then
		  goto 100
	       end if

c	 fem una verificacio dels resultats obtinguts per validar la ona P
	       if (ipbeg(n).ge.ipend(n).or.ippos1(n).le.ipbeg(n)
     &             .or.ippos1(n).ge.ipend(n).or.ipbeg(n).lt.itend(n-1)
     &             .or.ipend(n)-ipbeg(n).gt.180/samp) then
c    &             .or.ipend(n)-ipbeg(n).gt.150/samp) then
		  go to 100
	       else
		  go to 110
	       end if
	    end if
 100	    ippos1(n)=0
	    ippos2(n)=0
	    ipbeg(n)=0
	    ipend(n)=0
 110	    iew=iew-50./samp
	    ibw=ibw-50./samp
c       eliminem el bucle de buscar la ona P mes enrera si no la trobem a la
c       primera finestra
c       nofi=.false.
	 end do
	 return
	 end


c----------------------------------------------------------------------------
	 subroutine calcula_rmedio (rrmed, n, samp, irpos)

c	 aquesta subrutina calcula el la distancia mitja entre les ones R
c	 consequtives. Aquells intervals que es desviin molt de la mitja son
c	 rebutjats

	 dimension irpos(8000)

	 perc=0.5
	 if (n.le.1) then
	    if (rrmed.eq.0) then
		rrmed=irpos(2)-irpos(1) 
	        if (rrmed.gt.1000/samp.or.rrmed.lt.500/samp)
     &                   rrmed=nint(650/samp)
	      else
		return
	    end if
	  else
	    if (n.lt.7) then	
	       if (irpos(n)-irpos(n-1).gt.rrmed*(1-perc).and.
     &		irpos(n)-irpos(n-1).lt.rrmed*(1+perc)) then
     	       rrmed=rrmed*(n-2)/(n-1)+(irpos(n)-irpos(n-1))/(n-1)
	       end if
	     else
	       if (irpos(n)-irpos(n-1).gt.rrmed*(1-perc).and.
     &		irpos(n)-irpos(n-1).lt.rrmed*(1+perc)) then
                  rrmed=rrmed*4/5+(irpos(n)-irpos(n-1))/5
	       end if
	    end if
	 end if
         	
c        write(6,1) rrmed
c 1 	 format('rrmed=',f8)
	 
	 return
	 end

c---------------------------------------------------------------------------


	 subroutine test_pic (seny, lpic, samp, t)

c	 aquesta subrutina ens fa un test per comprobar si el pic trobat 
c	 amb el criteri de la derivada nula correspon exactament amb el
c	 pic del senyal 'seny'.
c	 't' ens diu si tenim un pic 'p' o una vall 'v'
	
	 dimension seny(100000)
	 character*1 t

	 kpos=nint(10/samp)
c 	 busquem kpos posicions a dreta i esquerra de lpic
	 laux=lpic
	 if (t.eq.'p') then	 
	   do i=lpic-kpos, lpic+kpos
	   if(seny(i).gt.seny(laux)) laux=i
	   end do
	 else
	   do i=lpic-kpos, lpic+kpos
	   if(seny(i).lt.seny(laux)) laux=i
	   end do
	 end if
	 lpic= laux
	 return 
	 end
	 
c----------------------------------------------------------------------------


 	 subroutine cal_base(n,ippos,ecg,basel,samp,ipend,iqbeg)     
  
 	 dimension basel(8000),ipend(8000), iqbeg(8000),ippos(8000)
	 dimension ecg(100000)

c	 busquem la linea de base com la mitjana dels punts del segment PR

	 nqui=nint(15./samp)
	 ntre=nint(30./samp)
         ntre_q = nint(10./samp)
         nqui_q = nint(5./samp)

	 if (iqbeg(n).ne.0) then
	   if (ippos(n).ne.0) then   
	 	sum=0.
		if (iqbeg(n)-ipend(n).gt.33/samp) then
			do k=ipend(n)+nqui, iqbeg(n)-nqui
			   sum=sum+ecg(k)
			end do
			baselin=sum/(iqbeg(n)-nqui-nqui-ipend(n)+1)
		else
			if(iqbeg(n).eq.ipend(n)) then
			  baselin=ecg(iqbeg(n))
			else
			  do k=ipend(n), iqbeg(n)
			   sum =sum+ecg(k)
			  end do
			  baselin= sum/(iqbeg(n)-ipend(n)+1)
			end if
		end if
	   else
		sum=0.
	 	if(iqbeg(n)-nqui_q.gt.1) then
		  k=iqbeg(n)-nqui_q
		  do while(k.ge.iqbeg(n)-ntre_q-nqui_q.and.k.gt.0)
			sum=sum+ecg(k)
			k=k-1
		  end do
		  if(k.eq.0)then
			baselin=sum/(iqbeg(n)-nqui_q)
		    else
			baselin=sum/(ntre_q+1)
		  end if
	        else
		  do k=1,iqbeg(n)
			sum=sum+ecg(k)
		  end do
		  baselin=sum/iqbeg(n)
		end if
	 end if
	 basel(n)=baselin
	 end if
 	 return 
	 end

c        --------------------------------------------------------------
         subroutine busca_soroll(ibw,iew,ecg,soroll)

         dimension ecg(100000)

         if (ibw.lt.0) ibw=0
         if (iebw.lt.0) iew=0

         soroll=0.
         inici=ibw

         ncont=0
         do while (inici.lt.iew)
            ncont=ncont+1
            ifinal=inici+5
            if (ifinal.gt.iew) then 
                ifinal=iew
            end if

            call buscamaxmin(inici,ifinal,ecg,imi2,ymin2,ima2,ymax2)
            soroll=(ymax2-ymin2)+soroll
            inici=ifinal
         end do 

         if (ncont.gt.0) then
            soroll=abs(soroll/ncont)
         end if

        return 
       end