Updates and additions to wsprcpm files.

CMakeLists.txt

@@ -596,6 +596,13 @@ set (wsjt_FSRCS
   lib/wqencode.f90
   lib/wspr_downsample.f90
   lib/zplot9.f90
+  lib/fsk4hf/wsprcpmsim.f90
+  lib/fsk4hf/wsprcpm_wav.f90
+  lib/fsk4hf/wsprcpmd.f90
+  lib/fsk4hf/bpdecode204.f90
+  lib/fsk4hf/osd204.f90
+  lib/fsk4hf/genwsprcpm.f90
+  lib/fsk4hf/encode204.f90
   )
 
 # temporary workaround for a gfortran v7.3 ICE on Fedora 27 64-bit
@@ -1327,6 +1334,12 @@ target_link_libraries (ft4sim_mult wsjt_fort wsjt_cxx)
 add_executable (record_time_signal Audio/tools/record_time_signal.cpp)
 target_link_libraries (record_time_signal wsjt_cxx wsjt_qtmm wsjt_qt)
 
+add_executable (wsprcpmsim lib/fsk4hf/wsprcpmsim.f90 wsjtx.rc)
+target_link_libraries (wsprcpmsim wsjt_fort wsjt_cxx)
+
+add_executable (wsprcpmd lib/fsk4hf/wsprcpmd.f90 wsjtx.rc)
+target_link_libraries (wsprcpmd wsjt_fort wsjt_cxx)
+
 endif(WSJT_BUILD_UTILS)
 
 # build the main application

lib/fsk4hf/genwsprcpm.f90

@@ -40,10 +40,10 @@ subroutine genwsprcpm(msg,msgsent,itone)
    id7=idat(7)
    if(id7.lt.0) id7=id7+256
    id7=id7/64
-write(*,*) 'idat ',idat
+   write(*,*) 'idat ',idat
    icrc=crc14(c_loc(idat),9)
-write(*,*) 'icrc: ',icrc
+   write(*,*) 'icrc: ',icrc
-write(*,'(a6,b16.16)') 'icrc: ',icrc
+   write(*,'(a6,b16.16)') 'icrc: ',icrc
    call wqdecode(idat,msgsent,itype)
    print*,msgsent,itype
    write(cbits,1004) idat(1:6),id7,iand(icrc,z'3FFF')
@@ -52,7 +52,7 @@ write(*,'(a6,b16.16)') 'icrc: ',icrc
    read(cbits,1006) msgbits(1:64)
 1006 format(64i1)
 
-write(*,'(50i1,1x,14i1,1x,4i1)') msgbits
+   write(*,'(50i1,1x,14i1,1x,4i1)') msgbits
 
    call encode204(msgbits,codeword)      !Encode the test message
 

lib/fsk4hf/wsprcpmd.f90

@@ -25,6 +25,7 @@ program wsprcpmd
    complex cd(0:216*10-1)                   !Complex waveform
    complex c1(0:9,0:1),c0(0:9,0:1)
    complex ccor(0:1,216)
+   complex cp(0:1,0:1)
    complex csum,cterm
    real*8 fMHz
    real rxdata(ND),llr(204)               !Soft symbols
@@ -57,8 +58,34 @@ program wsprcpmd
    tt=NSPS*dt                             !Duration of "itone" symbols (s)
    baud=1.0/tt                            !Keying rate for "itone" symbols (baud)
    txt=NZ*dt                              !Transmission length (s)
-  h=1.00                                !h=0.8 seems to be optimum for AWGN sensitivity (not for fading)
+   h=0.50                                !h=0.8 seems to be optimum for AWGN sensitivity (not for fading)
    twopi=8.0*atan(1.0)
