Add subtraction routines - not fully tested or tweaked yet.

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

CMakeLists.txt

@@ -413,6 +413,7 @@ set (wsjt_FSRCS
   lib/fsk4hf/getfc1w.f90
   lib/fsk4hf/getfc2w.f90
   lib/genqra64.f90
+  lib/fsk4hf/genft8refsig.f90
   lib/genwspr.f90
   lib/geodist.f90
   lib/getlags.f90
@@ -495,6 +496,7 @@ set (wsjt_FSRCS
   lib/spec9f.f90
   lib/stdmsg.f90
   lib/subtract65.f90
+  lib/fsk4hf/subtractft8.f90
   lib/sun.f90
   lib/symspec.f90
   lib/symspec2.f90

lib/fsk4hf/ft8b.f90

@@ -129,6 +129,7 @@ subroutine ft8b(dd0,newdat,nfqso,f1,xdt,nharderrors,dmin,nbadcrc,message,xsnr)
   if(nbadcrc.eq.0) then
      call extractmessage174(decoded,message,ncrcflag,recent_calls,nrecent)
      call genft8(message,msgsent,itone)
+     call subtractft8(dd0,itone,f1,xdt2)
      xsig=0.0
      xnoi=0.0
      do i=1,79

lib/fsk4hf/genft8refsig.f90

@@ -0,0 +1,22 @@
+subroutine genft8refsig(itone,cref,f0)
+  complex cref(79*1920)
+  integer itone(79)
+  real*8 twopi,phi,dphi,dt,xnsps
+  data twopi/0.d0/
+  save twopi
+  if( twopi .lt. 0.1 ) twopi=8.d0*atan(1.d0)
+
+  xnsps=1920.d0
+  dt=1.d0/12000.d0
+  phi=0.d0
+  k=1
+  do i=1,79
+    dphi=twopi*(f0*dt+itone(i)/xnsps)
+    do is=1,1920
+      cref(k)=cmplx(cos(phi),sin(phi))
+      phi=mod(phi+dphi,twopi)
+      k=k+1
+    enddo
+  enddo
+  return
+end subroutine genft8refsig

lib/fsk4hf/subtractft8.f90

@@ -0,0 +1,67 @@
+subroutine subtractft8(dd,itone,f0,dt)
+
+! Subtract an ft8 signal
+!
+! Measured signal  : dd(t)    = a(t)cos(2*pi*f0*t+theta(t))
+! Reference signal : cref(t)  = exp( j*(2*pi*f0*t+phi(t)) )
+! Complex amp      : cfilt(t) = LPF[ dd(t)*CONJG(cref(t)) ]
+! Subtract         : dd(t)    = dd(t) - 2*REAL{cref*cfilt}
+
+  use timer_module, only: timer
+
+  parameter (NMAX=15*12000,NFRAME=1920*79)
+  parameter (NFFT=NMAX,NFILT=1000)
+  real*4  dd(NMAX),dds(NMAX), window(-NFILT/2:NFILT/2)
+  complex cref,camp,cfilt,cw
+  integer itone(79)
+  logical first
+  data first/.true./
+  common/heap8/cref(NFRAME),camp(NMAX),cfilt(NMAX),cw(NMAX)
+  save first
+
+  nstart=dt*12000+1
+  call genft8refsig(itone,cref,f0)
+write(*,*) 'dt ',dt,'nstart ',nstart
+  camp=0.
+  do i=1,nframe
+    id=nstart-1+i 
+    if(id.ge.1.and.id.le.NMAX) camp(i)=dd(id)*conjg(cref(i))
+  enddo
+
+  if(first) then
+! Create and normalize the filter
+     pi=4.0*atan(1.0)
+     fac=1.0/float(nfft)
+     sum=0.0
+     do j=-NFILT/2,NFILT/2
+        window(j)=cos(pi*j/NFILT)**2
+        sum=sum+window(j)
+     enddo
+     cw=0.
+     cw(1:NFILT+1)=window/sum
+     cw=cshift(cw,NFILT/2+1)
+     call four2a(cw,nfft,1,-1,1)
+     cw=cw*fac
+     first=.false.
+  endif
+
+  cfilt=0.0
+  cfilt(1:nframe)=camp(1:nframe)
+  call four2a(cfilt,nfft,1,-1,1)
+  cfilt(1:nfft)=cfilt(1:nfft)*cw(1:nfft)
+  call four2a(cfilt,nfft,1,1,1)
+
+dds=dd
+! Subtract the reconstructed signal
+  do i=1,nframe
+     j=nstart+i-1
+     if(j.ge.1 .and. j.le.NMAX) dd(j)=dd(j)-2*REAL(cfilt(i)*cref(i))
+  enddo
+
+!do i=1,NMAX
+!write(34,*) i,dds(i),dd(i)
+!enddo
+
+  return
+end subroutine subtractft8
+

lib/ft8_decode.f90

@@ -24,9 +24,10 @@ contains
 
   subroutine decode(this,callback,iwave,nfqso,newdat,nutc,nfa,    &
        nfb,nagain,ndepth,nsubmode)
-
+!use wavhdr
     use timer_module, only: timer
     include 'fsk4hf/ft8_params.f90'
+!type(hdr) h
 
     class(ft8_decoder), intent(inout) :: this
     procedure(ft8_decode_callback) :: callback
@@ -69,7 +70,11 @@ contains
 !3051   format(4f9.1,3i5,2x,a22)
 !       flush(51)
     enddo
-
+!h=default_header(12000,NMAX)
+!open(10,file='subtract.wav',status='unknown',access='stream')
+!iwave=nint(dd)
+!write(10) h,iwave
+!close(10)
     return
   end subroutine decode