More work on FT2. New frame format is 16sync + 128codeword. Data chunk is assumed to be 2.5s long. A rudimentary sync routine is implemented which finds the single strongest signal.

CMakeLists.txt

@@ -546,6 +546,7 @@ set (wsjt_FSRCS
   lib/sync4.f90
   lib/sync64.f90
   lib/sync65.f90
+  lib/fsk4hf/getcandidates2.f90
   lib/ft8/sync8.f90
   lib/ft8/sync8d.f90
   lib/sync9.f90

lib/fsk4hf/ft2_params.f90

@@ -5,8 +5,8 @@ parameter (NS=16)                     !Sync symbols (2x8)
 parameter (NN=NS+ND)                  !Total channel symbols (144)
 parameter (NSPS=160)                  !Samples per symbol at 12000 S/s
 parameter (NZ=NSPS*NN)                !Samples in full 1.92 s waveform (23040)
-parameter (NMAX=3*12000)              !Samples in iwave (36,000)
+parameter (NMAX=2.5*12000)              !Samples in iwave (36,000)
-parameter (NFFT1=2*NSPS, NH1=NFFT1/2) !Length of FFTs for symbol spectra
+parameter (NFFT1=400, NH1=NFFT1/2)    !Length of FFTs for symbol spectra
 parameter (NSTEP=NSPS/4)              !Rough time-sync step size
 parameter (NHSYM=NMAX/NSTEP-3)        !Number of symbol spectra (1/4-sym steps)
-parameter (NDOWN=10)                  !Downsample factor
+parameter (NDOWN=16)                  !Downsample factor

lib/fsk4hf/ft2d.f90

@@ -7,23 +7,28 @@ program ft2d
   character*37 decodes(100)
   character*120 data_dir
   character*90 dmsg
-  complex c2(0:3*1200-1)                  !Complex waveform
+  complex c2(0:NMAX/16-1)                  !Complex waveform
+  complex cb(0:NMAX/16-1)
   complex cd(0:144*10-1)                  !Complex waveform
   complex c1(0:9),c0(0:9)
   complex ccor(0:1,144)
-  complex csum,cterm,cc0,cc1
+  complex csum,cterm,cc0,cc1,csync1,csync2
   real*8 fMHz
 
+  real a(5)
   real rxdata(128),llr(128)               !Soft symbols
   real llr2(128)
   real sbits(144),sbits1(144),sbits3(144)
   real ps(0:8191),psbest(0:8191)
   real candidates(100,2)
+  real savg(NH1),sbase(NH1)
   integer ihdr(11)
   integer*2 iwave(NMAX)                 !Generated full-length waveform  
   integer*1 message77(77),apmask(128),cw(128)
   integer*1 hbits(144),hbits1(144),hbits3(144)
+  integer*1 s16(16)
   logical unpk77_success
+  data s16/0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0/
 
   fs=12000.0/NDOWN                       !Sample rate
   dt=1/fs                                !Sample interval after downsample (s)
@@ -81,7 +86,8 @@ program ft2d
   fr2=0.0
   nav=0
   ngood=0
-
+fsum=0.0
+fsum2=0.0
   do ifile=iarg,nargs
      call getarg(ifile,infile)
      j2=index(infile,'.wav')
@@ -91,16 +97,55 @@ program ft2d
      datetime=infile(j2-11:j2-1)
      close(10)
 
+     call getcandidates2(iwave,100.0,3000.0,0.2,2200.0,100,savg,candidates,ncand,sbase)
+
      ndecodes=0
      ncand=1
      do icand=1,ncand
-        fc0=1500.0
+        f0=candidates(icand,1)
         xsnr=1.0
-        istart=6000+8
+        call ft2_downsample(iwave,f0,c2) ! downsample from 160s/Symbol to 10s/Symbol
-        call ft2_downsample(iwave,c2) ! downsample from 160s/Symbol to 10s/Symbol
+ibest=-1
-
+sybest=-99.
-        ib=istart/16
+dfbest=-1.
-        cd=c2(ib:ib+144*10-1) 
+do if=-15,+15
+        df=if
+        a=0.
+        a(1)=-df
+call twkfreq1(c2,NMAX/16,fs,a,cb)
+! 750 samples/second here
+        do is=0,374
+           csync1=0.
+           cterm=1
+           do ib=1,16
+              i1=(ib-1)*10+is
+              i2=i1+136*10
+              if(s16(ib).eq.1) then
+                csync1=csync1+sum(cb(i1:i1+9)*conjg(c1(0:9)))*cterm
+                cterm=cterm*cc1
+              else
+                csync1=csync1+sum(cb(i1:i1+9)*conjg(c0(0:9)))*cterm
+                cterm=cterm*cc0
+              endif
+           enddo
+           if(abs(csync1).gt.sybest) then
+             ibest=is
+             sybest=abs(csync1)
+             dfbest=df
+           endif
+        enddo 
+enddo
+freq=f0+dfbest
+fsum=fsum+freq
+fsum2=fsum+freq*freq
+a=0.
+a(1)=-dfbest
+!write(*,*) 'dfbest ',dfbest
+!dfbest=0.0
+!ibest=187
+call twkfreq1(c2,NMAX/16,fs,a,cb)
+        ib=ibest
+        cd=cb(ib:ib+144*10-1) 
         s2=sum(cd*conjg(cd))/(10*144)
         cd=cd/sqrt(s2)
         do nseq=1,4
@@ -153,8 +198,7 @@ program ft2d
               sbits=sbits3
               hbits=hbits3
            endif
-           rxdata(1:48)=sbits(9:56)
+           rxdata=sbits(17:144)
-           rxdata(49:128)=sbits(65:144)
            rxav=sum(rxdata(1:128))/128.0
            rx2av=sum(rxdata(1:128)*rxdata(1:128))/128.0
            rxsig=sqrt(rx2av-rxav*rxav)
@@ -175,6 +219,7 @@ program ft2d
            if( nharderror.ge.0 ) then
               write(c77,'(77i1)') message77(1:77)
               call unpack77(c77,message,unpk77_success)
+              idupe=0
               do i=1,ndecodes
                  if(decodes(i).eq.message) idupe=1 
               enddo
@@ -182,10 +227,9 @@ program ft2d
               ndecodes=ndecodes+1 
               decodes(ndecodes)=message
               nsnr=nint(xsnr)
-              freq=fMHz + 1.d-6*(fc1+fbest)
 1210          format(a11,2i4,f6.2,f12.7,2x,a22,i3)
-              write(*,1212) datetime(8:11),nsnr,xdt,freq,message,'*',idf,nseq,ijitter,nharderror,nhardmin
+              write(*,1212) datetime(8:11),nsnr,ibest/750.0,freq,message,'*',idf,nseq,ijitter,nharderror,nhardmin
-1212          format(a4,i4,f5.1,f11.6,2x,a22,a1,i5,i5,i5,i5,i5)
+1212          format(a4,i4,f5.1,f11.1,2x,a22,a1,i5,i5,i5,i5,i5)
               goto 888
            endif
         enddo ! nseq
@@ -195,7 +239,9 @@ program ft2d
 
   write(*,1120)
 1120 format("<DecodeFinished>")
-
+favg=fsum/1000.0
+fstd=sqrt(fsum2/1000.0-favg*favg)
+write(*,*) "Mean, std frequency: ",favg,fstd
 999 end program ft2d
 
 subroutine getbitmetric(ib,ps,ns,xmet)
@@ -234,66 +280,7 @@ subroutine downsample2(ci,f0,co)
   return
 end subroutine downsample2
 
-subroutine getcandidate2(c,npts,fs,fa,fb,ncand,candidates)
+subroutine ft2_downsample(iwave,f0,c)
-  parameter(NDAT=200,NFFT1=120*12000/32,NH1=NFFT1/2,NFFT2=120*12000/320,NH2=NFFT2/2)
-  complex c(0:npts-1)                   !Complex waveform
-  complex cc(0:NFFT1-1)
-  complex csfil(0:NFFT2-1)
-  complex cwork(0:NFFT2-1)
-  real bigspec(0:NFFT2-1)
-  complex c2(0:NFFT1-1)                 !Short spectra
-  real s(-NH1+1:NH1)                    !Coarse spectrum
-  real ss(-NH1+1:NH1)                   !Smoothed coarse spectrum
-  real candidates(100,2)
-  integer indx(NFFT2-1)
-  logical first
-  data first/.true./
-  save first,w,df,csfil
-
-  if(first) then
-    df=10*fs/NFFT1
-    csfil=cmplx(0.0,0.0)
-    do i=0,NFFT2-1
-       csfil(i)=exp(-((i-NH2)/20.0)**2)
-    enddo
-    csfil=cshift(csfil,NH2)
-    call four2a(csfil,NFFT2,1,-1,1)
-    first=.false.
-  endif
-
-  cc=cmplx(0.0,0.0)
-  cc(0:npts-1)=c;
-  call four2a(cc,NFFT1,1,-1,1)
-  cc=abs(cc)**2
-  call four2a(cc,NFFT1,1,-1,1)
-  cwork(0:NH2)=cc(0:NH2)*conjg(csfil(0:NH2))
-  cwork(NH2+1:NFFT2-1)=cc(NFFT1-NH2+1:NFFT1-1)*conjg(csfil(NH2+1:NFFT2-1))
-
-  call four2a(cwork,NFFT2,1,+1,1)
-  bigspec=cshift(real(cwork),-NH2)
-  il=NH2+fa/df
-  ih=NH2+fb/df 
-  nnl=ih-il+1
-  call indexx(bigspec(il:il+nnl-1),nnl,indx)
-  xn=bigspec(il-1+indx(nint(0.3*nnl)))
-  bigspec=bigspec/xn
-  ncand=0
-  do i=il,ih
-    if((bigspec(i).gt.bigspec(i-1)).and. &
-       (bigspec(i).gt.bigspec(i+1)).and. &
-       (bigspec(i).gt.1.15).and.ncand.lt.100) then 
-         ncand=ncand+1
-         candidates(ncand,1)=df*(i-NH2)
-         candidates(ncand,2)=10*log10(bigspec(i))-30.0
-    endif
-  enddo
-!  do i=1,ncand
-!    write(*,*) i,candidates(i,1),candidates(i,2)
-!  enddo 
-  return
-end subroutine getcandidate2
-
-subroutine ft2_downsample(iwave,c)
 
 ! Input: i*2 data in iwave() at sample rate 12000 Hz
 ! Output: Complex data in c(), sampled at 1200 Hz
@@ -310,7 +297,7 @@ subroutine ft2_downsample(iwave,c)
   df=12000.0/NMAX
   x=iwave
   call four2a(x,NMAX,1,-1,0)             !r2c FFT to freq domain
-  i0=nint(1500.0/df)
+  i0=nint(f0/df)
   c1(0)=cx(i0)
   do i=1,NFFT2/2
      c1(i)=cx(i0+i)

lib/fsk4hf/ft2sim.f90

@@ -89,7 +89,7 @@ program ft2sim
   call sgran()
 
   do ifile=1,nfiles
-     k=nint((xdt+0.5)/dt)
+     k=nint((xdt+0.25)/dt)
      ia=k
      phi=0.0 
      c0=0.0

lib/fsk4hf/genft2.f90

@@ -9,13 +9,8 @@ subroutine genft2(msg0,ichk,msgsent,i4tone,itype)
 !   - msgsent  message as it will be decoded
 !   - i4tone   array of audio tone values, 0 or 1
 !   - itype    message type 
-!                 1 = standard message  "Call_1 Call_2 Grid/Rpt"
+!                 1 = 77 bit message 
-!                 2 = type 1 prefix
+!                 7 = 16 bit message     "<Call_1 Call2> Rpt"
-!                 3 = type 1 suffix
-!                 4 = type 2 prefix
-!                 5 = type 2 suffix
-!                 6 = free text (up to 13 characters)
-!                 7 = short message     "<Call_1 Call2> Rpt"
 
   use iso_c_binding, only: c_loc,c_size_t
   use packjt77
@@ -27,24 +22,15 @@ subroutine genft2(msg0,ichk,msgsent,i4tone,itype)
   integer*1 codeword(128)
   integer*1 msgbits(77) 
   integer*1 bitseq(144)                   !Tone #s, data and sync (values 0-1)
-  integer*1 s8(8)
+  integer*1 s16(16)
-  real*8 pp(12)
   real*8 xi(864),xq(864),pi,twopi
-  data s8/0,1,1,1,0,0,1,0/
+  data s16/0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0/
   equivalence (ihash,i1hash)
-  logical first,unpk77_success
+  logical unpk77_success
-  data first/.true./
-  save
 
-  if(first) then
+  nsym=128
-    first=.false.
+  pi=4.0*atan(1.0)
-    nsym=128
+  twopi=8.*atan(1.0)
-    pi=4.0*atan(1.0)
-    twopi=8.*atan(1.0)
-    do i=1,12
-      pp(i)=sin((i-1)*pi/12)
-    enddo
-  endif
 
   message(1:37)=' ' 
   itype=1
@@ -89,36 +75,12 @@ subroutine genft2(msg0,ichk,msgsent,i4tone,itype)
      call encode_128_90(msgbits,codeword)
 
 !Create 144-bit channel vector:
-!8-bit sync word + 48 bits + 8-bit sync word + 80 bits
      bitseq=0 
-     bitseq(1:8)=s8
+     bitseq(1:16)=s16
-     bitseq(9:56)=codeword(1:48)
+     bitseq(17:144)=codeword
-     bitseq(57:64)=s8
-     bitseq(65:144)=codeword(49:128)
 
      i4tone=bitseq
-
-!     bitseq=2*bitseq-1
-!     xq(1:6)=bitseq(1)*pp(7:12)   !first bit is mapped to 1st half-symbol on q
-!     do i=1,71
-!       is=(i-1)*12+7
-!       xq(is:is+11)=bitseq(2*i+1)*pp
-!     enddo 
-!     xq(864-5:864)=bitseq(1)*pp(1:6)   !last half symbol
-!     do i=1,72                                    
-!       is=(i-1)*12+1
-!       xi(is:is+11)=bitseq(2*i)*pp
-!     enddo
-! Map I and Q  to tones. 
-!    i4tone=0 
-!    do i=1,72
-!      i4tone(2*i-1)=(bitseq(2*i)*bitseq(2*i-1)+1)/2;
-!      i4tone(2*i)=-(bitseq(2*i)*bitseq(mod(2*i,144)+1)-1)/2;
-!    enddo
   endif
 
-! Flip polarity
-!  i4tone=-i4tone+1
-
 999 return
 end subroutine genft2

lib/fsk4hf/getcandidates2.f90

@@ -0,0 +1,51 @@
+subroutine getcandidates2(id,nfa,nfb,syncmin,nfqso,maxcand,savg,candidate,   &
+     ncand,sbase)
+
+! For now, hardwired to find the largest peak in the average spectrum
+
+  include 'ft2_params.f90'
+  real s(NH1,NHSYM)
+  real savg(NH1),savsm(NH1)
+  real sbase(NH1)
+  real x(NFFT1)
+  complex cx(0:NH1)
+  real candidate(3,maxcand)
+  integer*2 id(NMAX)
+  integer*1 s8(8)
+  data s8/0,1,1,1,0,0,1,0/
+  equivalence (x,cx)
+
+! Compute symbol spectra, stepping by NSTEP steps.  
+  savg=0.
+  tstep=NSTEP/12000.0                         
+  df=12000.0/NFFT1                            !3.125 Hz
+  fac=1.0/300.0
+  do j=1,NHSYM
+     ia=(j-1)*NSTEP + 1
+     ib=ia+NSPS-1
+     x(1:NSPS)=fac*id(ia:ib)
+     x(NSPS+1:)=0.
+     call four2a(x,NFFT1,1,-1,0)              !r2c FFT
+     do i=1,NH1
+        s(i,j)=real(cx(i))**2 + aimag(cx(i))**2
+     enddo
+     savg=savg + s(1:NH1,j)                   !Average spectrum
+  enddo
+  savsm=0.
+  do i=2,NH1-1
+    savsm(i)=sum(savg(i-1:i+1))/3.
+  enddo
+
+  imax=-1
+  xmax=-99.
+  do i=2,NH1-1
+    if(savsm(i).gt.savsm(i-1).and.savsm(i).gt.savsm(i+1).and.savsm(i).gt.xmax) then
+      xmax=savsm(i) 
+      imax=i
+    endif
+  enddo
+  ncand=1
+  candidate(1,1)=imax*df 
+
+return
+end subroutine getcandidates2