+   pi=4.0*atan(1.0)
+
+   nargs=iargc()
+   if(nargs.lt.1) then
+      print*,'Usage:   wsprcpmd [-a <data_dir>] [-f fMHz] [-c ncoh] [-h h] file1 [file2 ...]'
+      go to 999
+   endif
+   iarg=1
+   data_dir="."
+   call getarg(iarg,arg)
+   if(arg(1:2).eq.'-a') then
+      call getarg(iarg+1,data_dir)
+      iarg=iarg+2
+   endif
+   call getarg(iarg,arg)
+   if(arg(1:2).eq.'-f') then
+      call getarg(iarg+1,arg)
+      read(arg,*) fMHz
+      iarg=iarg+2
+   endif
+   ncoh=1
+   if(arg(1:2).eq.'-h') then
+      call getarg(iarg+1,arg)
+      read(arg,*) h
+      iarg=iarg+2
+   endif
 
    isync2(1:100)=isync(1:100)
    isync2(101:104)=0  ! This is *not* backwards.
@@ -71,6 +98,7 @@ program wsprcpmd
 !   0    1   1
 !   1    0   2
 !   1    1   3
+
    dphi=twopi*baud*(h/2.0)*dt*20  ! dt*10 is samp interval after downsample
    do j=0,1
       if(j.eq.0) then
@@ -88,36 +116,10 @@ program wsprcpmd
          phi1=mod(phi1+dphi1,twopi)
          phi0=mod(phi0+dphi0,twopi)
       enddo
+      cp(1,j)=cmplx(cos(phi1),sin(phi1))
+      cp(0,j)=cmplx(cos(phi0),sin(phi0))
    enddo
 
-  nargs=iargc()
-  if(nargs.lt.1) then
-     print*,'Usage:   wsprcpmd [-a <data_dir>] [-f fMHz] [-c ncoh] file1 [file2 ...]'
-     go to 999
-  endif
-  iarg=1
-  data_dir="."
-  call getarg(iarg,arg)
-  if(arg(1:2).eq.'-a') then
-     call getarg(iarg+1,data_dir)
-     iarg=iarg+2
-  endif
-  call getarg(iarg,arg)
-  if(arg(1:2).eq.'-f') then
-     call getarg(iarg+1,arg)
-     read(arg,*) fMHz
-     iarg=iarg+2
-  endif
-  ncoh=1
-  npdi=16
-  if(arg(1:2).eq.'-c') then
-     call getarg(iarg+1,arg)
-     read(arg,*) ncoh
-     iarg=iarg+2
-     npdi=16/ncoh
-  endif
-!  write(*,*) 'ncoh: ',ncoh,' npdi: ',npdi
-  
    open(13,file=trim(data_dir)//'/ALL_WSPR.TXT',status='unknown',   &
       position='append')
 
@@ -154,35 +156,42 @@ program wsprcpmd
       npts=120*12000.0/32.0
       nsync=16
       call getcandidate2(c2,npts,fs,fa,fb,ncand,candidates)         !First approx for freq
-ndecodes=0
+      ndecodes=0
-do icand=1,ncand
+      do icand=1,ncand
          fc0=candidates(icand,1)
          xsnr=candidates(icand,2)
          xmax=-1e32
          do i=-5,5
             ft=fc0+i*0.2
-      call noncoherent_frame_sync(c2,ft,isync2,is,xf1)
+            call noncoherent_frame_sync(c2,h,ft,isync2,is,xf1)
             if(xf1.gt.xmax) then
                xmax=xf1
                fc1=ft
                is0=is
             endif
          enddo
-    call coherent_preamble_fsync(c2,ipreamble,nsync,NSPS,is0,fc1,fcest,xp1)
+         fcest=fc1
-    call noncoherent_frame_sync(c2,fcest,isync2,istart,xf2)
+         imode=0 ! refine freq
-    write(*,'(i5,i5,i5,4(f11.5,2x))') ifile,iii,istart,fc0,fc1,fcest
+         call coherent_preamble_fsync(c2,h,ipreamble,nsync,NSPS,is0,fcest,imode,xp0)
+         imode=1 ! refine istart
+         istart=is0
+         call coherent_preamble_fsync(c2,h,ipreamble,nsync,NSPS,istart,fcest,imode,xp1)
+         write(*,'(i5,i5,i5,6(f11.5,2x))') ifile,is0,istart,fc0,fc1,fcest,xf1,xp0,xp1
 
