Much better algorithm for detecting and removing birdies in Q65.

lib/qra/q65/q65.f90

@@ -58,10 +58,8 @@ subroutine q65_dec0(iavg,nutc,iwave,ntrperiod,nfqso,ntol,ndepth,lclearave,  &
   parameter (LN=2176*63)           !LN=LL*NN; LL=64*(mode_q65+2), NN=63
   integer*2 iwave(0:12000*ntrperiod-1)   !Raw data
   integer dat4(13)
-  integer ipk1(1)
   character*37 decoded
   logical first,lclearave
-  integer, allocatable :: hist(:)
   real, allocatable :: s1(:,:)           !Symbol spectra, 1/8-symbol steps
   real, allocatable :: s3(:,:)           !Data-symbol energies s3(LL,63)
   real, allocatable :: ccf1(:)           !CCF(freq) at fixed lag (red)
@@ -101,7 +99,6 @@ subroutine q65_dec0(iavg,nutc,iwave,ntrperiod,nfqso,ntol,ndepth,lclearave,  &
   allocate(s3(-64:LL-65,63))
   allocate(ccf1(-ia2:ia2))
   allocate(ccf2(iz))
-  allocate(hist(LL))
   if(LL.ne.LL0 .or. iz.ne.iz0 .or. jz.ne.jz0 .or. lclearave) then
      if(allocated(s1a)) deallocate(s1a)
      allocate(s1a(iz,jz,0:1))
@@ -134,18 +131,6 @@ subroutine q65_dec0(iavg,nutc,iwave,ntrperiod,nfqso,ntol,ndepth,lclearave,  &
   call pctile(s1(i0-64:i0-65+LL,1:jz),LL*jz,45,base)
   s1=s1/base
 
-! Detect and remove strong non-drifting birdies
-  hist=0
-  do j=1,jz
-     ipk1=maxloc(s1(i0-64:i0+LL-65,j))
-     i=ipk1(1)
-     hist(i)=hist(i)+1
-  enddo
-  nbirdie=80
-  do i=1,LL
-     if(hist(i).gt.nbirdie) s1(i0-65+i,1:jz)=1.0
-  enddo
-  
 ! Apply fast AGC to the symbol spectra
   s1max=20.0                                  !Empirical choice
   do j=1,jz                                   !### Maybe wrong way? ###
@@ -559,6 +544,7 @@ subroutine q65_s1_to_s3(s1,iz,jz,ipk,jpk,LL,mode_q65,sync,s3)
         if(j.ge.1 .and. j.le.jz) s3(-64:LL-65,n)=s1(i1:i2,j)
      enddo
   endif
+  call q65_bzap(s3,LL)                   !Zap birdies
   
   return
 end subroutine q65_s1_to_s3
@@ -611,4 +597,27 @@ subroutine q65_sync_curve(ccf1,ia,ib,rms1)
   return
 end subroutine q65_sync_curve
 
+subroutine q65_bzap(s3,LL)
+
+  parameter (NBZAP=15)
+  real s3(-64:LL-65,63)
+  integer ipk1(1)
+  integer, allocatable :: hist(:)
+
+  allocate(hist(-64:LL-65))
+  hist=0
+  do j=1,63
+     ipk1=maxloc(s3(:,j))
+     i=ipk1(1) - 65
+     hist(i)=hist(i)+1
+  enddo
+  if(maxval(hist).gt.NBZAP) then
+     do i=-64,LL-65
+        if(hist(i).gt.NBZAP) s3(i,1:63)=1.0
+     enddo
+  endif
+
+  return
+end subroutine q65_bzap
+
 end module q65

lib/qra/q65/q65_loops.f90

@@ -56,6 +56,7 @@ subroutine q65_loops(c00,npts2,nsps2,nsubmode,ndepth,jpk0,    &
         call pctile(s3,LL*NN,40,base)
         s3=s3/base
         where(s3(1:LL*NN)>s3lim) s3(1:LL*NN)=s3lim
+        call q65_bzap(s3,LL)                   !Zap birdies
         do ibw=ibwa,ibwb
            ndist=ndf**2 + ndt**2 + (ibw-ibw0)**2
            if(ndist.gt.maxdist) cycle