FT4: Some groundwork for subtraction.

CMakeLists.txt

@@ -555,6 +555,7 @@ set (wsjt_FSRCS
   lib/stdmsg.f90
   lib/subtract65.f90
   lib/ft8/subtractft8.f90
+  lib/ft4/subtractft4.f90
   lib/sun.f90
   lib/symspec.f90
   lib/symspec2.f90

lib/ft4/ft4sim_mult.f90

@@ -10,6 +10,7 @@ program ft4sim_mult
   type(hdr) h                            !Header for .wav file
   character arg*12,fname*17,cjunk*4
   character msg37*37,msgsent37*37,c77*77
+  complex cwave0((NN+2)*NSPS)
   real wave0((NN+2)*NSPS)
   real wave(NZZ)
   real tmp(NZZ)
@@ -63,7 +64,8 @@ program ft4sim_mult
         call pack77(msg37,i3,n3,c77)
         call genft4(msg37,0,msgsent37,itone)
         nwave0=(NN+2)*NSPS
-        call gen_ft4wave(itone,NN,NSPS,12000.0,f0,wave0,nwave0)
+        icmplx=0
+        call gen_ft4wave(itone,NN,NSPS,12000.0,f0,cwave0,wave0,icmplx,nwave0)
 
         k0=nint((xdt+0.5)/dt)
         if(k0.lt.1) k0=1

lib/ft4/gen_ft4wave.f90

@@ -1,6 +1,7 @@
-subroutine gen_ft4wave(itone,nsym,nsps,fsample,f0,wave,nwave)
+subroutine gen_ft4wave(itone,nsym,nsps,fsample,f0,cwave,wave,icmplx,nwave)
 
   real wave(nwave)
+  complex cwave(nwave)
   real pulse(6144)              !512*4*3
   real dphi(0:240000-1)
   integer itone(nsym)
@@ -34,19 +35,32 @@ subroutine gen_ft4wave(itone,nsym,nsps,fsample,f0,wave,nwave)
   phi=0.0
   dphi = dphi + twopi*f0*dt                          !Shift frequency up by f0
   wave=0.
+  cwave=0.
   k=0
   do j=0,nwave-1
      k=k+1
+     if(icmplx.eq.0) then
         wave(k)=sin(phi)
+     else
+        cwave(k)=cmplx(cos(phi),sin(phi))
+     endif
      phi=mod(phi+dphi(j),twopi)
   enddo
 
 ! Compute the ramp-up and ramp-down symbols
+  if(icmplx.eq.0) then
      wave(1:nsps)=wave(1:nsps) *                                          &
           (1.0-cos(twopi*(/(i,i=0,nsps-1)/)/(2.0*nsps)))/2.0
      k1=(nsym+1)*nsps+1
      wave(k1:k1+nsps-1)=wave(k1:k1+nsps-1) *                              &
           (1.0+cos(twopi*(/(i,i=0,nsps-1)/)/(2.0*nsps)))/2.0
+  else
+     cwave(1:nsps)=cwave(1:nsps) *                                          &
+          (1.0-cos(twopi*(/(i,i=0,nsps-1)/)/(2.0*nsps)))/2.0
+     k1=(nsym+1)*nsps+1
+     cwave(k1:k1+nsps-1)=cwave(k1:k1+nsps-1) *                              &
+          (1.0+cos(twopi*(/(i,i=0,nsps-1)/)/(2.0*nsps)))/2.0
+  endif
 
   return
 end subroutine gen_ft4wave

lib/ft4/subtractft4.f90

@@ -0,0 +1,66 @@
+subroutine subtractft4(dd,itone,f0,dt)
+
+! Subtract an ft4 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=6*12000,NFRAME=(103+2)*79)
+  parameter (NFFT=NMAX,NFILT=1400)
+  real*4  dd(NMAX), window(-NFILT/2:NFILT/2), xjunk
+  complex cref,camp,cfilt,cw
+  integer itone(79)
+  logical first
+  data first/.true./
+  common/heap8/cref(NFRAME),camp(NMAX),cfilt(NMAX),cw(NMAX),xjunk(NFRAME)
+  save first
+
+  nstart=dt*12000+1
+  nsym=79
+  nsps=1920
+  fs=12000.0
+  icmplx=1
+  bt=1.0      ! Temporary compromise?
+  call gen_ft4wave(itone,nsym,nsps,fs,f0,cref,xjunk,icmplx,NFRAME)
+  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)
+
+! 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
+
+  return
+end subroutine subtractft4
+

lib/ft8/gen_ft8wave.f90

@@ -42,6 +42,7 @@ subroutine gen_ft8wave(itone,nsym,nsps,bt,fsample,f0,cwave,wave,icmplx,nwave)
   phi=0.0
   dphi = dphi + twopi*f0*dt                          !Shift frequency up by f0
   wave=0.
+  cwave=0.
   k=0
   do j=nsps,nsps+nwave-1                             !Don't include dummy symbols
      k=k+1

widgets/mainwindow.cpp

@@ -106,7 +106,7 @@ extern "C" {
                     float xjunk[], float wave[], int* icmplx, int* nwave);
 
   void gen_ft4wave_(int itone[], int* nsym, int* nsps, float* fsample, float* f0,
-                    float wave[], int* nwave);
+                    float xjunk[], float wave[], int* icmplx, int* nwave);
 
   void gen4_(char* msg, int* ichk, char* msgsent, int itone[],
                int* itext, fortran_charlen_t, fortran_charlen_t);
@@ -3745,7 +3745,9 @@ void MainWindow::guiUpdate()
             float fsample=48000.0;
             float f0=ui->TxFreqSpinBox->value() - m_XIT;
             int nwave=(nsym+2)*nsps;
-            gen_ft4wave_(const_cast<int *>(itone),&nsym,&nsps,&fsample,&f0,foxcom_.wave,&nwave);
+            int icmplx=0;
+            gen_ft4wave_(const_cast<int *>(itone),&nsym,&nsps,&fsample,&f0,foxcom_.wave,
+                           foxcom_.wave,&icmplx,&nwave);
           }
 
           if(SpecOp::EU_VHF==m_config.special_op_id()) {