-do ijitter=0,4
+!genie sync
+!istart=375
+!fcest=0.0
+         do ijitter=0,4
             io=-10*(ijitter/2+1)
             if(mod(ijitter,2).eq.0) io=10*(ijitter/2)
             ib=max(0,istart+io)
             cframe=c2(ib:ib+216*200-1)
-    call downsample2(cframe,fcest,cd)
+            call downsample2(cframe,fcest,h,cd)
 
             s2=sum(cd*conjg(cd))/(10*216)
             cd=cd/sqrt(s2)
 
-    do nseq=1,7
+            do nseq=1,5
                if( nseq.eq.1 ) then  ! noncoherent single-symbol detection
                   sbits1=0.0
                   do ibit=1,216
@@ -215,9 +224,10 @@ do ijitter=0,4
                         cterm=1.0
                         k=1
                         do i=nbit-1,0,-1
+                           j=isync2(ibit-(nbit/2+1)+k)
                            ibb=iand(iseq/(2**i),1)
                            csum=csum+ccor(ibb,ibit-(nbit/2+1)+k)*cterm
-              cterm=-cterm
+                           cterm=cterm*conjg(cp(ibb,j))
                            k=k+1
                         enddo
                         ps(iseq)=abs(csum)
@@ -237,6 +247,8 @@ do ijitter=0,4
                   sbits=sbits3
                   hbits=hbits3
                endif
+               if( count(hbits(101:116).ne.ipreamble) .gt.7 ) cycle
+
                rxdata(1:100)=sbits(1:100)
                rxdata(101:200)=sbits(117:216);
                rxav=sum(rxdata(1:200))/200.0
@@ -285,14 +297,14 @@ do ijitter=0,4
                   nfdot=0
                   write(13,1210) datetime,0,nsnr,xdt,freq,message,nfdot
 1210              format(a11,2i4,f6.2,f12.7,2x,a22,i3)
-        write(*,1212) datetime(8:11),nsnr,xdt,freq,nfdot,message,'*',idf,nseq,ijitter,nharderror,nhardmin
+                  write(*,1212) datetime(8:11),nsnr,xdt,freq,nfdot,message,'*',nseq,ijitter,nharderror,nhardmin
-1212    format(a4,i4,f5.1,f11.6,i3,2x,a22,a1,i5,i5,i5,i5,i5)
+1212              format(a4,i4,f5.1,f11.6,i3,2x,a22,a1,i5,i5,i5,i5)
                   goto 888
                endif
             enddo ! nseq
-enddo !jitter
+         enddo !jitter
 888      continue
-enddo !candidate list
+      enddo !candidate list
    enddo !files
 
    write(*,1120)
@@ -301,7 +313,9 @@ enddo !candidate list
 999 end program wsprcpmd
 
 
-subroutine coherent_preamble_fsync(c2,ipreamble,nsync,nsps,istart,fc0,fc1,xmax)
+subroutine coherent_preamble_fsync(c2,h,ipreamble,nsync,nsps,istart,fc,imode,xmax)
+! imode=0: refine fc using given istart 
+! imode=1: refine istart using given fc
    complex c2(0:120*12000/32-1)
    complex cpreamble(0:16*200-1)
    complex ctmp1(0:4*16*200-1)
@@ -309,12 +323,11 @@ subroutine coherent_preamble_fsync(c2,ipreamble,nsync,nsps,istart,fc0,fc1,xmax)
    complex ccohp(0:15)
    integer ipreamble(nsync)
    logical first/.true./
-    save dt,first,h,twopi,cpreamble
+   save dt,first,twopi,cpreamble
 
    if(first) then
       baud=12000.0/6400.0
       dt=32.0/12000.0
-       h=1.00    
       twopi=8.0*atan(1.0)
       k=0
       phi=0.0
@@ -331,7 +344,7 @@ subroutine coherent_preamble_fsync(c2,ipreamble,nsync,nsps,istart,fc0,fc1,xmax)
       first=.false.
    endif
 
-    dphi=twopi*fc0*dt
+   dphi=twopi*fc*dt
    ctwkp=cmplx(0.0,0.0)
    phi=0
    do i=0,nsync*nsps-1
