Add QRA64 birdie-zapping capability. Needs further testing!

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/map65@7509 ab8295b8-cf94-4d9e-aec4-7959e3be5d79

libm65/CMakeLists.txt

@@ -114,6 +114,7 @@ set (FSRCS
   pfxdump.f90
   qra64b.f90
   qra64c.f90
+  qra64zap.f90
   recvpkt.f90
   rfile3a.f90
   s3avg.f90

libm65/qra64b.f90

@@ -3,12 +3,13 @@ subroutine qra64b(nutc,nqd,ikhz,mousedf,ntol,xpol,mycall_12,hiscall_12,   &
 
   parameter (NFFT1=5376000)              !56*96000
   parameter (NFFT2=336000)               !56*6000 (downsampled by 1/16)
-  complex cx(0:NFFT2-1),cy(0:NFFT2-1)
   complex ca(NFFT1),cb(NFFT1)            !FFTs of raw x,y data
+  complex cx(0:NFFT2-1),cy(0:NFFT2-1)
   logical xpol
   character*12 mycall_12,hiscall_12
   character*6 hisgrid_6
   common/cacb/ca,cb
+  data nzap/3/
 
   open(17,file='red.dat',status='unknown')
   df=96000.0/NFFT1
@@ -22,6 +23,8 @@ subroutine qra64b(nutc,nqd,ikhz,mousedf,ntol,xpol,mycall_12,hiscall_12,   &
   cy(nh+1:NFFT2-1)=cb(k0-nh+1:k0-1)
   cy=fac*cy
 
+  if(nzap.gt.0) call qra64zap(cx,cy,nzap)
+
 ! Transform back to time domain with sample rate 6000 Hz.
   call four2a(cx,NFFT2,1,-1,1)
   call four2a(cy,NFFT2,1,-1,1)

libm65/qra64zap.f90

@@ -0,0 +1,60 @@
+subroutine qra64zap(cx,cy,nzap)
+
+  parameter (NFFT1=5376000)              !56*96000
+  parameter (NFFT2=336000)               !56*6000 (downsampled by 1/16)
+  complex cx(0:NFFT2-1),cy(0:NFFT2-1)
+  real s(-1312:1312)
+  integer iloc(1)
+
+  slimit=3.0
+  sbottom=1.5
+  nadd=128
+  nblks=NFFT2/nadd
+  nbh=nblks/2
+  k=-1
+  s=0.
+  df=nadd*96000.0/NFFT1
+  do i=1,nblks
+     j=i
+     if(j.gt.nblks/2) j=j-nblks
+     do n=1,nadd
+        k=k+1
+        s(j)=s(j) + real(cx(k))**2 + aimag(cx(k))**2 +                 &
+             real(cy(k))**2 + aimag(cy(k))**2
+     enddo
+  enddo
+  call pctile(s,nblks,45,base)
+  s=s/base
+  do nzap=1,3
+     iloc=maxloc(s)
+     ipk=iloc(1)-1313
+     smax=s(ipk)
+     nw=3
+     do n=1,3
+        nw=2*nw
+        s1=maxval(s(ipk-2*nw:ipk-nw))
+        s2=maxval(s(ipk+nw:ipk+2*nw))
+        if(smax.gt.slimit .and. s1.lt.sbottom .and. s2.lt.sbottom) then
+           s(ipk-nw:ipk+nw)=1.0
+           i0=ipk
+           if(i0.lt.0) i0=i0+2625
+           ia=(i0-nw)*nadd
+           ib=(i0+nw)*nadd
+           cx(ia:ib)=0.
+           cy(ia:ib)=0.
+!           print*,'!',nzap,ipk*df,nw
+           exit
+        endif
+     enddo
+  enddo
+     
+!  rewind 75
+!  do i=-nbh,nbh   
+!     freq=i*df
+!     write(75,3001) freq,s(i)
+!3001 format(2f12.3)
+!  enddo
+!  flush(75)
+
+  return
+end subroutine qra64zap

libm65/sync64.f90

@@ -135,9 +135,9 @@ subroutine sync64(c0,nf1,nf2,nfqso,ntol,mode64,emedelay,dtx,f0,jpk,sync,  &
   rms2=sqrt(sq/40.0)
   sync2=10.0*log10(a(2)/rms2)
 
-  slimit=2.0
+  slimit=3.0
   rewind 17
-  rewind 76
+!  rewind 76
   do i=2,iz-2*nskip-1,3
      x=i
      z=(x-a(3))/(0.5*a(4))
@@ -150,11 +150,11 @@ subroutine sync64(c0,nf1,nf2,nfqso,ntol,mode64,emedelay,dtx,f0,jpk,sync,  &
      freq=j*df3-3000.0
      ss=(s0a(j-1)+s0a(j)+s0a(j+1))/3.0
      if(ss.gt.slimit) write(17,1110) freq,ss
-     write(76,1110) freq,ss,yfit
 1110 format(3f10.3)
+!     write(76,1110) freq,ss,yfit
   enddo
   flush(17)
-  flush(76)
+!  flush(76)
 
   return
 end subroutine sync64