@@ -340,6 +353,20 @@ subroutine coherent_preamble_fsync(c2,ipreamble,nsync,nsps,istart,fc0,fc1,xmax)
    enddo
    ipstart=istart+100*200
    ctmp1=0.0
+   xmax=0.0
+   if(imode.eq.1) then !refine DT with given fc
+      do iii=-40,40,5
+         ctmp1(0:16*200-1)=c2(ipstart+iii:ipstart+iii+16*200-1)*conjg(ctwkp)
+         xx=abs(sum(ctmp1))
+         if(xx.gt.xmax) then
+            iiibest=iii
+            xmax=xx
+         endif
+      enddo
+      istart=istart+iiibest
+      return
+   endif
+! else refine fc with given DT
    ctmp1(0:16*200-1)=c2(ipstart:ipstart+16*200-1)*conjg(ctwkp)
    call four2a(ctmp1,4*16*200,1,-1,1)             !c2c FFT to freq domain
    xmax=0.0
@@ -358,24 +385,22 @@ subroutine coherent_preamble_fsync(c2,ipreamble,nsync,nsps,istart,fc0,fc1,xmax)
    xp1=abs(ctmp1(ishift+1))
    xint=(log(xm1)-log(xp1))/(log(xm1)+log(xp1)-2*log(x0))
    delta2=delta+xint*dfp/2.0
-    fc1=fc0+delta2
+   fc=fc+delta2
    return
 end subroutine coherent_preamble_fsync
 
-
+subroutine noncoherent_frame_sync(c2,h,fc,isync2,istart,ssmax)
-subroutine noncoherent_frame_sync(c2,fc,isync2,istart,ssmax)
    complex c2(0:120*12000/32-1)
    complex ct0(0:199),ct1(0:199),ct2(0:199),ct3(0:199)
    integer isync2(216)
 
    twopi=8.0*atan(1.0)
-   h=1.0
    dt=32.0/12000.0
    baud=12000.0/6400.0
    imax=370  ! defines dt search range (375 samples/s)
    ssmax=-1e32
    izero=375
-   do it = -imax,imax,10
+   do it = -imax,imax,5
 ! noncoherent wspr-type dt estimation
       dp0=twopi*(fc-1.5*h*baud)*dt
       dp1=twopi*(fc-0.5*h*baud)*dt
@@ -397,6 +422,8 @@ subroutine noncoherent_frame_sync(c2,fc,isync2,istart,ssmax)
       enddo
       ss=0.0
       totp=0.0
+      avp=0.0
+      xc=0.0
       do is=1,216
          i0=izero+it+(is-1)*200
          p0=abs(sum(c2(i0:i0+199)*conjg(ct0)))
@@ -409,15 +436,21 @@ subroutine noncoherent_frame_sync(c2,fc,isync2,istart,ssmax)
          p3=sqrt(p3)
 
          totp=totp+p0+p1+p2+p3
+         avp=avp+(p0+p1+p2+p3)/4.0
 !        cmet=(p1+p3)-(p0+p2)
          cmet=max(p1,p3)-max(p0,p2)  ! This works better near threshold SNR
          if(isync2(is).eq.0) ss=ss-cmet
          if(isync2(is).eq.1) ss=ss+cmet
+         if(isync2(is).eq.0) xc=xc+max(p0,p2)
+         if(isync2(is).eq.1) xc=xc+max(p1,p3)
       enddo
       ss=ss/totp
-      if(ss.gt.ssmax) then
+      sy=xc/avp
+!      if(ss.gt.ssmax) then
+      if(sy.gt.ssmax) then
          ioffset=it
-        ssmax=ss
+!        ssmax=ss
+         ssmax=sy
       endif
    enddo
    istart=izero+ioffset
@@ -436,7 +469,7 @@ subroutine getmetric2(ib,ps,ns,xmet)
    return
 end subroutine getmetric2
 
-subroutine downsample2(ci,f0,co)
+subroutine downsample2(ci,f0,h,co)
    parameter(NI=216*200,NH=NI/2,NO=NI/20)  ! downsample from 200 samples per symbol to 10
    complex ci(0:NI-1),ct(0:NI-1)
    complex co(0:NO-1)
@@ -448,7 +481,7 @@ subroutine downsample2(ci,f0,co)
    ct=cshift(ct,i0)
    co=0.0
    co(0)=ct(0)
-  b=8.0
+   b=max(1.0,h)*8.0
    do i=1,NO/2
       arg=(i*df/b)**2
       filt=exp(-arg)
@@ -480,7 +513,7 @@ subroutine getcandidate2(c,npts,fs,fa,fb,ncand,candidates)
       df=10*fs/NFFT1
       csfil=cmplx(0.0,0.0)
       do i=0,NFFT2-1
-       csfil(i)=exp(-((i-NH2)/20.0)**2)
+         csfil(i)=exp(-((i-NH2)/20.0)**2)  ! revisit this
       enddo
       csfil=cshift(csfil,NH2)
       call four2a(csfil,NFFT2,1,-1,1)
@@ -505,17 +538,15 @@ subroutine getcandidate2(c,npts,fs,fa,fb,ncand,candidates)
    bigspec=bigspec/xn
    ncand=0
    do i=il,ih
+      write(21,*) i*df,bigspec(i)
       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
+         candidates(ncand,2)=10*log10(bigspec(i)-1)-30.0
       endif
    enddo
-!  do i=1,ncand
-!    write(*,*) i,candidates(i,1),candidates(i,2)
-!  enddo 
    return
 end subroutine getcandidate2
 

lib/fsk4hf/wsprcpmsim.f90

@@ -17,9 +17,9 @@ program wsprcpmsim
 
 ! Get command-line argument(s)
    nargs=iargc()
-  if(nargs.ne.8) then
+   if(nargs.ne.9) then
-     print*,'Usage:   wsprmsksim "message"       f0  DT fsp del  nwav nfiles snr'
+      print*,'Usage:   wsprcpmsim "message"       f0  DT fsp del  nwav nfiles snr  h'
-     print*,'Example: wsprmsksim "K1ABC FN42 30" 50 0.0 0.1 1.0  1      10   -33'
+      print*,'Example: wsprcpmsim "K1ABC FN42 30" 50 0.0 0.1 1.0  1      10   -33  1.0'
       go to 999
    endif
    call getarg(1,msg)                     !Message to be transmitted
@@ -37,6 +37,8 @@ program wsprcpmsim
    read(arg,*) nfiles                     !Number of files
    call getarg(8,arg)
    read(arg,*) snrdb                      !SNR_2500
+   call getarg(9,arg)
+   read(arg,*) h                          !h
 
    twopi=8.0*atan(1.0)
    fs=12000.0/NDOWN                       !
@@ -54,7 +56,6 @@ program wsprcpmsim
 1000 format('f0:',f9.3,'   DT:',f6.2,'   txt:',f6.1,'   SNR:',f6.1,    &
       '   fspread:',f6.1,'   delay:',f6.1,'  nfiles:',i3,2x,a22)
 
-  h=1.00
    c0=0.
    k=-1 + nint(xdt/dt)
    do j=1,NN
@@ -70,11 +71,10 @@ program wsprcpmsim
       c=c0
       if(nwav.eq.0) then
          if( fspread .ne. 0.0 .or. delay .ne. 0.0 ) then
-           call watterson(c,NMAX/NDOWN,fs,delay,fspread)
+            ntot=NMAX/NDOWN
+            nsig=NN*NSPS
+            call watterson(c,ntot,nsig,fs,delay,fspread)
          endif
-!do i=0,NMAX/NDOWN-1
-!write(23,*) i,real(c(i)),imag(c(i))
-!enddo
          c=c*sig
          if(snrdb.lt.90) then
             do i=0,NMAX/NDOWN-1                   !Add gaussian noise at specified SNR
@@ -90,9 +90,6 @@ program wsprcpmsim
          nmin=2
          write(10) fname,nmin,fMHz,c      !Save to *.c2 file
          close(10)
-!do i=0,NMAX/NDOWN-1
-!write(57,*) i,real(c(i)),imag(c(i))
-!enddo
       else
          call wsprcpm_wav(baud,xdt,h,f0,itone,snrdb,iwave)
          hwav=default_header(12000,NMAX)