code:
subroutine tweak1(ca,jz,f0,cb)

! Shift frequency of analytic signal ca, with output to cb

  complex ca(jz),cb(jz)
  real*8 twopi
  complex*16 w,wstep
  complex w4
  data twopi/0.d0/
  save twopi

  if(twopi.eq.0.d0) twopi=8.d0*atan(1.d0)
  w=1.d0
  dphi=twopi*f0/12000.d0
  wstep=cmplx(cos(dphi),sin(dphi))
  do i=1,jz
     w=w*wstep
     w4=w
     cb(i)=w4*ca(i)
  enddo

  return
end subroutine tweak1


code:
module wavhdr
  type hdr
     character*4 ariff
     integer*4 lenfile
     character*4 awave
     character*4 afmt
     integer*4 lenfmt
     integer*2 nfmt2
     integer*2 nchan2
     integer*4 nsamrate
     integer*4 nbytesec
     integer*2 nbytesam2
     integer*2 nbitsam2
     character*4 adata
     integer*4 ndata
  end type hdr

  contains

    function default_header(nsamrate,npts)
      type(hdr) default_header,h
      h%ariff='RIFF'
      h%awave='WAVE'
      h%afmt='fmt '
      h%lenfmt=16
      h%nfmt2=1
      h%nchan2=1
      h%nsamrate=nsamrate
      h%nbitsam2=16
      h%nbytesam2=h%nbitsam2 * h%nchan2 / 8
      h%adata='data'
      h%nbytesec=h%nsamrate * h%nbitsam2 * h%nchan2 / 8
      h%ndata=2*npts
      h%lenfile=h%ndata + 44 - 8
      default_header=h
    end function default_header

    subroutine set_wsjtx_wav_params(fMHz,mode,nsubmode,ntrperiod,id2)

      parameter (NBANDS=23,NMODES=11)
      character*8 mode,modes(NMODES)
      integer*2 id2(4)
      integer iperiod(7)
      real fband(NBANDS)
      data fband/0.137,0.474,1.8,3.5,5.1,7.0,10.14,14.0,18.1,21.0,24.9,  &
           28.0,50.0,144.0,222.0,432.0,902.0,1296.0,2304.0,3400.0,       &
           5760.0,10368.0,24048.0/
      data modes/'Echo','FSK441','ISCAT','JT4','JT65','JT6M','JT9',      &
           'JT9+JT65','JTMS','JTMSK','WSPR'/
      data iperiod/5,10,15,30,60,120,900/

      dmin=1.e30
      iband=0
      do i=1,NBANDS
         if(abs(fMHz-fband(i)).lt.dmin) then
            dmin=abs(fMHz-fband(i))
            iband=i
         endif
      enddo

      imode=0
      do i=1,NMODES
         if(mode.eq.modes(i)) imode=i
      enddo

      ip=0
      do i=1,7
         if(ntrperiod.eq.iperiod(i)) ip=i
      enddo

      id2(1)=iband
      id2(2)=imode
      id2(3)=nsubmode
      id2(4)=ip

      return
    end subroutine set_wsjtx_wav_params

    subroutine get_wsjtx_wav_params(id2,band,mode,nsubmode,ntrperiod,ok)

      parameter (NBANDS=23,NMODES=11)
      character*8 mode,modes(NMODES)
      character*6 band,bands(NBANDS)
      integer*2 id2(4)
      integer iperiod(7)
      logical ok
      data modes/'Echo','FSK441','ISCAT','JT4','JT65','JT6M','JT9',    &
           'JT9+JT65','JTMS','JTMSK','WSPR'/
      data iperiod/5,10,15,30,60,120,900/
      data bands/'2190m','630m','160m','80m','60m','40m','30m','20m',  &
           '17m','15m','12m','10m','6m','2m','1.25m','70cm','33cm',    &
           '23cm','13cm','9cm','6cm','3cm','1.25cm'/

      ok=.true.
      if(id2(1).lt.1 .or. id2(1).gt.NBANDS) ok=.false.
      if(id2(2).lt.1 .or. id2(2).gt.NMODES) ok=.false.
      if(id2(3).lt.1 .or. id2(3).gt.8) ok=.false.
      if(id2(4).lt.1 .or. id2(4).gt.7) ok=.false.

      if(ok) then
         band=bands(id2(1))
         mode=modes(id2(2))
         nsubmode=id2(3)
         ntrperiod=iperiod(id2(4))
      endif

      return
    end subroutine get_wsjtx_wav_params

end module wavhdr

license: mit


code:
subroutine wspr5_downsample(iwave,c)

! Input: i*2 data in iwave() at sample rate 12000 Hz
! Output: Complex data in c(), sampled at 400 Hz

  include 'wsprlf_params.f90'
  parameter (NMAX=300*12000,NFFT2=NMAX/30)
  integer*2 iwave(NMAX)
  complex c(0:NZ-1)
  complex c1(0:NFFT2-1)
  complex cx(0:NMAX/2)
  real x(NMAX)
  equivalence (x,cx)

  df=12000.0/NMAX
  x=iwave
  call four2a(x,NMAX,1,-1,0)             !r2c FFT to freq domain
  i0=nint(1500.0/df)
  c1(0)=cx(i0)
  do i=1,NFFT2/2
     c1(i)=cx(i0+i)
     c1(NFFT2-i)=cx(i0-i)
  enddo
  c1=c1/NFFT2
  call four2a(c1,NFFT2,1,1,1)            !c2c FFT back to time domain
  c=c1(0:NZ-1)

  return
end subroutine wspr5_downsample


code:
subroutine wspr5_wav(baud,xdt,f0,itone,snrdb,iwave)

! Generate iwave() from itone().

  include 'wsprlf_params.f90'
  parameter (NMAX=300*12000)
  integer itone(NN)
  integer*2 iwave(NMAX)
  real*8 twopi,dt,dphi0,dphi1,dphi,phi
  real dat(NMAX)

  twopi=8.d0*atan(1.d0)
  dt=1.d0/12000.d0

  dat=0.
  if(snrdb.lt.90) then
     do i=1,NMAX
        dat(i)=gran()          !Generate gaussian noise
     enddo
     bandwidth_ratio=2500.0/6000.0
     sig=sqrt(2*bandwidth_ratio)*10.0**(0.05*snrdb)
  else
     sig=1.0
  endif

  dphi0=twopi*(f0-0.25d0*baud)*dt
  dphi1=twopi*(f0+0.25d0*baud)*dt
  phi=0.d0
  k=nint(xdt/dt)
  do j=1,NN
     dphi=dphi0
     if(itone(j).eq.1) dphi=dphi1
     if(k.eq.0) phi=-dphi
     do i=1,NSPS0
        k=k+1
        phi=phi+dphi
        if(phi.gt.twopi) phi=phi-twopi
        xphi=phi
        if(k.gt.0 .and. k.le.NMAX) dat(k)=dat(k) + sig*sin(xphi)
     enddo
  enddo
  fac=32767.0
  rms=100.0
  if(snrdb.ge.90.0) iwave=nint(fac*dat)
  if(snrdb.lt.90.0) iwave=nint(rms*dat)

  return
end subroutine wspr5_wav


code:
program wspr5d

! Decode WSPR-LF data read from *.c5 or *.wav files.

! WSPR-LF is a potential WSPR-like mode intended for use at LF and MF.
! It uses an LDPC (300,60) code, OQPSK modulation, and 5 minute T/R sequences.

! Reception and Demodulation algorithm:
!   1. Compute coarse spectrum; find fc1 = approx carrier freq
!   2. Mix from fc1 to 0; LPF at +/- 0.75*R
!   3. Square, FFT; find peaks near -R/2 and +R/2 to get fc2
!   4. Mix from fc2 to 0
!   5. Fit cb13 (central part of csync) to c -> lag, phase
!   6. Fit complex ploynomial for channel equalization
!   7. Get soft bits from equalized data

! Still to do: find and decode more than one signal in the specified passband.

  include 'wsprlf_params.f90'
  parameter (NMAX=300*12000)
  character arg*8,message*22,cbits*50,infile*80,fname*16,datetime*11
  character*120 data_dir
  complex csync(0:NZ-1)                 !Sync symbols only, from cbb
  complex c(0:NZ-1)                     !Complex waveform
  complex c1(0:NZ-1)                    !Complex waveform
  complex zz(NS+ND)                     !Complex symbol values (intermediate)
  complex z
  real*8 fMHz
  real rxdata(ND),llr(ND)               !Soft symbols
  real pp(2*NSPS)                       !Shaped pulse for OQPSK
  real a(5)                             !For twkfreq1
  real aa(20),bb(20)                    !Fitted polyco's
  real fpks(20)
  integer id(NS+ND)                     !NRZ values (+/-1) for Sync and Data
  integer ierror(NS+ND)
  integer isync(48)                     !Long sync vector
  integer ib13(13)                      !Barker 13 code
  integer ihdr(11)
  integer*8 n8
  integer*2 iwave(NMAX)                 !Generated full-length waveform
  integer*1 idat(7)
  integer*1 decoded(KK),apmask(ND),cw(ND)
  data ib13/1,1,1,1,1,-1,-1,1,1,-1,1,-1,1/

  nargs=iargc()
  if(nargs.lt.2) then
     print*,'Usage:   wspr5d [-a <data_dir>] [-f fMHz] 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

  open(13,file=trim(data_dir)//'/ALL_WSPR.TXT',status='unknown',   &
       position='append')
!  maxn=4                                 !Default value
  maxn=2
  twopi=8.0*atan(1.0)
  fs=NSPS*12000.0/NSPS0                  !Sample rate
  dt=1.0/fs                              !Sample interval (s)
  tt=NSPS*dt                             !Duration of "itone" symbols (s)
  ts=2*NSPS*dt                           !Duration of OQPSK symbols (s)
  baud=1.0/tt                            !Keying rate for "itone" symbols (baud)
  txt=NZ*dt                              !Transmission length (s)

  do i=1,N2                              !Half-sine pulse shape
     pp(i)=sin(0.5*(i-1)*twopi/(2*NSPS))
  enddo
  n8=z'cbf089223a51'
  do i=1,48
     isync(i)=-1
     if(iand(n8,1).eq.1) isync(i)=1
     n8=n8/2
  enddo

! Define array id() for sync symbols
  id=0
  do j=1,48                             !First group of 48
     id(2*j-1)=2*isync(j)
  enddo
  do j=1,13                             !Barker 13 code
     id(j+96)=2*ib13(j)
  enddo
  do j=1,48                             !Second group of 48
     id(2*j+109)=2*isync(j)
  enddo

  csync=0.
  do j=1,205
     if(abs(id(j)).eq.2) then
        ia=nint((j-0.5)*N2)
        ib=ia+N2-1
        csync(ia:ib)=pp*id(j)/abs(id(j))
     endif
  enddo

write(*,*) 'iarg, nargs ',iarg,nargs
  do ifile=iarg,nargs
     call getarg(ifile,infile)
     open(10,file=infile,status='old',access='stream')
     j1=index(infile,'.c5')
     j2=index(infile,'.wav')
     if(j1.gt.0) then
        read(10,end=999) fname,ntrmin,fMHz,c
        read(fname(8:11),*) nutc
        write(datetime,'(i11)') nutc
     else if(j2.gt.0) then
        read(10,end=999) ihdr,iwave
        read(infile(j2-4:j2-1),*) nutc
        datetime=infile(j2-11:j2-1)
        call wspr5_downsample(iwave,c)
     else
        print*,'Wrong file format?'
        go to 999
     endif
     close(10)
     fa=100.0
     fb=150.0
     call getfc1w(c,fs,fa,fb,fc1,xsnr)         !First approx for freq
     npeaks=20
     call getfc2w(c,csync,npeaks,fs,fc1,fpks)      !Refined freq

     a(1)=-fc1
     a(2:5)=0.
     call twkfreq1(c,NZ,fs,a,c)       !Mix c down by fc1+fc2
! Find time offset xdt
     amax=0.
     jpk=0
     iaa=0
     ibb=NZ-1
     jmax=1260
     do j=-jmax,jmax,NSPS/8
        ia=j
        ib=NZ-1+j
        if(ia.lt.0) then
           ia=0
           iaa=-j
        else
           iaa=0
        endif
        if(ib.gt.NZ-1) then
           ib=NZ-1
           ibb=NZ-1-j
        endif
        z=sum(c(ia:ib)*conjg(csync(iaa:ibb)))
        if(abs(z).gt.amax) then
           amax=abs(z)
           jpk=j
        endif
     enddo
     xdt=jpk/fs
xdt=1.0
jpk=fs*xdt
     do i=0,NZ-1
        j=i+jpk
        if(j.ge.0 .and. j.lt.NZ) c1(i)=c(j)
     enddo

     nterms=maxn
     do itry=1,npeaks
        nhard0=0
        nhardsync0=0
        ifer=1
        a(1)=-fpks(itry)
        a(2:5)=0.
        call twkfreq1(c1,NZ,fs,a,c)       !Mix c1 into c
        call cpolyfitw(c,pp,id,maxn,aa,bb,zz,nhs)
        call msksoftsymw(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhardsync0)
        if(nhardsync0.gt.35) cycle
        rxav=sum(rxdata)/ND
        rx2av=sum(rxdata*rxdata)/ND
        rxsig=sqrt(rx2av-rxav*rxav)
        rxdata=rxdata/rxsig
        ss=0.84
        llr=2.0*rxdata/(ss*ss)
        apmask=0
        max_iterations=40
        ifer=0
        call bpdecode300(llr,apmask,max_iterations,decoded,niterations,cw)
        nhardmin=0
        if(niterations.lt.0) call osd300(llr,apmask,5,decoded,cw,nhardmin,dmin)
        nbadcrc=0
        call chkcrc10(decoded,nbadcrc)
        if(nbadcrc.ne.0) ifer=1
        if(ifer.eq.0) exit
     enddo                                !Freq dither loop
     message='                      '
     if(ifer.eq.0) then
        write(cbits,1100) decoded(1:50)
1100    format(50i1)
        read(cbits,1102) idat
1102    format(6b8,b2)
        idat(7)=ishft(idat(7),6)
        call wqdecode(idat,message,itype)
        nsnr=nint(xsnr)
!        freq=fMHz + 1.d-6*(fc1+fc2)
!        freq=fMHz + 1.d-6*(fc1+fpks(itry))
        freq=fc1+fpks(itry)
        nfdot=0
        write(13,1110) datetime,0,nsnr,xdt,freq,message,nfdot
1110    format(a11,2i4,f6.2,f12.7,2x,a22,i3)
        write(*,1112) datetime(8:11),nsnr,xdt,freq,nfdot,message,itry,nhardmin
!1112    format(a4,i4,f5.1,f11.6,i3,2x,a22,i4)
1112    format(a4,i4,f8.3,f8.3,i3,2x,a22,i4,i4)
     endif
  enddo                                   ! ifile loop
  write(*,1120)
1120 format("<DecodeFinished>")

999 end program wspr5d


code:
program wspr5d

! Decode WSPR-LF data read from *.c5 or *.wav files.

! WSPR-LF is a potential WSPR-like mode intended for use at LF and MF.
! It uses an LDPC (300,60) code, OQPSK modulation, and 5 minute T/R sequences.
!
! Still to do: find and decode more than one signal in the specified passband.

!  include 'wsprlf_params.f90'

  parameter (NDOWN=30)
  parameter (KK=60)
  parameter (ND=300)
  parameter (NS=109)
  parameter (NR=3)
  parameter (NN=NR+NS+ND)
  parameter (NSPS0=8640)
  parameter (NSPS=16)
  parameter (N2=2*NSPS)
  parameter (NZ=NSPS*NN)
  parameter (NZ400=288*NN)
  parameter (NMAX=300*12000)

  character arg*8,message*22,cbits*50,infile*80,fname*16,datetime*11
  character*120 data_dir
  complex csync(0:NZ-1)                 !Sync symbols only, from cbb
  complex c400(0:NZ400-1)                     !Complex waveform
  complex c(0:NZ-1)                     !Complex waveform
  complex cd(0:NZ-1)                    !Complex waveform
  complex ca(0:NZ-1)                    !Complex waveform
  complex zz,zzsum
  real*8 fMHz
  real rxdata(ND),llr(ND)               !Soft symbols
  real pp(32)                       !Shaped pulse for OQPSK
  real sbits(412),softbits(9)
  real fpks(20)
  integer id(NS+ND)                     !NRZ values (+/-1) for Sync and Data
  integer isync(48)                     !Long sync vector
  integer ib13(13)                      !Barker 13 code
  integer ihdr(11)
  integer*8 n8
  integer*2 iwave(NMAX)                 !Generated full-length waveform
  integer*1 idat(7)
  integer*1 decoded(KK),apmask(ND),cw(ND)
  integer*1 hbits(412),bits(13)
  logical reset
  data ib13/1,1,1,1,1,-1,-1,1,1,-1,1,-1,1/

  nargs=iargc()
  if(nargs.lt.2) then
     print*,'Usage:   wspr5d [-a <data_dir>] [-f fMHz] 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

  open(13,file=trim(data_dir)//'/ALL_WSPR.TXT',status='unknown',   &
       position='append')
  maxn=8                                 !Default value
  twopi=8.0*atan(1.0)
  fs=NSPS*12000.0/NSPS0                  !Sample rate
  dt=1.0/fs                              !Sample interval (s)
  tt=NSPS*dt                             !Duration of "itone" symbols (s)
  ts=2*NSPS*dt                           !Duration of OQPSK symbols (s)
  baud=1.0/tt                            !Keying rate for "itone" symbols (baud)
  txt=NZ*dt                              !Transmission length (s)

  do i=1,32                              !Half-sine pulse shape
     pp(i)=sin(0.5*(i-1)*twopi/(32))
  enddo
  n8=z'cbf089223a51'
  do i=1,48
     isync(i)=-1
     if(iand(n8,1).eq.1) isync(i)=1
     n8=n8/2
  enddo

! Define array id() for sync symbols
  id=0
  do j=1,48                             !First group of 48
     id(2*j-1)=2*isync(j)
  enddo
  do j=1,13                             !Barker 13 code
     id(j+96)=2*ib13(j)
  enddo
  do j=1,48                             !Second group of 48
     id(2*j+109)=2*isync(j)
  enddo

  csync=0.
  do j=1,205
     if(abs(id(j)).eq.2) then
        ia=nint((j-0.5)*N2)
        ib=ia+N2-1
        csync(ia:ib)=pp*id(j)/abs(id(j))
     endif
  enddo

  do ifile=iarg,nargs
    call getarg(ifile,infile)
    open(10,file=infile,status='old',access='stream')
    j1=index(infile,'.c5')
    j2=index(infile,'.wav')
    if(j1.gt.0) then
       read(10,end=999) fname,ntrmin,fMHz,c400
       read(fname(8:11),*) nutc
       write(datetime,'(i11)') nutc
    else if(j2.gt.0) then
       read(10,end=999) ihdr,iwave
       read(infile(j2-4:j2-1),*) nutc
       datetime=infile(j2-11:j2-1)
       call wspr5_downsample(iwave,c400)
    else
       print*,'Wrong file format?'
       go to 999
    endif
    close(10)

    fa=100.0
    fb=150.0
    fs400=400.0
    call getfc1(c400,fs400,fa,fb,fc1,xsnr)         !First approx for freq
    npeaks=5
    call getfc2(c400,npeaks,fs400,fc1,fpks)      !Refined freq

!    do idf=1,npeaks ! consider the top npeak peaks
    do idf=1,1  ! for genie-aided sync
      fc1=125.0 ! genie provided
      fc2=0.0   ! from the genie
!      fc2=fpks(idf)
      call downsample(c400,fc1+fc2,cd)
      s2=sum(cd*conjg(cd))/(16*412)
      cd=cd/sqrt(s2)
      do is=0,0 ! dt search range is zeroed for genie-aided sync
        idt=is/2
        if( mod(is,2).eq. 1 ) idt=-(is+1)/2
        xdt=real(22+idt)/22.222 - 1.0
        ca=cshift(cd,22+idt)
        zzsum=0.0
        do iseq=3,4
          if(iseq.eq.4) then
            k=1-2*3
            nseq=9
            istep=3*4
          else
            k=1-2*iseq
            nseq=iseq*3
            istep=iseq*4
          endif
          do i=1,408,istep
            j=(i+1)*16
            if(iseq.eq.4) then
! phase=-1.18596900
! For now, average complex corr. coeffs over the entire frame to
! estimate phase
               phase=atan2(imag(zzsum),real(zzsum))
               k=k+3*2
               call mskcohdet(nseq,ca(j),pp,id(k),softbits,phase)
            else
               k=k+iseq*2
               call mskseqdet(nseq,ca(j),pp,id(k),softbits,1,zz)
               zzsum=zzsum+zz
            endif
            sbits(i+1)=softbits(1)
            sbits(i+2)=softbits(2)
            if( id(k+1) .ne. 0 ) sbits(i+2)=id(k+1)*25
              sbits(i+3)=softbits(3)
            if( iseq .ge. 2 ) then
              sbits(i+5)=softbits(4)
              sbits(i+6)=softbits(5)
            if( id(k+3) .ne. 0 ) sbits(i+6)=id(k+3)*25
              sbits(i+7)=softbits(6)
              if( iseq .ge. 3 ) then
                sbits(i+9)=softbits(7)
                sbits(i+10)=softbits(8)
                if( id(k+5) .ne. 0 ) sbits(i+10)=id(k+5)*25
                sbits(i+11)=softbits(9)
              endif
            endif
          enddo
          j=1
          do i=1,205
            if( abs(id(i)) .ne. 2 ) then
              rxdata(j)=sbits(2*i-1)
              j=j+1
            endif
          enddo
          do i=1,204
            rxdata(j)=sbits(2*i)
            j=j+1
          enddo
          rxav=sum(rxdata)/ND
          rx2av=sum(rxdata*rxdata)/ND
          rxsig=sqrt(rx2av-rxav*rxav)
          rxdata=rxdata/rxsig
          sigma=1.20
          llr=2*rxdata/(sigma*sigma)
          apmask=0
          max_iterations=40

          ifer=0
          call bpdecode300(llr,apmask,max_iterations,decoded,niterations,cw)
! niterations will be equal to the Hamming distance between hard received word and the codeword
          nhardmin=0
          if(niterations.lt.0) call osd300(llr,apmask,5,decoded,cw,nhardmin,dmin)
          if(nhardmin.gt.0) niterations=nhardmin
          nbadcrc=0
          call chkcrc10(decoded,nbadcrc)
          if(nbadcrc.ne.0) ifer=1

          if( ifer.eq.0 ) then
            write(cbits,1200) decoded(1:50)
1200        format(50i1)
            read(cbits,1202) idat
1202        format(6b8,b2)
            idat(7)=ishft(idat(7),6)
            call wqdecode(idat,message,itype)
            nsnr=nint(xsnr)
!            freq=fMHz + 1.d-6*(fc1+fc2)
            freq=fc1+fc2
            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,is,iseq,niterations
!1212        format(a4,i4,f5.1,f11.6,i3,2x,a22,a1,i3,i3,i3,i4)
1212        format(a4,i4,f8.3,f8.3,i3,2x,a22,a1,i3,i3,i3,i3,i4)
            goto 888
          endif
        enddo !iseq
      enddo
    enddo
888 continue
  enddo

  write(*,1120)
1120 format("<DecodeFinished>")

999 end program wspr5d

subroutine getmetric(ib,ps,xmet)
  real ps(0:511)
  xm1=0
  xm0=0
  do i=0,511
    if( iand(i/ib,1) .eq. 1 .and. ps(i) .gt. xm1 ) xm1=ps(i)
    if( iand(i/ib,1) .eq. 0 .and. ps(i) .gt. xm0 ) xm0=ps(i)
  enddo
  xmet=xm1-xm0
  return
end subroutine getmetric

subroutine mskseqdet(ns,cdat,pp,bsync,softbits,ncoh,zz)
!
! Detect sequences of 3, 6, or 9 bits (ns).
! Sync bits are assumed to be known.
!
complex cdat(16*12),cbest(16*12),cideal(16*12)
complex cdf(16*12),cfac,zz
real cm(0:511),cmbest(0:511)
real pp(32),softbits(9)
integer bit(13),bestbits(13),sgn(13)
integer bsync(7)

twopi=8.0*atan(1.0)
dt=30.0*18.0/12000.0
cmax=0;
fbest=0.0;
np=2**ns-1
idfmax=40
if( ncoh .eq. 1 ) idfmax=0
do idf=0,idfmax
  if( mod(idf,2).eq.0 ) deltaf=idf/2*0.02
  if( mod(idf,2).eq.1 ) deltaf=-(idf+1)/2*0.02
  dphi=twopi*deltaf*dt
  cfac=cmplx(cos(dphi),sin(dphi))
  cdf=1.0
  do i=2,16*(ns-1)
    cdf(i)=cdf(i-1)*cfac
  enddo

  cm=0
  ibflag=0
  do i=0,np
    bit(1)=(bsync(1)+2)/4
    bit(2)=iand(i/(2**(ns-1)),1)
    bit(3)=iand(i/(2**(ns-2)),1)
    if( bsync(2).ne.0 ) then ! force the barker bits
      bit(3)=(bsync(2)+2)/4
    endif
    bit(4)=iand(i/(2**(ns-3)),1)
    bit(5)=(bsync(3)+2)/4

    if( ns .ge. 6 ) then
      bit(6)=iand(i/(2**(ns-4)),1)
      bit(7)=iand(i/(2**(ns-5)),1)
      if( bsync(4).ne.0 ) then ! force the barker bits
        bit(7)=(bsync(4)+2)/4
      endif
      bit(8)=iand(i/(2**(ns-6)),1)
      bit(9)=(bsync(5)+2)/4
      if( ns .eq. 9 ) then
        bit(10)=iand(i/4,1)
        bit(11)=iand(i/2,1)
        if( bsync(6).ne.0 ) then ! force the barker bits
          bit(11)=(bsync(6)+2)/4
        endif
        bit(12)=iand(i/1,1)
        bit(13)=(bsync(7)+2)/4
      endif
    endif

    sgn=2*bit-1
    cideal(1:16)   =cmplx(sgn(1)*pp(17:32),sgn(2)*pp(1:16))
    cideal(17:32)  =cmplx(sgn(3)*pp(1:16),sgn(2)*pp(17:32))
    cideal(33:48)  =cmplx(sgn(3)*pp(17:32),sgn(4)*pp(1:16))
    cideal(49:64)  =cmplx(sgn(5)*pp(1:16),sgn(4)*pp(17:32))
    if( ns .ge. 6 ) then
      cideal(65:80)  =cmplx(sgn(5)*pp(17:32),sgn(6)*pp(1:16))
      cideal(81:96)  =cmplx(sgn(7)*pp(1:16),sgn(6)*pp(17:32))
      cideal(97:112) =cmplx(sgn(7)*pp(17:32),sgn(8)*pp(1:16))
      cideal(113:128)=cmplx(sgn(9)*pp(1:16),sgn(8)*pp(17:32))
      if( ns .eq. 9 ) then
        cideal(129:144)  =cmplx(sgn(9)*pp(17:32),sgn(10)*pp(1:16))
        cideal(145:160)  =cmplx(sgn(11)*pp(1:16),sgn(10)*pp(17:32))
        cideal(161:176) =cmplx(sgn(11)*pp(17:32),sgn(12)*pp(1:16))
        cideal(177:192)=cmplx(sgn(13)*pp(1:16),sgn(12)*pp(17:32))
      endif
    endif
    cideal=cideal*cdf
    cm(i)=abs(sum(cdat(1:64*ns/3)*conjg(cideal(1:64*ns/3))))/1.e3
    if( cm(i) .gt. cmax ) then
      ibflag=1
      cmax=cm(i)
      bestbits=bit
      cbest=cideal
      fbest=deltaf
      zz=sum(cdat*conjg(cbest))/1.e3
    endif
  enddo
  if( ibflag .eq. 1 ) then ! new best found
    cmbest=cm
  endif
enddo
softbits=0.0
call getmetric(1,cmbest,softbits(ns))
call getmetric(2,cmbest,softbits(ns-1))
call getmetric(4,cmbest,softbits(ns-2))
if( ns .ge. 6 ) then
  call getmetric(8,cmbest,softbits(ns-3))
  call getmetric(16,cmbest,softbits(ns-4))
  call getmetric(32,cmbest,softbits(ns-5))
  if( ns .eq. 9 ) then
    call getmetric(64,cmbest,softbits(3))
    call getmetric(128,cmbest,softbits(2))
    call getmetric(256,cmbest,softbits(1))
  endif
endif
end subroutine mskseqdet

subroutine mskcohdet(ns,cdat,pp,bsync,softbits,phase)
!
! Coherent demodulate blocks of 9 bits (ns).
!
  complex cdat(16*12),crot(16*12)
  real pp(32),softbits(9)

  np=2**ns-1

  softbits=0.0
  crot=cdat*cmplx(cos(phase),-sin(phase))
  softbits(1)=sum(imag(crot(1:32)*pp))
  softbits(2)=sum(real(crot(17:48)*pp))
  softbits(3)=sum(imag(crot(33:64)*pp))
  softbits(4)=sum(imag(crot(65:96)*pp))
  softbits(5)=sum(real(crot(81:112)*pp))
  softbits(6)=sum(imag(crot(97:128)*pp))
  softbits(7)=sum(imag(crot(129:160)*pp))
  softbits(8)=sum(real(crot(145:176)*pp))
  softbits(9)=sum(imag(crot(161:192)*pp))
  softbits=softbits/64.
end subroutine mskcohdet

subroutine downsample(ci,f0,co)
  parameter(NI=412*288,NO=NI/18)
  complex ci(0:NI-1),ct(0:NI-1)
  complex co(0:NO-1)

  df=400.0/NI
  ct=ci
  call four2a(ct,NI,1,-1,1)             !c2c FFT to freq domain
  i0=nint(f0/df)
  co=0.0
  co(0)=ct(i0)
!  b=3.0 !optimized for sequence detection
  b=6.0
  do i=1,NO/2
     arg=(i*df/b)**2
     filt=exp(-arg)
     co(i)=ct(i0+i)*filt
     co(NO-i)=ct(i0-i)*filt
  enddo
  co=co/NO
  call four2a(co,NO,1,1,1)            !c2c FFT back to time domain
  return
end subroutine downsample

subroutine getfc1(c,fs,fa,fb,fc1,xsnr)

!  include 'wsprlf_params.f90'
  parameter (NZ=288*412)
  parameter (NSPS=288)
  parameter (N2=2*NSPS)
  parameter (NFFT1=16*NSPS)
  parameter (NH1=NFFT1/2)

  complex c(0:NZ-1)                     !Complex waveform
  complex c2(0:NFFT1-1)                 !Short spectra
  real s(-NH1+1:NH1)                    !Coarse spectrum
  nspec=NZ/N2
  df1=fs/NFFT1
  s=0.
  do k=1,nspec
     ia=(k-1)*N2
     ib=ia+N2-1
     c2(0:N2-1)=c(ia:ib)
     c2(N2:)=0.
     call four2a(c2,NFFT1,1,-1,1)
     do i=0,NFFT1-1
        j=i
        if(j.gt.NH1) j=j-NFFT1
        s(j)=s(j) + real(c2(i))**2 + aimag(c2(i))**2
     enddo
  enddo
!        call smo121(s,NFFT1)
  smax=0.
  ipk=0
  fc1=0.
  ia=nint(fa/df1)
  ib=nint(fb/df1)
  do i=ia,ib
     f=i*df1
     if(s(i).gt.smax) then
        smax=s(i)
        ipk=i
        fc1=f
     endif
!            write(51,3001) f,s(i),db(s(i))
! 3001       format(f10.3,e12.3,f10.3)
  enddo

! The following is for testing SNR calibration:
  sp3n=(s(ipk-1)+s(ipk)+s(ipk+1))               !Sig + 3*noise
  base=(sum(s)-sp3n)/(NFFT1-3.0)                !Noise per bin
  psig=sp3n-3*base                              !Sig only
  pnoise=(2500.0/df1)*base                      !Noise in 2500 Hz
  xsnr=db(psig/pnoise)
  xsnr=xsnr+5.0
  return
end subroutine getfc1

subroutine getfc2(c,npeaks,fs,fc1,fpks)

!  include 'wsprlf_params.f90'
  parameter (NZ=288*412)
  parameter (NSPS=288)
  parameter (N2=2*NSPS)
  parameter (NFFT1=16*NSPS)
  parameter (NH1=NFFT1/2)

  complex c(0:NZ-1)                     !Complex waveform
  complex cs(0:NZ-1)                    !For computing spectrum
  real a(5)
  real freqs(413),sp2(413),fpks(npeaks)
  integer pkloc(1)

  df=fs/NZ
  baud=fs/NSPS
  a(1)=-fc1
  a(2:5)=0.
  call twkfreq1(c,NZ,fs,a,cs)         !Mix down by fc1

! Filter, square, then FFT to get refined carrier frequency fc2.
  call four2a(cs,NZ,1,-1,1)          !To freq domain

  ia=nint(0.75*baud/df)
  cs(ia:NZ-1-ia)=0.                  !Save only freqs around fc1
!  do i=1,NZ/2
!    filt=1/(1+((i*df)**2/(0.50*baud)**2)**8)
!    cs(i)=cs(i)*filt
!    cs(NZ+1-i)=cs(NZ+1-i)*filt
!  enddo
  call four2a(cs,NZ,1,1,1)           !Back to time domain
  cs=cs/NZ
  cs=cs*cs                           !Square the data
  call four2a(cs,NZ,1,-1,1)          !Compute squared spectrum
! Find two peaks separated by baud
  pmax=0.
  fc2=0.
!  ja=nint(0.3*baud/df)
  ja=nint(0.5*baud/df)
  k=1
  sp2=0.0
  do j=-ja,ja
     f2=j*df
     ia=nint((f2-0.5*baud)/df)
     if(ia.lt.0) ia=ia+NZ
     ib=nint((f2+0.5*baud)/df)
     p=real(cs(ia))**2 + aimag(cs(ia))**2 +                        &
          real(cs(ib))**2 + aimag(cs(ib))**2
     if(p.gt.pmax) then
        pmax=p
        fc2=0.5*f2
     endif
     freqs(k)=0.5*f2
     sp2(k)=p
     k=k+1
!           write(52,1200) f2,p,db(p)
!1200       format(f10.3,2f15.3)
  enddo

  do i=1,npeaks
    pkloc=maxloc(sp2)
    ipk=pkloc(1)
    fpks(i)=freqs(ipk)
    ipk0=max(1,ipk-2)
    ipk1=min(413,ipk+2)
!    ipk0=ipk
!    ipk1=ipk
    sp2(ipk0:ipk1)=0.0
  enddo
  return
end subroutine getfc2


code:
program wspr5sim

! Generate simulated data for a 5-minute "WSPR-LF" mode.  Output is saved
! to a *.c5 or *.wav file.

  use wavhdr
  include 'wsprlf_params.f90'            !Set various constants
  parameter (NMAX=300*12000)
  type(hdr) h                            !Header for .wav file
  character arg*12,fname*16
  character msg*22,msgsent*22
  complex c0(0:NZ-1)
  complex c(0:NZ-1)
  real*8 fMHz
  integer itone(NN)
  integer*2 iwave(NMAX)                  !Generated full-length waveform

! Get command-line argument(s)
  nargs=iargc()
  if(nargs.ne.8) then
     print*,'Usage:   wspr5sim "message"       f0  DT fsp del  nwav nfiles snr'
     print*,'Example: wspr5sim "K1ABC FN42 30" 50 0.0 0.1 1.0  1      10   -33'
     go to 999
  endif
  call getarg(1,msg)                     !Message to be transmitted
  call getarg(2,arg)
  read(arg,*) f0                         !Freq relative to WSPR-band center (Hz)
  call getarg(3,arg)
  read(arg,*) xdt                        !Time offset from nominal (s)
  call getarg(4,arg)
  read(arg,*) fspread                    !Watterson frequency spread (Hz)
  call getarg(5,arg)
  read(arg,*) delay                      !Watterson delay (ms)
  call getarg(6,arg)
  read(arg,*) nwav                       !1 for *.wav file, 0 for *.c5 file
  call getarg(7,arg)
  read(arg,*) nfiles                     !Number of files
  call getarg(8,arg)
  read(arg,*) snrdb                      !SNR_2500

  twopi=8.0*atan(1.0)
  fs=12000.0/NDOWN                       !Sample rate
  dt=1.0/fs                              !Sample interval (s)
  tt=NSPS*dt                             !Duration of "itone" symbols (s)
  ts=2*NSPS*dt                           !Duration of OQPSK symbols (s)
  baud=1.0/tt                            !Keying rate for "itone" symbols (baud)
  txt=NZ*dt                              !Transmission length (s)
  bandwidth_ratio=2500.0/(fs/2.0)
  sig=sqrt(bandwidth_ratio) * 10.0**(0.05*snrdb)
  if(snrdb.gt.90.0) sig=1.0
  txt=NN*NSPS0/12000.0

  call genwspr5(msg,msgsent,itone)       !Encode the message, get itone
  write(*,1000) f0,xdt,txt,snrdb,fspread,delay,nfiles,msgsent
1000 format('f0:',f9.3,'   DT:',f6.2,'   txt:',f6.1,'   SNR:',f6.1,    &
          '   fspread:',f6.1,'   delay:',f6.1,'  nfiles:',i3,2x,a22)

  dphi0=twopi*(f0-0.25*baud)*dt
  dphi1=twopi*(f0+0.25*baud)*dt
  phi=0.0
  c0=0.
  k=-1 + nint(xdt/dt)
  do j=1,NN                              !Generate OQPSK waveform from itone
     dphi=dphi0
     if(itone(j).eq.1) dphi=dphi1
     if(k.eq.0) phi=-dphi
     do i=1,NSPS
        k=k+1
        phi=phi+dphi
        if(phi.gt.twopi) phi=phi-twopi
        xphi=phi
        if(k.ge.0 .and. k.lt.NZ) c0(k)=cmplx(cos(xphi),sin(xphi))
     enddo
  enddo

  call sgran()
  do ifile=1,nfiles
     c=c0
     if(nwav.eq.0) then
        if( fspread .ne. 0.0 .or. delay .ne. 0.0 ) then
           call watterson(c,NZ,fs,delay,fspread)
        endif
        c=c*sig
        if(snrdb.lt.90) then
           do i=0,NZ-1                   !Add gaussian noise at specified SNR
              xnoise=gran()
              ynoise=gran()
              c(i)=c(i) + cmplx(xnoise,ynoise)
           enddo
        endif
        write(fname,1100) ifile
1100    format('000000_',i4.4,'.c5')
        open(10,file=fname,status='unknown',access='stream')
        fMHz=10.1387d0
        nmin=5
        write(10) fname,nmin,fMHz,c      !Save to *.c5 file
        close(10)
     else
        call wspr5_wav(baud,xdt,f0,itone,snrdb,iwave)
        h=default_header(12000,NMAX)
        write(fname,1102) ifile
1102    format('000000_',i4.4,'.wav')
        open(10,file=fname,status='unknown',access='stream')
        write(10) h,iwave                !Save to *.wav file
        close(10)
     endif
     write(*,1110) ifile,xdt,f0,snrdb,fname
1110 format(i4,f7.2,f8.2,f7.1,2x,a16)
  enddo

999 end program wspr5sim


code:
subroutine wspr_fsk8_downsample(iwave,c)

! Input: i*2 data in iwave() at sample rate 12000 Hz
! Output: Complex data in c(), sampled at 12000/24=500 Hz

  include 'wspr_fsk8_params.f90'
  integer*2 iwave(NMAX)
  complex c(0:NMAXD-1)
  complex c1(0:NMAXD-1)
  complex cx(0:NMAX/2)
  real x(NMAX)
  equivalence (x,cx)

  df=12000.0/NMAX
  x=iwave
  call four2a(x,NMAX,1,-1,0)             !r2c FFT to freq domain
  i0=nint(1500.0/df)
  c1(0)=cx(i0)
  do i=1,NMAXD/2
     c1(i)=cx(i0+i)
     c1(NMAXD-i)=cx(i0-i)
  enddo
  c=c1/NMAXD
  call four2a(c,NMAXD,1,1,1)            !c2c FFT back to time domain

  return
end subroutine wspr_fsk8_downsample


code:
! LDPC (300,60) code
parameter (NDOWN=24)                  !Downsample factor
parameter (KK=60)                     !Information bits (50 + CRC10)
parameter (ND=100)                    !Data symbols
parameter (NS=14)                     !Sync symbols (2 @ Costas 7x7)
parameter (NN=NS+ND)                  !Total symbols (114)
parameter (NSPS0=24576)               !Samples per symbol at 12000 S/s
parameter (NSPS=NSPS0/NDOWN)          !Sam/sym, downsampled (1024)
parameter (N7=7*NSPS)                 !Samples in Costas 7x7 array (7168)
parameter (NZ=NSPS*NN)                !Samples in downsampled waveform (116,736)
parameter (NMAX=240*12000)            !Samples in iwave()
parameter (NMAXD=NMAX/24)             !Samples in c(), after downsampling
parameter (NFFT1=4*NSPS,NH1=NFFT1/2)
parameter (NH2=NSPS/2)


code:
program wspr_fsk8_sim

! Generate simulated data for a 4-minute "WSPR-LF" mode using 8-FSK.
! Output is saved to a *.wav file.

  use wavhdr
  include 'wspr_fsk8_params.f90'         !Set various constants
  type(hdr) h                            !Header for .wav file
  character arg*12,fname*16
  character msg*22,msgsent*22
  complex c0(0:NZ-1)
  complex c(0:NZ-1)
  real*8 fMHz
  integer itone(NN)
  integer*2 iwave(NMAX)                  !Generated full-length waveform

! Get command-line argument(s)
  nargs=iargc()
  if(nargs.ne.8) then
     print*,'Usage:   wspr5sim "message"       f0  DT fsp del  nwav nfiles snr'
     print*,'Example: wspr5sim "K1ABC FN42 30" 50 0.0 0.1 1.0  1      10   -33'
     go to 999
  endif
  call getarg(1,msg)                     !Message to be transmitted
  call getarg(2,arg)
  read(arg,*) f0                         !Freq relative to WSPR-band center (Hz)
  call getarg(3,arg)
  read(arg,*) xdt                        !Time offset from nominal (s)
  call getarg(4,arg)
  read(arg,*) fspread                    !Watterson frequency spread (Hz)
  call getarg(5,arg)
  read(arg,*) delay                      !Watterson delay (ms)
  call getarg(6,arg)
  read(arg,*) nwav                       !1 for *.wav file, 0 for *.c4 file
  call getarg(7,arg)
  read(arg,*) nfiles                     !Number of files
  call getarg(8,arg)
  read(arg,*) snrdb                      !SNR_2500

  twopi=8.0*atan(1.0)
  fs=12000.0/NDOWN                       !Sample rate after downsampling
  dt=1.0/fs                              !Sample interval (s)
  tt=NSPS*dt                             !Duration of symbols (s)
  baud=1.0/tt                            !Keying rate
  bw=8*baud
  txt=NZ*dt                              !Transmission length (s)
  bandwidth_ratio=2500.0/(fs/2.0)
  sig=sqrt(bandwidth_ratio) * 10.0**(0.05*snrdb)
  if(snrdb.gt.90.0) sig=1.0
  txt=NN*NSPS0/12000.0

  call genwspr_fsk8(msg,msgsent,itone)       !Encode the message, get itone
  write(*,1000) f0,xdt,txt,snrdb,bw,msgsent
1000 format('f0:',f9.3,'   DT:',f6.2,'   TxT:',f6.1,'   SNR:',f6.1,    &
          '  BW:',f4.1,2x,a22)

  phi=0.0
  c0=0.
  k=-1 + nint(xdt/dt)
  do j=1,NN                              !Generate OQPSK waveform from itone
     dphi=twopi*(f0+itone(j)*baud)*dt
     if(k.eq.0) phi=-dphi
     do i=1,NSPS
        k=k+1
        phi=phi+dphi
        if(phi.gt.twopi) phi=phi-twopi
        xphi=phi
        if(k.ge.0 .and. k.lt.NZ) c0(k)=cmplx(cos(xphi),sin(xphi))
     enddo
  enddo

  call sgran()
  do ifile=1,nfiles
    if(nwav.eq.0) then
      c=c0
      if( fspread .ne. 0.0 .or. delay .ne. 0.0 ) then
        call watterson(c,NZ,fs,delay,fspread)
      endif
      c=c*sig
      if( snrdb.lt.90) then
        do i=0,NZ-1
          xnoise=gran()
          ynoise=gran()
          c(i)=c(i)+cmplx(xnoise,ynoise)
        enddo
      endif
      write(fname,1100) ifile
1100  format('000000_',i4.4,'.c4')
      open(10,file=fname,status='unknown',access='stream')
      fMHz=1.866d0
      nmin=4
      write(10) fname,nmin,fMHz,c
    else
      call wspr_fsk8_wav(baud,xdt,f0,itone,snrdb,iwave)
      h=default_header(12000,NMAX)
      write(fname,1102) ifile
1102  format('000000_',i4.4,'.wav')
      open(10,file=fname,status='unknown',access='stream')
      write(10) h,iwave                !Save to *.wav file
      close(10)
    endif
    write(*,1110) ifile,xdt,f0,snrdb,fname
1110 format(i4,f7.2,f8.2,f7.1,2x,a16)
  enddo

999 end program wspr_fsk8_sim


code:
subroutine wspr_fsk8_wav(baud,xdt,f0,itone,snrdb,iwave)

! Generate iwave() from itone().

  include 'wspr_fsk8_params.f90'
  integer itone(NN)
  integer*2 iwave(NMAX)
  real*8 twopi,dt,dphi,phi
  real dat(NMAX)

  twopi=8.d0*atan(1.d0)
  dt=1.d0/12000.d0

  dat=0.
  if(snrdb.lt.90) then
     do i=1,NMAX
        dat(i)=gran()          !Generate gaussian noise, rms = 1.0
     enddo
     bandwidth_ratio=2500.0/6000.0
     sig=sqrt(2*bandwidth_ratio)*10.0**(0.05*snrdb)
  else
     sig=1.0
  endif

  phi=0.d0
  k=nint((xdt+1.0)/dt)
  do j=1,NN
     dphi=twopi*(f0+ itone(j)*baud)*dt
     if(k.eq.0) phi=-dphi
     do i=1,NSPS0
        k=k+1
        phi=phi+dphi
        if(phi.gt.twopi) phi=phi-twopi
        xphi=phi
        if(k.gt.0 .and. k.le.NMAX) dat(k)=dat(k) + sig*sin(xphi)
     enddo
  enddo
  fac=32767.0
  rms=100.0
  if(snrdb.ge.90.0) iwave=nint(fac*dat)
  if(snrdb.lt.90.0) iwave=nint(rms*dat)

  return
end subroutine wspr_fsk8_wav


code:
program wspr_fsk8d

! WSPR-LF is a potential WSPR-like mode intended for use at LF and MF.
! This version uses 4-minute T/R sequences, an LDPC (300,60) code,
! 8-FSK modulation, and noncoherent demodulation.  This decoder reads
! data from *.wav files.

! Reception and Demodulation algorithm:
!   1. Compute coarse spectrum; find fc1 = approx carrier freq
!   2. Mix from fc1 to 0; LPF at +/- 0.75*R
!   3. Find two 7x7 Costas arrays to get xdt and fc2
!   4. Mix from fc2 to 0, compute aligned symbol spectra
!   5. Get soft bits from symbol spectra

! Still to do: find and decode more than one signal in the specified passband.

  include 'wspr_fsk8_params.f90'
  character arg*8,message*22,cbits*50,infile*80,fname*16,datetime*11
  character*120 data_dir
  complex csync(0:N7-1)                 !Sync symbols for Costas 7x7 array
  complex c1(0:2*N7-1)
  complex c2(0:2*N7-1)
  complex c(0:NMAXD-1)                     !Complex waveform
  real*8 fMHz
  real rxdata(3*ND),llr(3*ND)           !Soft symbols
  real a(5)                             !For twkfreq1
  real s(0:NH2,NN)
  real savg(0:NH2)
  real ss(0:N7)
  real ss0(0:N7)
  real ps(0:7)
  integer ihdr(11)
  integer*2 iwave(NMAX)                 !Generated full-length waveform
  integer*1 idat(7)
  integer*1 decoded(KK),apmask(3*ND),cw(3*ND)
  integer icos7(0:6)
  data icos7/2,5,6,0,4,1,3/             !Costas 7x7 tone pattern

  nargs=iargc()
  if(nargs.lt.7) then
     print*,'Usage:   wspr_fsk8d -d db -f fMHz -a data_dir file1 [file2 ...]'
     go to 999
  endif
  call getarg(1,arg)
  if(arg(1:3).ne.'-d ') go to 999
  call getarg(2,arg)
  read(arg,*) degrade
  rxbw=3000.

  call getarg(3,arg)
  if(arg(1:3).ne.'-f ') go to 999
  call getarg(4,arg)
  read(arg,*) fMHz

  call getarg(5,arg)
  if(arg(1:3).ne.'-a ') go to 999
  call getarg(6,data_dir)

  open(13,file=trim(data_dir)//'/ALL_WSPR.TXT',status='unknown',   &
       position='append')

  twopi=8.0*atan(1.0)
  fs=NSPS*12000.0/NSPS0                  !Sample rate after downsampling (Hz)
  dt=1.0/fs                              !Sample interval (s)
  ts=NSPS*dt                             !Symbol duration (s)
  baud=1.0/ts                            !Keying rate (Hz)
  txt=NZ*dt                              !Transmission length (s)

  phi=0.
  k=-1
  do j=0,6
     dphi=twopi*baud*icos7(j)*dt
     do i=1,NSPS
        phi=phi+dphi
        if(phi.gt.twopi) phi=phi-twopi
        k=k+1
        csync(k)=cmplx(cos(phi),sin(phi))
     enddo
  enddo

  do ifile=7,nargs
     call getarg(ifile,infile)
     open(10,file=infile,status='old',access='stream')
     j1=index(infile,'.c4')
     j2=index(infile,'.wav')
     if(j1.gt.0) then
       read(10,end=999) fname,ntrmin,fMHz,c(0:NZ-1)
       read(fname(8:11),*) nutc
       write(datetime,'(i11)') nutc
     else if(j2.gt.0) then
       read(10,end=999) ihdr,iwave
       read(infile(j2-4:j2-1),*) nutc
       datetime=infile(j2-11:j2-1)
       if(degrade.gt.0.0) call degrade_snr(iwave,NMAX,degrade,rxbw)
       call wspr_fsk8_downsample(iwave,c)
     else
       print*,'Wrong file format?'
       go to 999
     endif
     close(10)
     pmax=0.
     df1=fs/(2*N7)
     ia=nint(100.0/df1)
     ib=nint(150.0/df1)
     ipk=0
     jpk=0
! Get xdt and f0 from the sync arrays
     do j0=0,1000,50
        j0b=j0+107*NSPS
        c1(0:N7-1)=c(j0:j0+N7-1)*conjg(csync)
        c1(N7:)=0.
        c2(0:N7-1)=c(j0b:j0b+N7-1)*conjg(csync)
        c2(N7:)=0.
        call four2a(c1,2*N7,1,-1,1)
        call four2a(c2,2*N7,1,-1,1)
        do i=0,N7
           p=1.e-9*(real(c1(i))**2 + aimag(c1(i))**2 +                        &
                real(c2(i))**2 + aimag(c2(i))**2)
           ss(i)=p
        enddo
        do i=ia,ib
           p=ss(i)
           if(p.gt.pmax) then
              pmax=p
              ipk=i
              jpk=j0
           endif
        enddo
        if(jpk.eq.j0) ss0=ss
     enddo
     f0=ipk*df1
     xdt=jpk*dt - 1.0

     sp3n=(ss0(ipk-1)+ss0(ipk)+ss0(ipk+1))      !Sig + 3*noise
     call pctile(ss0,N7,45,base)
     psig=sp3n-3*base                           !Sig only
     pnoise=(2500.0/df1)*base                   !Noise in 2500 Hz
     xsnr=db(psig/pnoise)                       !SNR from sync tones

     if(jpk.ge.0) c(0:NMAXD-1-jpk)=c(jpk:NMAXD-1)

     a(1)=-f0
     a(2:5)=0.
     call twkfreq1(c,NZ,fs,a,c)              !Mix from f0 down to 0
     call spec8(c,s,savg)                    !Get symbol spectra

     do j=1,ND
        k=j+7
        ps=s(0:7,k)
        ps=sqrt(ps)                                !### ??? ###
!        ps=log(ps)
        r1=max(ps(1),ps(3),ps(5),ps(7))-max(ps(0),ps(2),ps(4),ps(6))
        r2=max(ps(2),ps(3),ps(6),ps(7))-max(ps(0),ps(1),ps(4),ps(5))
        r4=max(ps(4),ps(5),ps(6),ps(7))-max(ps(0),ps(1),ps(2),ps(3))
        rxdata(3*j-2)=r4
        rxdata(3*j-1)=r2
        rxdata(3*j)=r1
     enddo

     rxav=sum(rxdata)/ND
     rx2av=sum(rxdata*rxdata)/ND
     rxsig=sqrt(rx2av-rxav*rxav)
     rxdata=rxdata/rxsig
     s0=1.1
     llr=2.0*rxdata/(s0*s0)
     apmask=0
     max_iterations=40
     ifer=0
     call bpdecode300(llr,apmask,max_iterations,decoded,niterations,cw)
     if(niterations.lt.0) call osd300(llr,apmask,5,decoded,cw,nhardmin,dmin)
     if(nhardmin.ge.0) niterations=nhardmin
     nbadcrc=0
     if(niterations.ge.0) call chkcrc10(decoded,nbadcrc)
     if(niterations.lt.0 .or. nbadcrc.ne.0) ifer=1
     nsnr=nint(xsnr)
!     freq=fMHz + 1.d-6*f0
     freq=1.d-6*(f0+1500)
     nfdot=0
     message='                      '
     if(ifer.eq.0) then
        write(cbits,1100) decoded(1:50)
1100    format(50i1)
        read(cbits,1102) idat
1102    format(6b8,b2)
        idat(7)=ishft(idat(7),6)
        call wqdecode(idat,message,itype)
        write(*,1112) datetime(8:11),nsnr,xdt,freq,nfdot,message
1112    format(a4,i4,f5.1,f11.6,i3,2x,a22)
     endif
     write(13,1110) datetime,0,nsnr,xdt,freq,message,nfdot
1110 format(a11,2i4,f6.2,f12.7,2x,a22,i3)
  enddo                                   ! ifile loop
  write(*,1120)
1120 format("<DecodeFinished>")

999 end program wspr_fsk8d


code:
subroutine wspr_wav(baud,xdt,h,f0,itone,snrdb,iwave)

! Generate iwave() from itone().

  include 'wspr_params.f90'
  integer itone(NN)
  integer*2 iwave(NMAX)
  real*8 twopi,dt,dphi0,dphi1,dphi,phi
  real dat(NMAX)

  twopi=8.d0*atan(1.d0)
  dt=1.d0/12000.d0
  baud=375.0/256.0

  dat=0.
  if(snrdb.lt.90) then
     do i=1,NMAX
        dat(i)=gran()          !Generate gaussian noise
     enddo
     bandwidth_ratio=2500.0/6000.0
     sig=sqrt(2*bandwidth_ratio)*10.0**(0.05*snrdb)
  else
     sig=1.0
  endif

  phi=0.d0
  k=nint(xdt/dt)
  do j=1,NN
     dphi=twopi*(f0+h*(itone(j)-1.5)*baud)*dt
     do i=1,NSPS0
        k=k+1
        phi=mod(phi+dphi,twopi)
        if(k.gt.0 .and. k.le.NMAX) dat(k)=dat(k) + sig*sin(phi)
     enddo
  enddo

  rms=100.0
  if(snrdb.lt.90.0) then
    dat=rms*dat;
    if(maxval(abs(dat)).gt.32767.0) print*,"Warning - data will be clipped."
  else
    datpk=maxval(abs(dat))
    fac=32767.9/datpk
    dat=fac*dat
  endif
  iwave=nint(dat)

  return
end subroutine wspr_wav


code:
parameter (KK=64)                     !Information bits (50 + CRC14) ?
parameter (ND=200)                    !Data symbols: LDPC (204,68), r=1/3, don't send last 4 bits
parameter (NS=16)                     !Sync symbols (16)
parameter (NN=NS+ND)                  !Total symbols (216)

parameter (NSPS0=6400)                !Samples per symbol at 12000 S/s

parameter (NDOWN=32)                  !Downsample to 200 sa/symbol (375 Hz) for candidate selection
parameter (NSPS=NSPS0/NDOWN)          !Samples per symbol
parameter (NZ=NSPS*NN)                !Samples in baseband waveform

parameter (NZ0=NSPS0*NN)              !Samples in waveform at 12000 S/s
parameter (NFFT1=4*NSPS,NH1=NFFT1/2)


code:
subroutine wsprcpm_wav(baud,xdt,h,f0,itone,snrdb,iwave)

! Generate iwave() from itone().

  include 'wsprcpm_params.f90'
  parameter (NMAX=120*12000)
  integer itone(NN)
  integer*2 iwave(NMAX)
  real*8 twopi,dt,dphi0,dphi1,dphi,phi
  real dat(NMAX)

  twopi=8.d0*atan(1.d0)
  dt=1.d0/12000.d0

  dat=0.
  if(snrdb.lt.90) then
     do i=1,NMAX
        dat(i)=gran()          !Generate gaussian noise
     enddo
     bandwidth_ratio=2500.0/6000.0
     sig=sqrt(2*bandwidth_ratio)*10.0**(0.05*snrdb)
  else
     sig=1.0
  endif

  phi=0.d0
  k=nint(xdt/dt)
  do j=1,NN
     dphi=twopi*(f0+itone(j)*(h/2.0d0)*baud)*dt
     do i=1,NSPS0
        k=k+1
        phi=phi+dphi
        if(phi.gt.twopi) phi=phi-twopi
        xphi=phi
        if(k.gt.0 .and. k.le.NMAX) dat(k)=dat(k) + sig*sin(xphi)
     enddo
  enddo
  fac=32767.0
  rms=100.0
  if(snrdb.ge.90.0) iwave=nint(fac*dat)
  if(snrdb.lt.90.0) iwave=nint(rms*dat)

  return
end subroutine wsprcpm_wav


code:
program wsprcpmd

! Decode WSPRCPM data read from *.c2 or *.wav files.

! WSPRCPM is a WSPR-like mode based on full-response CPM.
!
! Currently configured to use (204,68) r=1/3 LDPC code, regular column weight 3.
!   50 data bits + 14 bit CRC + 4 "0" bits. The 4 "0" bits are unused bits that
!   are not transmitted. At the decoder, these bits are treated as "AP" bits.
!   This shortens the code to (200,64) r=0.32, slightly decreasing the code rate.
!
! Frame format is:
! d100 p32 d100 (232)  channel symbols
!
  use crc
  include 'wsprcpm_params.f90'
  parameter(NMAX=120*12000)
  character arg*8,message*22,cbits*50,infile*80,fname*16,datetime*11
  character*22 decodes(100)
  character*120 data_dir
  character*32 uwbits
  character*68 dmsg
  complex c2(0:120*12000/32-1)              !Complex waveform
  complex cframe(0:216*200-1)               !Complex waveform
  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 csum,cterm
  real*8 fMHz
  real rxdata(ND),llr(204)               !Soft symbols
  real sbits(216),sbits1(216),sbits3(216)
  real ps(0:8191),psbest(0:8191)
  real candidates(100,2)
  integer iuniqueword0
  integer isync(200)                     !Unique word
  integer isync2(216)
  integer ipreamble(16)                 !Preamble vector
  integer ihdr(11)
  integer*2 iwave(NMAX)                 !Generated full-length waveform
  integer*1,target ::  idat(9)
  integer*1 decoded(68),apmask(204),cw(204)
  integer*1 hbits(216),hbits1(216),hbits3(216)
  data ipreamble/1,1,1,1,0,0,0,0,0,0,0,0,1,1,1,1/
  data iuniqueword0/z'30C9E8AD'/

  write(uwbits,'(b32.32)') iuniqueword0
  read(uwbits,'(32i1)') isync(1:32)
  read(uwbits,'(32i1)') isync(33:64)
  read(uwbits,'(32i1)') isync(65:96)
  read(uwbits,'(32i1)') isync(97:128)
  read(uwbits,'(32i1)') isync(129:160)
  read(uwbits,'(32i1)') isync(161:192)
  read(uwbits,'(8i1)') isync(193:200)

  fs=12000.0/NDOWN                       !Sample rate
  dt=1.0/fs                              !Sample interval (s)
  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)
  twopi=8.0*atan(1.0)

  isync2(1:100)=isync(1:100)
  isync2(101:104)=0  ! This is *not* backwards.
  isync2(105:112)=1
  isync2(113:116)=0
  isync2(117:216)=isync(101:200)

! data sync tone
!   0    0   0
!   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
       dphi0=-3*dphi
       dphi1=+1*dphi
    else
       dphi0=-1*dphi
       dphi1=+3*dphi
    endif
    phi0=0.0
    phi1=0.0
    do i=0,9
      c1(i,j)=cmplx(cos(phi1),sin(phi1))
      c0(i,j)=cmplx(cos(phi0),sin(phi0))
      phi1=mod(phi1+dphi1,twopi)
      phi0=mod(phi0+dphi0,twopi)
    enddo
  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')

  xs1=0.0
  xs2=0.0
  fr1=0.0
  fr2=0.0
  nav=0
  ngood=0

  do ifile=iarg,nargs
    call getarg(ifile,infile)
    open(10,file=infile,status='old',access='stream')
    j1=index(infile,'.c2')
    j2=index(infile,'.wav')
    if(j1.gt.0) then
       read(10,end=999) fname,ntrmin,fMHz,c2
       read(fname(8:11),*) nutc
       write(datetime,'(i11)') nutc
    else if(j2.gt.0) then
       read(10,end=999) ihdr,iwave
       read(infile(j2-4:j2-1),*) nutc
       datetime=infile(j2-11:j2-1)
       call wsprcpm_downsample(iwave,c2)
    else
       print*,'Wrong file format?'
       go to 999
    endif
    close(10)

    fa=-100.0
    fb=100.0
    fs=12000.0/32.0
    npts=120*12000.0/32.0
    nsync=16
    call getcandidate2(c2,npts,fs,fa,fb,ncand,candidates)         !First approx for freq
ndecodes=0
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)
      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)
    call noncoherent_frame_sync(c2,fcest,isync2,istart,xf2)
    write(*,'(i5,i5,i5,4(f11.5,2x))') ifile,iii,istart,fc0,fc1,fcest

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)

    s2=sum(cd*conjg(cd))/(10*216)
    cd=cd/sqrt(s2)

    do nseq=1,7
      if( nseq.eq.1 ) then  ! noncoherent single-symbol detection
        sbits1=0.0
        do ibit=1,216
          j=isync2(ibit)
          ib=(ibit-1)*10
          ccor(1,ibit)=sum(cd(ib:ib+9)*conjg(c1(0:9,j)))
          ccor(0,ibit)=sum(cd(ib:ib+9)*conjg(c0(0:9,j)))
          sbits1(ibit)=abs(ccor(1,ibit))-abs(ccor(0,ibit))
          hbits1(ibit)=0
          if(sbits1(ibit).gt.0) hbits1(ibit)=1
        enddo
        sbits=sbits1
        hbits=hbits1
        sbits3=sbits1
        hbits3=hbits1
      elseif( nseq.ge.2 ) then
        ps=0
        if( nseq.eq. 2 ) nbit=3
        if( nseq.eq. 3 ) nbit=5
        if( nseq.eq. 4 ) nbit=7
        if( nseq.eq. 5 ) nbit=9
        if( nseq.eq. 6 ) nbit=11
        if( nseq.eq. 7 ) nbit=13
        numseq=2**(nbit)
        do ibit=nbit/2+1,216-nbit/2
          ps=0.0
          pmax=0.0
          do iseq=0,numseq-1
            csum=0.0
            cterm=1.0
            k=1
            do i=nbit-1,0,-1
              ibb=iand(iseq/(2**i),1)
              csum=csum+ccor(ibb,ibit-(nbit/2+1)+k)*cterm
              cterm=-cterm
              k=k+1
            enddo
            ps(iseq)=abs(csum)
            if( ps(iseq) .gt. pmax ) then
              pmax=ps(iseq)
              ibflag=1
            endif
          enddo
          if( ibflag .eq. 1 ) then
            psbest=ps
            ibflag=0
          endif
          call getmetric2(2**(nbit/2),psbest,numseq,sbits3(ibit))
          hbits3(ibit)=0
          if(sbits3(ibit).gt.0) hbits3(ibit)=1
        enddo
        sbits=sbits3
        hbits=hbits3
      endif
      rxdata(1:100)=sbits(1:100)
      rxdata(101:200)=sbits(117:216);
      rxav=sum(rxdata(1:200))/200.0
      rx2av=sum(rxdata(1:200)*rxdata(1:200))/200.0
      rxsig=sqrt(rx2av-rxav*rxav)
      rxdata=rxdata/rxsig
      sigma=0.90
      llr(201:204)=-5.0
      llr(1:200)=2*rxdata/(sigma*sigma)
      apmask=0
      apmask(201:204)=1
      max_iterations=40
      ifer=0
      call bpdecode204(llr,apmask,max_iterations,decoded,cw,nharderror,niterations)
      nhardmin=-1
      if(nharderror.lt.0) call osd204(llr,apmask,5,decoded,cw,nhardmin,dmin)
      if(sum(decoded).eq.0) cycle
      if(nhardmin.ge.0 .or. nharderror.ge.0) then
        idat=0
        write(dmsg,'(68i1)') decoded
        read(dmsg(1:50),'(6b8.8,b2.2)') idat(1:7)
        idat(7)=idat(7)*64
        read(dmsg(51:64),'(b14.14)') ndec_crc
        ncalc_crc=iand(crc14(c_loc(idat),9),z'FFFF')
        nbadcrc=1
        if(ncalc_crc .eq. ndec_crc) nbadcrc=0
      else
        cycle
      endif
      if( nbadcrc.eq.0 ) then
        write(cbits,1200) decoded(1:50)
1200    format(50i1)
        read(cbits,1202) idat
1202    format(8b8,b4)
        idat(7)=ishft(idat(7),6)
        call wqdecode(idat,message,itype)
        idupe=0
        do i=1,ndecodes
          if(decodes(i).eq.message) idupe=1
        enddo
        if(idupe.eq.1) goto 888
        ndecodes=ndecodes+1
        decodes(ndecodes)=message
        nsnr=nint(xsnr)
        freq=fMHz + 1.d-6*(fc1+fbest)
        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
1212    format(a4,i4,f5.1,f11.6,i3,2x,a22,a1,i5,i5,i5,i5,i5)
        goto 888
      endif
    enddo ! nseq
enddo !jitter
888 continue
enddo !candidate list
  enddo !files

  write(*,1120)
1120 format("<DecodeFinished>")

999 end program wsprcpmd

subroutine coherent_preamble_fsync(c2,ipreamble,nsync,nsps,istart,fc0,fc1,xmax)
    complex c2(0:120*12000/32-1)
    complex cpreamble(0:16*200-1)
    complex ctmp1(0:4*16*200-1)
    complex ctwkp(0:16*200-1)
    complex ccohp(0:15)
    integer ipreamble(nsync)
    logical first/.true./
    save dt,first,h,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
       dphi=twopi*baud*0.5*h*dt
       do i=1,16
          dp=dphi
          if(ipreamble(i).eq.0) dp=-dphi
          do j=1,200
             cpreamble(k)=cmplx(cos(phi),sin(phi))
             phi=mod(phi+dp,twopi)
             k=k+1
          enddo
       enddo
       first=.false.
    endif

    dphi=twopi*fc0*dt
    ctwkp=cmplx(0.0,0.0)
    phi=0
    do i=0,nsync*nsps-1
      ctwkp(i)=cpreamble(i)*cmplx(cos(phi),sin(phi))
      phi=mod(phi+dphi,twopi)
    enddo
    ipstart=istart+100*200
    ctmp1=0.0
    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
    ctmp1=cshift(ctmp1,-200)
    dfp=1/(4*6400.0/12000.0*16)
    do i=150,250
      xa=abs(ctmp1(i))
      if(xa.gt.xmax) then
        ishift=i
        xmax=xa
      endif
    enddo
    delta=(ishift-200)*dfp
    xm1=abs(ctmp1(ishift-1))
    x0=abs(ctmp1(ishift))
    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
    return
end subroutine coherent_preamble_fsync

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
! noncoherent wspr-type dt estimation
      dp0=twopi*(fc-1.5*h*baud)*dt
      dp1=twopi*(fc-0.5*h*baud)*dt
      dp2=twopi*(fc+0.5*h*baud)*dt
      dp3=twopi*(fc+1.5*h*baud)*dt
      th0=0.0
      th1=0.0
      th2=0.0
      th3=0.0
      do i=0,199
        ct0(i)=cmplx(cos(th0),sin(th0))
        ct1(i)=cmplx(cos(th1),sin(th1))
        ct2(i)=cmplx(cos(th2),sin(th2))
        ct3(i)=cmplx(cos(th3),sin(th3))
        th0=mod(th0+dp0,twopi)
        th1=mod(th1+dp1,twopi)
        th2=mod(th2+dp2,twopi)
        th3=mod(th3+dp3,twopi)
      enddo
      ss=0.0
      totp=0.0
      do is=1,216
        i0=izero+it+(is-1)*200
        p0=abs(sum(c2(i0:i0+199)*conjg(ct0)))
        p1=abs(sum(c2(i0:i0+199)*conjg(ct1)))
        p2=abs(sum(c2(i0:i0+199)*conjg(ct2)))
        p3=abs(sum(c2(i0:i0+199)*conjg(ct3)))
        p0=sqrt(p0)
        p1=sqrt(p1)
        p2=sqrt(p2)
        p3=sqrt(p3)

        totp=totp+p0+p1+p2+p3
!        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
      enddo
      ss=ss/totp
      if(ss.gt.ssmax) then
        ioffset=it
        ssmax=ss
      endif
   enddo
   istart=izero+ioffset
   return
end subroutine noncoherent_frame_sync

subroutine getmetric2(ib,ps,ns,xmet)
  real ps(0:ns-1)
  xm1=0
  xm0=0
  do i=0,ns-1
    if( iand(i/ib,1) .eq. 1 .and. ps(i) .gt. xm1 ) xm1=ps(i)
    if( iand(i/ib,1) .eq. 0 .and. ps(i) .gt. xm0 ) xm0=ps(i)
  enddo
  xmet=xm1-xm0
  return
end subroutine getmetric2

subroutine downsample2(ci,f0,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)
  fs=12000.0/32.0
  df=fs/NI
  ct=ci
  call four2a(ct,NI,1,-1,1)             !c2c FFT to freq domain
  i0=nint(f0/df)
  ct=cshift(ct,i0)
  co=0.0
  co(0)=ct(0)
  b=8.0
  do i=1,NO/2
     arg=(i*df/b)**2
     filt=exp(-arg)
     co(i)=ct(i)*filt
     co(NO-i)=ct(NI-i)*filt
  enddo
  co=co/NO
  call four2a(co,NO,1,1,1)            !c2c FFT back to time domain
  return
end subroutine downsample2

subroutine getcandidate2(c,npts,fs,fa,fb,ncand,candidates)
  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 wsprcpm_downsample(iwave,c)

! Input: i*2 data in iwave() at sample rate 12000 Hz
! Output: Complex data in c(), sampled at 400 Hz

  include 'wsprcpm_params.f90'
  parameter (NMAX=120*12000,NFFT2=NMAX/32)
  integer*2 iwave(NMAX)
  complex c(0:NMAX/32-1)
  complex c1(0:NFFT2-1)
  complex cx(0:NMAX/2)
  real x(NMAX)
  equivalence (x,cx)

  df=12000.0/NMAX
  x=iwave
  call four2a(x,NMAX,1,-1,0)             !r2c FFT to freq domain
  i0=nint(1500.0/df)
  c1(0)=cx(i0)
  do i=1,NFFT2/2
     c1(i)=cx(i0+i)
     c1(NFFT2-i)=cx(i0-i)
  enddo
  c1=c1/NFFT2
  call four2a(c1,NFFT2,1,1,1)            !c2c FFT back to time domain
  c=c1(0:NMAX/32-1)
  return
end subroutine wsprcpm_downsample


code:
program wsprcpmsim

! Generate simulated data for a 2-minute "WSPR-LF" mode.  Output is saved
! to a *.c2 or *.wav file.

  use wavhdr
  include 'wsprcpm_params.f90'            !Set various constants
  parameter (NMAX=120*12000)
  type(hdr) hwav                          !Header for .wav file
  character arg*12,fname*16
  character msg*22,msgsent*22
  complex c0(0:NMAX/NDOWN-1)
  complex c(0:NMAX/NDOWN-1)
  real*8 fMHz
  integer itone(NN)
  integer*2 iwave(NMAX)                  !Generated full-length waveform

! Get command-line argument(s)
  nargs=iargc()
  if(nargs.ne.8) then
     print*,'Usage:   wsprmsksim "message"       f0  DT fsp del  nwav nfiles snr'
     print*,'Example: wsprmsksim "K1ABC FN42 30" 50 0.0 0.1 1.0  1      10   -33'
     go to 999
  endif
  call getarg(1,msg)                     !Message to be transmitted
  call getarg(2,arg)
  read(arg,*) f0                         !Freq relative to WSPR-band center (Hz)
  call getarg(3,arg)
  read(arg,*) xdt                        !Time offset from nominal (s)
  call getarg(4,arg)
  read(arg,*) fspread                    !Watterson frequency spread (Hz)
  call getarg(5,arg)
  read(arg,*) delay                      !Watterson delay (ms)
  call getarg(6,arg)
  read(arg,*) nwav                       !1 for *.wav file, 0 for *.c2 file
  call getarg(7,arg)
  read(arg,*) nfiles                     !Number of files
  call getarg(8,arg)
  read(arg,*) snrdb                      !SNR_2500

  twopi=8.0*atan(1.0)
  fs=12000.0/NDOWN                       !
  dt=1.0/fs                              !Sample interval (s)
  tt=NSPS*dt                             !Duration of "itone" symbols (s)
  baud=1.0/tt                            !Keying rate for "itone" symbols (baud)
  txt=NZ*dt                              !Transmission length (s)
  bandwidth_ratio=2500.0/(fs/2.0)
  sig=sqrt(bandwidth_ratio) * 10.0**(0.05*snrdb)
  if(snrdb.gt.90.0) sig=1.0
  txt=NN*NSPS0/12000.0

  call genwsprcpm(msg,msgsent,itone)       !Encode the message, get itone
  write(*,1000) f0,xdt,txt,snrdb,fspread,delay,nfiles,msgsent
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
     dp=twopi*(f0+itone(j)*(h/2.0)*baud)*dt
     do i=1,NSPS
        k=k+1
        phi=mod(phi+dp,twopi)
        if(k.ge.0 .and. k.lt.NMAX/NDOWN) c0(k)=cmplx(cos(phi),sin(phi))
     enddo
  enddo
  call sgran()
  do ifile=1,nfiles
     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)
        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
              xnoise=gran()
              ynoise=gran()
              c(i)=c(i) + cmplx(xnoise,ynoise)
           enddo
        endif
        write(fname,1100) ifile
1100    format('000000_',i4.4,'.c2')
        open(10,file=fname,status='unknown',access='stream')
        fMHz=10.1387d0
        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)
        write(fname,1102) ifile
1102    format('000000_',i4.4,'.wav')
        open(10,file=fname,status='unknown',access='stream')
        write(10) hwav,iwave                !Save to *.wav file
        close(10)
     endif
     write(*,1110) ifile,xdt,f0,snrdb,fname
1110 format(i4,f7.2,f8.2,f7.1,2x,a16)
  enddo

999 end program wsprcpmsim


code:
program wsprdpskd

! Decode WSPRDPSK data read from *.c2 or *.wav files.

! Currently configured to use (204,68) r=1/3 LDPC code, regular column weight 3.
!   50 data bits + 14 bit CRC + 4 "0" bits. The 4 "0" bits are unused bits that
!   are not transmitted. At the decoder, these bits are treated as "AP" bits.
!   This shortens the code to (200,64) r=0.32, slightly decreasing the code rate.
! Frame format is:
! d100 p32 d100 (232)  channel symbols
!
  use crc
  include 'wsprdpsk_params.f90'
  parameter(NMAX=120*12000)
  character arg*8,message*22,cbits*50,infile*80,fname*16,datetime*11
  character*22 decodes(100)
  character*120 data_dir
  character*32 uwbits
  character*68 dmsg
  complex c2(0:120*12000/30-1)              !Complex waveform
  complex cframe(0:232*200-1)               !Complex waveform
  complex cd(0:240*10-1)                   !Complex waveform
  complex cs(0:240)
  complex c1(0:9,0:1),c0(0:9,0:1)
  complex ccor(0:1,232)
  complex csum,cterm
  real*8 fMHz
  real rxdata(ND),llr(204)               !Soft symbols
  real sbits(232),sbits1(232),sbits3(232)
  real ps(0:8191),psbest(0:8191)
  real candidates(100,2)
  integer iuniqueword0
  integer isync(200)                     !Unique word
  integer isync2(232)
  integer ipreamble(16)                 !Preamble vector
  integer ihdr(11)
  integer*2 iwave(NMAX)                 !Generated full-length waveform
  integer*1,target ::  idat(9)
  integer*1 decoded(68),apmask(204),cw(204)
  integer*1 hbits(232),hbits1(232),hbits3(232)
  integer*1 b(14),bbest(14)
  data ipreamble/1,1,1,1,0,0,0,0,0,0,0,0,1,1,1,1/
  data iuniqueword0/z'30C9E8AD'/

  write(uwbits,'(b32.32)') iuniqueword0
  read(uwbits,'(32i1)') isync(1:32)
  read(uwbits,'(32i1)') isync(33:64)
  read(uwbits,'(32i1)') isync(65:96)
  read(uwbits,'(32i1)') isync(97:128)
  read(uwbits,'(32i1)') isync(129:160)
  read(uwbits,'(32i1)') isync(161:192)
  read(uwbits,'(8i1)') isync(193:200)

  fs=12000.0/NDOWN                       !Sample rate
  dt=1.0/fs                              !Sample interval (s)
  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)
  twopi=8.0*atan(1.0)

  isync2(1:100)=isync(1:100)
  isync2(101:104)=0  ! This is *not* backwards.
  isync2(105:112)=1
  isync2(113:116)=0
  isync2(117:216)=isync(101:200)

  dphi=twopi*baud*(h/2.0)*dt*20  ! dt*20 is samp interval after downsample
  do j=0,1
    if(j.eq.0) then
       dphi0=-3*dphi
       dphi1=+1*dphi
    else
       dphi0=-1*dphi
       dphi1=+3*dphi
    endif
    phi0=0.0
    phi1=0.0
    do i=0,9
      c1(i,j)=cmplx(cos(phi1),sin(phi1))
      c0(i,j)=cmplx(cos(phi0),sin(phi0))
      phi1=mod(phi1+dphi1,twopi)
      phi0=mod(phi0+dphi0,twopi)
    enddo
  enddo

  nargs=iargc()
  if(nargs.lt.1) then
     print*,'Usage:   wsprdpskd [-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')

  xs1=0.0
  xs2=0.0
  fr1=0.0
  fr2=0.0
  nav=0
  ngood=0

  do ifile=iarg,nargs
    call getarg(ifile,infile)
    open(10,file=infile,status='old',access='stream')
    j1=index(infile,'.c2')
    j2=index(infile,'.wav')
    if(j1.gt.0) then
       read(10,end=999) fname,ntrmin,fMHz,c2
       read(fname(8:11),*) nutc
       write(datetime,'(i11)') nutc
    else if(j2.gt.0) then
       read(10,end=999) ihdr,iwave
       read(infile(j2-4:j2-1),*) nutc
       datetime=infile(j2-11:j2-1)
       call wsprdpsk_downsample(iwave,c2)
    else
       print*,'Wrong file format?'
       go to 999
    endif
    close(10)

    fa=-10.0
    fb=10.0
    fs=12000.0/30.0
    npts=120*12000.0/30.0
!    call getcandidate2(c2,npts,fs,fa,fb,ncand,candidates)         !First approx for freq
    ncand=1
    candidates(1,1)=0.0
    candidates(1,2)=-28
    ndecodes=0
    do icand=1,ncand
      fc0=candidates(icand,1)
      xsnr=candidates(icand,2)
      call dsdpsk(c2,fc0,cd)
      i0=40
      do i=0,231
        cs(i)=cd(i0+10*i)/5e4
      enddo
! 2-bit differential detection
      do i=1,231
!      do i=1,232
        sbits(i)=-real(cs(i)*conjg(cs(i-1))) !2 symbol dpsk
!        sbits(i)=real(cs(i-1)) !for coherent dpsk
!        sbits(i)=real(cs(i))   !for coherent bpsk
      enddo
!      do i=1,231
!        sbits3(i)=-sbits(i+1)*sbits(i)  ! for coherent dpsk
!      enddo

! detect a differentially encoded block of symbols using the
! Divsalar and Simon approach, except that we estimate only
! the central symbol of the block and then step the block forward
! by one symbol.
!
      sbits3=sbits
!goto 100
      nbit=13  ! number of decoded bits to be derived from nbit+1 symbols
      numseq=2**nbit
      il=(nbit+1)/2
      ih=231-nbit/2
      do isym=il,ih
        rmax=-1e32
        b=0
        do iseq=0,numseq-1
          do i=1,nbit
            b(i)=merge(1,0,iand(iseq,2**(nbit-i))>0)
          enddo
          b(1:nbit)=2*b(1:nbit)-1
          i1=isym-(nbit+1)/2
          csum=cs(i1)
          do i=1,nbit
            bb=1
            do m=1,i
              bb=bb*b(m)
            enddo
            csum=csum+bb*cs(i1+i)
          enddo
!          ps(iseq)=abs(csum)**2
          ps(iseq)=abs(csum)
          if(ps(iseq).gt.rmax) then
            bbest=b
            rmax=ps(iseq)
          endif
        enddo
        if(isym.eq.il) then
          do i=1,isym-1
            call getmetric(2**(nbit-i),ps,numseq,xmet)
            sbits3(i)=-xmet
          enddo
        endif
        call getmetric(2**((nbit-1)/2),ps,numseq,xmet)
        sbits3(isym)=-xmet
        if(isym.eq.ih) then
          do i=ih+1,231
            call getmetric(2**(231-i),ps,numseq,xmet)
            sbits3(i)=-xmet
          enddo
        endif
      enddo
100 continue
      rxdata(1:100)=sbits3(1:100)
      rxdata(101:200)=sbits3(132:231);
      rxav=sum(rxdata(1:200))/200.0
      rx2av=sum(rxdata(1:200)*rxdata(1:200))/200.0
      rxsig=sqrt(rx2av-rxav*rxav)
      rxdata=rxdata/rxsig
      sigma=0.90
      llr(201:204)=-5.0
      llr(1:200)=2*rxdata/(sigma*sigma)
      apmask=0
      apmask(201:204)=1
      max_iterations=40
      ifer=0
      call bpdecode204(llr,apmask,max_iterations,decoded,cw,nharderror,niterations)
      nhardmin=-1
      if(nharderror.lt.0) call osd204(llr,apmask,5,decoded,cw,nhardmin,dmin)
      if(sum(decoded).eq.0) cycle
      if(nhardmin.ge.0 .or. nharderror.ge.0) then
        idat=0
        write(dmsg,'(68i1)') decoded
        read(dmsg(1:50),'(6b8.8,b2.2)') idat(1:7)
        idat(7)=idat(7)*64
        read(dmsg(51:64),'(b14.14)') ndec_crc
        ncalc_crc=iand(crc14(c_loc(idat),9),z'FFFF')
        nbadcrc=1
        if(ncalc_crc .eq. ndec_crc) nbadcrc=0
      else
        cycle
      endif
      if( nbadcrc.eq.0 ) then
        write(cbits,1200) decoded(1:50)
1200    format(50i1)
        read(cbits,1202) idat
1202    format(8b8,b4)
        idat(7)=ishft(idat(7),6)
        call wqdecode(idat,message,itype)
        idupe=0
        do i=1,ndecodes
          if(decodes(i).eq.message) idupe=1
        enddo
        if(idupe.eq.1) goto 888
        ndecodes=ndecodes+1
        decodes(ndecodes)=message
        nsnr=nint(xsnr)
        freq=fMHz + 1.d-6*(fc1+fbest)
        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
1212    format(a4,i4,f5.1,f11.6,i3,2x,a22,a1,i5,i5,i5,i5,i5)
        goto 888
      endif
888 continue
    enddo !candidate list
  enddo !files

  write(*,1120)
1120 format("<DecodeFinished>")

999 end program wsprdpskd

subroutine getmetric(ib,ps,ns,xmet)
  real ps(0:ns-1)
  xm1=-1e32
  xm0=-1e32
  do i=0,ns-1
    if( iand(i/ib,1) .eq. 1 .and. ps(i) .gt. xm1 ) then
      xm1=ps(i)
    endif
    if( iand(i/ib,1) .eq. 0 .and. ps(i) .gt. xm0 ) then
      xm0=ps(i)
    endif
  enddo
  xmet=xm1-xm0
  return
end subroutine getmetric

subroutine getmetric3(ib,ps,ns,xmet)
  real ps(0:ns-1)
  xm1=0
  xm0=0
  do i=0,ns-1
    if( iand(i/ib,1) .eq. 1 ) then
      xm1=xm1+ps(i)
    endif
    if( iand(i/ib,1) .eq. 0 ) then
      xm0=xm0+ps(i)
    endif
  enddo
  xmet=xm1-xm0
  return
end subroutine getmetric3

subroutine dsdpsk(ci,f0,co)
  parameter(NI=240*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)

  pi=4.0*atan(1.0)
  fs=12000.0/30.0
  df=fs/NI
  ct=ci
  call four2a(ct,NI,1,-1,1)             !c2c FFT to freq domain
  i0=nint(f0/df)
  ct=cshift(ct,i0)
  co=0.0
  co(0)=ct(0)

  dt=20/fs
  beta=1.0
  tt=10*dt
  baud=1/tt
  bw=(1+beta)*baud/2.0
  bf=(1-beta)*baud/2.0
  iw=bw/df
  if=bf/df
  co=0.0
  co(0)=ct(0)
  do i=1,iw
     filt=((1.0+cos(pi*(i-if)/(iw-if)))/2.0)**0.5
     co(i)=ct(i)*filt
     co(NO-i)=ct(NI-i)*filt
  enddo
  co=co/NO
  call four2a(co,NO,1,1,1)            !c2c FFT back to time domain
  return
end subroutine dsdpsk

subroutine getcandidate2(c,npts,fs,fa,fb,ncand,candidates)
  parameter(NDAT=200,NFFT1=120*12000/30,NH1=NFFT1/2,NFFT2=120*12000/300,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 wsprdpsk_downsample(iwave,c)

! Input: i*2 data in iwave() at sample rate 12000 Hz
! Output: Complex data in c(), sampled at 400 Hz

  include 'wsprdpsk_params.f90'
  parameter (NMAX=120*12000,NFFT2=NMAX/30)
  integer*2 iwave(NMAX)
  complex c(0:NMAX/30-1)
  complex c1(0:NFFT2-1)
  complex cx(0:NMAX/2)
  real x(NMAX)
  equivalence (x,cx)

  df=12000.0/NMAX
  x=iwave
  call four2a(x,NMAX,1,-1,0)             !r2c FFT to freq domain
  i0=nint(1500.0/df)
  c1(0)=cx(i0)
  do i=1,NFFT2/2
     c1(i)=cx(i0+i)
     c1(NFFT2-i)=cx(i0-i)
  enddo
  c1=c1/NFFT2
  call four2a(c1,NFFT2,1,1,1)            !c2c FFT back to time domain
  c=c1(0:NMAX/30-1)
  return
end subroutine wsprdpsk_downsample


code:
program wsprdpsksim

! Generate simulated data for a 2-minute "WSPR-DPSK" mode.  Output is saved
! to a *.c2 or *.wav file.

  use wavhdr
  include 'wsprdpsk_params.f90'            !Set various constants
  parameter (NMAX=120*12000)
  type(hdr) hwav                          !Header for .wav file
  character arg*12,fname*16
  character msg*22,msgsent*22
  complex c0(0:NMAX/NDOWN-1)
  complex c(0:NMAX/NDOWN-1)
  complex c0wav(0:NMAX-1)
  complex cwav(0:NMAX-1)
  real*8 fMHz
  integer imessage(NN)
  integer*2 iwave(NMAX)                  !Generated full-length waveform

! Get command-line argument(s)
  nargs=iargc()
  if(nargs.ne.8) then
     print*,'Usage:   wsprdpsksim "message"       f0  DT fsp del  nwav nfiles snr'
     print*,'Example: wsprdpsksim "K1ABC FN42 30" 50 0.0 0.1 1.0  1      10   -33'
     go to 999
  endif
  call getarg(1,msg)                     !Message to be transmitted
  call getarg(2,arg)
  read(arg,*) f0                         !Freq relative to WSPR-band center (Hz)
  call getarg(3,arg)
  read(arg,*) xdt                        !Time offset from nominal (s)
  call getarg(4,arg)
  read(arg,*) fspread                    !Watterson frequency spread (Hz)
  call getarg(5,arg)
  read(arg,*) delay                      !Watterson delay (ms)
  call getarg(6,arg)
  read(arg,*) nwav                       !1 for *.wav file, 0 for *.c2 file
  call getarg(7,arg)
  read(arg,*) nfiles                     !Number of files
  call getarg(8,arg)
  read(arg,*) snrdb                      !SNR_2500

  twopi=8.0*atan(1.0)
  pi=twopi/2.0
  fs=12000.0/NDOWN
  dt=1.0/fs                              !Sample interval (s)
  tt=NSPS*dt                             !Duration of "itone" symbols (s)
  baud=1.0/tt                            !Keying rate for "itone" symbols (baud)
  txt=NZ*dt                              !Transmission length (s)
  bandwidth_ratio=2500.0/(fs/2.0)
  sig=sqrt(bandwidth_ratio) * 10.0**(0.05*snrdb)
  if(snrdb.gt.90.0) sig=1.0
  txt=NN*NSPS0/12000.0

  call genwsprdpsk(msg,msgsent,imessage)       !Encode the message, get itone
  imessage=2*imessage-1
  write(*,1000) f0,xdt,txt,snrdb,fspread,delay,nfiles,msgsent
1000 format('f0:',f9.3,'   DT:',f6.2,'   txt:',f6.1,'   SNR:',f6.1,    &
          '   fspread:',f6.1,'   delay:',f6.1,'  nfiles:',i3,2x,a22)

  beta=1.0   ! excess bandwidth
  if(nwav.eq.0) then
     df=fs/(NMAX/NDOWN)                      !
     c=0
     bw=(1+beta)*baud/2.0
     bf=(1-beta)*baud/2.0
     iw=bw/df
     if=bf/df
     c(0:if-1)=1.0
     if(iw.gt.if) then
       do i=if,iw
         c(i)=((1.0+cos(pi*(i-if)/(iw-if)))/2.0)**0.5
       enddo
     endif
     c(NMAX/NDOWN-1:NMAX/NDOWN-iw:-1)=c(1:iw)

     istart=xdt/dt
     c0=0.0
     do i=1,NN
       c0(istart+(i-1)*200)=imessage(i)
     enddo
     call four2a(c0,NMAX/NDOWN,1,1,1)
     c0=c0*conjg(c)
     ic=f0/df
     c0=cshift(c0,ic)
     call four2a(c0,NMAX/NDOWN,1,-1,1)
     xx=sum(abs(c0(istart:istart+NN*200-1)**2))/(NN*200)
     c0=c0/sqrt(xx)

     call sgran()
     do ifile=1,nfiles
       c=c0
       if( fspread .ne. 0.0 .or. delay .ne. 0.0 ) then
          call watterson(c,NMAX/NDOWN,fs,delay,fspread)
       endif
       c=c*sig
       if(snrdb.lt.90) then
          do i=0,NMAX/NDOWN-1                   !Add gaussian noise at specified SNR
             xnoise=gran()
             ynoise=gran()
             c(i)=c(i) + cmplx(xnoise,ynoise)
          enddo
       endif
snrtest=sum(abs(c(istart:istart+NN*200-1)**2))/(NN*200)/2.0-1.0
write(*,*) 'sample SNR: ',10*log10(snrtest)+10*log10(0.4/2.5)
       write(fname,1100) ifile
1100   format('000000_',i4.4,'.c2')
       open(10,file=fname,status='unknown',access='stream')
       fMHz=10.1387d0
       nmin=2
       write(10) fname,nmin,fMHz,c      !Save to *.c2 file
       close(10)
     enddo
  else
     fs=12000.0
     df=fs/NMAX
     dt=1/fs
     bandwidth_ratio=2500.0/(fs/2.0)
     sig=sqrt(2*bandwidth_ratio) * 10.0**(0.05*snrdb)
     if(snrdb.gt.90.0) sig=1.0
     cwav=0
     bw=(1+beta)*baud/2.0
     bf=(1-beta)*baud/2.0
     iw=bw/df
     if=bf/df
     cwav(0:if-1)=1.0
     if(iw.gt.if) then
       do i=if,iw
         cwav(i)=((1.0+cos(pi*(i-if)/(iw-if)))/2.0)**0.5
       enddo
     endif
     cwav(NMAX-1:NMAX-iw:-1)=cwav(1:iw)

     istart=xdt/dt
     c0wav=0.0
     do i=1,NN
       c0wav(istart+(i-1)*200*NDOWN)=imessage(i)
     enddo
     call four2a(c0wav,NMAX,1,1,1)
     c0wav=c0wav*conjg(cwav)
     ic=f0/df
     c0wav=cshift(c0wav,-ic)
     call four2a(c0wav,NMAX,1,-1,1)
     xx=sum(abs(c0wav(istart:istart+NN*200*NDOWN-1))**2)/(NN*200*NDOWN)
     c0wav=c0wav/sqrt(xx)
write(*,*) 'Peak power: ',maxval(abs(c0wav)**2)
write(*,*) 'Average power: ',sum(abs(c0wav(istart:istart+NN*200*NDOWN-1))**2)/(NN*200*NDOWN)
     call sgran()
     do ifile=1,nfiles
       cwav=c0wav
       if( fspread .ne. 0.0 .or. delay .ne. 0.0 ) then
          call watterson(cwav,NMAX,fs,delay,fspread)
       endif
       cwav=cwav*sig
       if(snrdb.lt.90) then
          do i=1,NMAX                   !Add gaussian noise at specified SNR
             xnoise=gran()
             iwave(i)=100*(real(cwav(i-1)) + xnoise)
          enddo
       endif
snrtest=sum(real(iwave(istart:istart+NN*200*NDOWN-1)**2)/(NN*200*NDOWN))/100.0**2-1
write(*,*) 'sample SNR: ',10*log10(snrtest)+10*log10(6.0/2.5)
       hwav=default_header(12000,NMAX)
       write(fname,1102) ifile
1102   format('000000_',i4.4,'.wav')
       open(10,file=fname,status='unknown',access='stream')
       write(10) hwav,iwave                !Save to *.wav file
       close(10)
     enddo
  endif
  write(*,1110) ifile,xdt,f0,snrdb,fname
1110 format(i4,f7.2,f8.2,f7.1,2x,a16)

999 end program wsprdpsksim


code:
program wsprlf

  parameter (NN=121)                    !Total symbols
  parameter (NSPS=28800)                  !Samples per symbol @ fs=12000 Hz
  parameter (NZ=NSPS*NN)                !Samples in waveform

  character*8 arg
  complex c(0:NZ-1)
  real*8 twopi,fs,f0,dt,phi,dphi
  real x(0:NZ-1)
  real p(0:NZ/2)
  real h0(0:NSPS/2)                     !Pulse shape, rising edge
  real h1(0:NSPS/2)                     !Pulse shape, trailing edge
  real tmp(NN)
  integer id(NN)                        !Generated data
  integer ie(NN)                        !Differentially encoded data
  data fs/12000.d0/

  nargs=iargc()
  if(nargs.ne.3) then
     print*,'Usage: wsprlf f0 t1 snr'
     goto 999
  endif
  call getarg(1,arg)
  read(arg,*) f0
  call getarg(2,arg)
  read(arg,*) t1
  call getarg(3,arg)
  read(arg,*) snrdb

  call random_number(tmp)          !Generate random bipolar data
  id=1
  where(tmp.lt.0.5) id=-1
  ie(1)=1
  do i=2,NN                        !Differentially encode
     ie(i)=id(i)*ie(i-1)
  enddo

  n1=nint(t1*NSPS)
  twopi=8.d0*atan(1.d0)

  do i=0,2*n1-1                    !Define the shape functions
     if(i.le.n1-1) then
        h0(i)=0.5*(1.0-cos(0.5*i*twopi/n1))
     else
        h1(i-n1)=0.5*(1.0-cos(0.5*i*twopi/n1))
     endif
  enddo
  if(t1.eq.0.0) h0=1
  if(t1.eq.0.0) h1=1

! Shape the channel pulses
  x=1.
  x(0:n1-1)=h0(0:n1-1)           !Leading edge of 1st pulse
  do j=2,NN                      !Leading edges
     if(ie(j).ne.ie(j-1)) then
        ia=(j-1)*NSPS + 1
        ib=ia+n1-1
        x(ia:ib)=h0(0:n1-1)
     endif
  enddo
  do j=1,NN-1                    !Trailing edges
     if(ie(j+1).ne.ie(j)) then
        ib=j*NSPS
        ia=ib-n1+1
        x(ia:ib)=h1(0:n1-1)
     endif
  enddo
  ib=NN*NSPS-1
  ia=ib-n1+1
  x(ia:ib)=h1(0:n1-1)           !Trailing edge of last pulse

  dt=1.d0/fs
  ts=dt*NSPS
  baud=fs/NSPS
  write(*,1000) baud,ts
1000 format('Baud:',f6.3,'  Tsym:',f6.3)

  dphi=twopi*f0*dt
  phi=0.d0
  i=-1
  do j=1,NN                     !Generate the baseband waveform
     a=ie(j)
     do k=1,NSPS
        i=i+1
        x(i)=a*x(i)
        phi=phi+dphi
        if(phi.gt.twopi) phi=phi-twopi
        xphi=phi
        c(i)=x(i)*cmplx(cos(xphi),sin(xphi))
        sym=i*dt/ts
        if(j.le.20) write(13,1010) sym,x(i),c(i)
1010    format(4f12.6)
     enddo
  enddo

  call four2a(c,NZ,1,-1,1)      !To freq domain
  df=fs/NZ
  nh=NZ/2
  do i=0,nh
     f=i*df
     p(i)=real(c(i))**2 + aimag(c(i))**2
  enddo
  p=p/maxval(p)
  do i=0,nh                      !Save spectrum for plotting
     write(14,1020) i*df,p(i),10.0*log10(p(i)+1.e-8)
1020 format(f10.3,2e12.3)
  enddo

999 end program wsprlf


code:
!parameter (NDOWN=540)                 !Downsample factor (default 540)
parameter (NDOWN=30)                  !Downsample factor (default 540)
parameter (KK=60)                     !Information bits (50 + CRC10)
parameter (ND=300)                    !Data symbols: LDPC (300,60), r=1/5
parameter (NS=109)                    !Sync symbols (2 x 48 + Barker 13)
parameter (NR=3)                      !Ramp up/down
parameter (NN=NR+NS+ND)               !Total symbols (412)
parameter (NSPS0=8640)                !Samples per symbol at 12000 S/s
parameter (NSPS=NSPS0/NDOWN)          !Samples per MSK symbol (16)
parameter (N2=2*NSPS)                 !Samples per OQPSK symbol (32)
parameter (N13=13*N2)                 !Samples in central sync vector (416)
parameter (NZ=NSPS*NN)                !Samples in baseband waveform (6592)
parameter (NZMAX=NSPS0*NN)
parameter (NFFT1=4*NSPS,NH1=NFFT1/2)


code:
program wsprlfsim

! Simulate characteristics of a potential "WSPR-LF" mode using LDPC (300,60)
! code, OQPSK modulation, and 5 minute T/R sequences.

! Reception and Demodulation algorithm:
!   1. Compute coarse spectrum; find fc1 = approx carrier freq
!   2. Mix from fc1 to 0; LPF at +/- 0.75*R
!   3. Square, FFT; find peaks near -R/2 and +R/2 to get fc2
!   4. Mix from fc2 to 0
!   5. Fit cb13 (central part of csync) to c -> lag, phase
!   6. Fit complex ploynomial for channel equalization
!   7. Get soft bits from equalized data

  include 'wsprlf_params.f90'

! Q: Would it be better for central Sync array to use both I and Q channels?

  character*8 arg
  complex cbb(0:NZ-1)                   !Complex baseband waveform
  complex csync(0:NZ-1)                 !Sync symbols only, from cbb
  complex c(0:NZ-1)                     !Complex waveform
  complex c0(0:NZ-1)                    !Complex waveform
  complex c1(0:NZ-1)                    !Complex waveform
  complex zz(NS+ND)                     !Complex symbol values (intermediate)
  complex z
  real xnoise(0:NZ-1)                   !Generated random noise
  real ynoise(0:NZ-1)                   !Generated random noise
  real rxdata(ND),llr(ND)               !Soft symbols
  real pp(2*NSPS)                       !Shaped pulse for OQPSK
  real a(5)                             !For twkfreq1
  real aa(20),bb(20)                    !Fitted polyco's
  real t(11)
  character*12 label(11)
  integer*8 count0,count1,count2,count3,clkfreq
  integer nc(11)
  integer id(NS+ND)                     !NRZ values (+/-1) for Sync and Data
  integer ierror(NS+ND)
  integer icw(NN)
  integer itone(NN)
  integer*1 msgbits(KK),decoded(KK),apmask(ND),cw(ND)
!  integer*1 codeword(ND)
  data msgbits/0,0,1,0,0,1,1,1,1,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,1,1,0,0,0,1, &
       1,1,1,0,1,1,1,1,1,1,1,0,0,1,0,0,1,1,0,1,1,0,1,0,1,1,0,0,1,1/
  data label/'genwsprlf','twkfreq1 a','watterson','noise gen','getfc1w',    &
       'getfc2w','twkfreq1 b','xdt loop','cpolyfitw','msksoftsym',          &
       'bpdecode300'/

  nargs=iargc()
  if(nargs.ne.6) then
     print*,'Usage:   wsprlfsim f0(Hz) delay(ms) fspread(Hz) maxn iters snr(dB)'
     print*,'Example: wsprlfsim 0 0 0 5 10 -20'
     print*,'Set snr=0 to cycle through a range'
     go to 999
  endif
  call getarg(1,arg)
  read(arg,*) f0                         !Generated carrier frequency
  call getarg(2,arg)
  read(arg,*) delay                      !Delta_t (ms) for Watterson model
  call getarg(3,arg)
  read(arg,*) fspread                    !Fspread (Hz) for Watterson model
  call getarg(4,arg)
  read(arg,*) maxn                       !Max nterms for polyfit
  call getarg(5,arg)
  read(arg,*) iters                      !Iterations at each SNR
  call getarg(6,arg)
  read(arg,*) snrdb                      !Specified SNR_2500

  nc=0
  twopi=8.0*atan(1.0)
  fs=NSPS*12000.0/NSPS0                  !Sample rate = 22.2222... Hz
  dt=1.0/fs                              !Sample interval (s)
  tt=NSPS*dt                             !Duration of "itone" symbols (s)
  ts=2*NSPS*dt                           !Duration of OQPSK symbols (s)
  baud=1.0/tt                            !Keying rate for "itone" symbols (baud)
  txt=NZ*dt                              !Transmission length (s)
  bandwidth_ratio=2500.0/(fs/2.0)
  write(*,1000) fs,f0,delay,fspread,maxn,baud,3*baud,txt,iters
1000 format('fs:',f10.3,'  f0:',f5.1,'  Delay:',f4.1,'  fSpread:',f5.2,   &
          '  maxn:',i3,/'Baud:',f8.3,'  BW:',f5.1,'  TxT:',f6.1,'  iters:',i4/)
  write(*,1004)
1004 format(/'  SNR     sync   data   ser     ber    fer   fsigma  tsigma',   &
          '    tsec'/68('-'))

  do i=1,N2                              !Half-sine pulse shape
     pp(i)=sin(0.5*(i-1)*twopi/(2*NSPS))
  enddo

  t=0.
  call system_clock(count0,clkfreq)
  call genwsprlf(msgbits,id,icw,cbb,csync,itone)!Generate baseband waveform
  call system_clock(count1,clkfreq)
  t(1)=float(count1-count0)/float(clkfreq)
  nc(1)=nc(1)+1
  do i=0,NZ-1
     write(40,4001) i,cbb(i),csync(i)
4001 format(i8,4f12.6)
  enddo

  call system_clock(count0,clkfreq)
  a=0.
  a(1)=f0
  call twkfreq1(cbb,NZ,fs,a,c0)          !Mix baseband to specified frequency
  call system_clock(count1,clkfreq)
  t(2)=float(count1-count0)/float(clkfreq)
  nc(2)=nc(2)+1

  isna=-20
  isnb=-40
  if(snrdb.ne.0.0) then
     isna=nint(snrdb)
     isnb=isna
  endif
  do isnr=isna,isnb,-1                   !Loop over SNR range
     if(isna.ne.isnb) snrdb=isnr
     sig=sqrt(bandwidth_ratio) * 10.0**(0.05*snrdb)
     if(snrdb.gt.90.0) sig=1.0
     nhard=0
     nhardsync=0
     nfe=0
     sqf=0.
     sqt=0.

     call system_clock(count2,clkfreq)
     do iter=1,iters                     !Loop over requested iterations
        c=c0
write(*,*) 'iter ',iter
        call system_clock(count0,clkfreq)
        if(delay.ne.0.0 .or. fspread.ne.0.0) then
           call watterson(c,NZ,fs,delay,fspread)
        endif
        call system_clock(count1,clkfreq)
        t(3)=t(3)+float(count1-count0)/float(clkfreq)
        nc(3)=nc(3)+1

        call system_clock(count0,clkfreq)
        c=sig*c                          !Scale to requested SNR
        if(snrdb.lt.90) then
           do i=0,NZ-1                   !Generate gaussian noise
              xnoise(i)=gran()
              ynoise(i)=gran()
           enddo
           c=c + cmplx(xnoise,ynoise)    !Add AWGN noise
        endif
        call system_clock(count1,clkfreq)
        t(4)=t(4)+float(count1-count0)/float(clkfreq)
        nc(4)=nc(4)+1

        call system_clock(count0,clkfreq)
        call getfc1w(c,fs,fc1)               !First approx for freq
        call system_clock(count1,clkfreq)
        t(5)=t(5)+float(count1-count0)/float(clkfreq)
        nc(5)=nc(5)+1
write(*,*) 'fc1 ',fc1
        call system_clock(count0,clkfreq)
        call getfc2w(c,csync,fs,fc1,fc2,fc3) !Refined freq
write(*,*) 'fc1,fc2,fc3 ',fc1,fc2,fc3
        call system_clock(count1,clkfreq)
        t(6)=t(6)+float(count1-count0)/float(clkfreq)
        nc(6)=nc(6)+1
        sqf=sqf + (fc1+fc2-f0)**2

        call system_clock(count0,clkfreq)
!NB: Measured performance is about equally good using fc2 or fc3 here:
        a(1)=-(fc1+fc2)
        a(2:5)=0.
        call twkfreq1(c,NZ,fs,a,c)       !Mix c down by fc1+fc2
        call system_clock(count1,clkfreq)
        t(7)=t(7)+float(count1-count0)/float(clkfreq)
        nc(7)=nc(7)+1

! The following may not be necessary?
!        z=sum(c(3088:3503)*cb13)/208.0     !Get phase from Barker 13 vector
!        z0=z/abs(z)
!        c=c*conjg(z0)

        call system_clock(count0,clkfreq)
!---------------------------------------------------------------- DT
! Not presently used:
        amax=0.
        jpk=0
        iaa=0
        ibb=NZ-1
        do j=-20*NSPS,20*NSPS,NSPS/8
           ia=j
           ib=NZ-1+j
           if(ia.lt.0) then
              ia=0
              iaa=-j
           else
              iaa=0
           endif
           if(ib.gt.NZ-1) then
              ib=NZ-1
              ibb=NZ-1-j
           endif
           z=sum(c(ia:ib)*conjg(csync(iaa:ibb)))
           if(abs(z).gt.amax) then
              amax=abs(z)
              jpk=j
           endif
        enddo
        xdt=jpk/fs
        sqt=sqt + xdt**2
        call system_clock(count1,clkfreq)
        t(8)=t(8)+float(count1-count0)/float(clkfreq)
        nc(8)=nc(8)+1

!-----------------------------------------------------------------

        nterms=maxn
        c1=c
        do itry=1,20
           idf=itry/2
           if(mod(itry,2).eq.0) idf=-idf
           nhard0=0
           nhardsync0=0
           ifer=1
           a(1)=idf*0.00085
           a(2:5)=0.
           call system_clock(count0,clkfreq)
           call twkfreq1(c1,NZ,fs,a,c)       !Mix c1 into c
           call cpolyfitw(c,pp,id,maxn,aa,bb,zz,nhs)
           call system_clock(count1,clkfreq)
           t(9)=t(9)+float(count1-count0)/float(clkfreq)
           nc(9)=nc(9)+1

           call system_clock(count0,clkfreq)
           call msksoftsymw(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhardsync0)
           call system_clock(count1,clkfreq)
           t(10)=t(10)+float(count1-count0)/float(clkfreq)
           nc(10)=nc(10)+1

           if(nhardsync0.gt.35) cycle
           rxav=sum(rxdata)/ND
           rx2av=sum(rxdata*rxdata)/ND
           rxsig=sqrt(rx2av-rxav*rxav)
           rxdata=rxdata/rxsig
           ss=0.84
           llr=2.0*rxdata/(ss*ss)
           apmask=0
           max_iterations=40
           ifer=0
           call system_clock(count0,clkfreq)
           call bpdecode300(llr,apmask,max_iterations,decoded,niterations,cw)
           call system_clock(count1,clkfreq)
           t(11)=t(11)+float(count1-count0)/float(clkfreq)
           nc(11)=nc(11)+1
           nbadcrc=0
           if(niterations.ge.0) call chkcrc10(decoded,nbadcrc)
           if(niterations.lt.0 .or. count(msgbits.ne.decoded).gt.0 .or.     &
                nbadcrc.ne.0) ifer=1
           if(ifer.eq.0) exit
        enddo                                !Freq dither loop
        nhard=nhard+nhard0
        nhardsync=nhardsync+nhardsync0
        nfe=nfe+ifer
        if(nhardsync0+nhard0+niterations+ifer.gt.0) write(42,1045) snrdb,  &
             nhardsync0,nhard0,niterations,ifer,xdt
1045    format(f6.1,4i6,f8.2)
     enddo
     call system_clock(count3,clkfreq)
     tsec=float(count3-count2)/float(clkfreq)

     fsigma=sqrt(sqf/iters)
     tsigma=sqrt(sqt/iters)
     ser=float(nhardsync)/(NS*iters)
     ber=float(nhard)/(ND*iters)
     fer=float(nfe)/iters
     write(*,1050)  snrdb,nhardsync,nhard,ser,ber,fer,fsigma,tsigma,tsec
1050 format(f6.1,2i7,2f8.4,f7.3,2f8.2f8.3)
  enddo

  write(*,1060) NS*iters,ND*iters
1060 format(68('-')/6x,2i7)

  write(*,1065)
1065 format(/'Timing           sec      frac    calls'/39('-'))
  do i=1,11
     write(*,1070) label(i),t(i),t(i)/sum(t),nc(i)
1070 format(a12,2f9.3,i8)
  enddo
  write(*,1072) sum(t),1.0
1072 format(39('-')/12x,2f10.3)

999 end program wsprlfsim


code:
!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
!-------------------------------------------------------------------------------

program wsprsim

  use wavhdr
  include 'wspr_params.f90'
  type(hdr) hwav
  character arg*12,fname14*14,fname15*15
  character*22 msg,msgsent
  complex   c0(0:NMAX/NDOWN-1)
  complex   c(0:NMAX/NDOWN-1)
  integer   itone(NN)
  integer*2 iwave(NMAX)
  real*8    fMHz

! Get command-line argument(s)
  nargs=iargc()
  if(nargs.ne.8) then
     print*,'Usage:   wsprsim "message"       f0  DT fsp del  nwav nfiles snr'
     print*,'Example: wsprsim "K1ABC FN42 30" 50 0.0 0.1 1.0  1      10   -33'
     go to 999
  endif
  call getarg(1,msg)                     !Message to be transmitted
  call getarg(2,arg)
  read(arg,*) f0                         !Freq relative to WSPR-band center (Hz)
  call getarg(3,arg)
  read(arg,*) xdt                        !Time offset from nominal (s)
  call getarg(4,arg)
  read(arg,*) fspread                    !Watterson frequency spread (Hz)
  call getarg(5,arg)
  read(arg,*) delay                      !Watterson delay (ms)
  call getarg(6,arg)
  read(arg,*) nwav                       !1 for *.wav file, 0 for *.c2 file
  call getarg(7,arg)
  read(arg,*) nfiles                     !Number of files
  call getarg(8,arg)
  read(arg,*) snrdb                      !SNR_2500

  twopi=8.0*atan(1.0)
  fs=12000.0/NDOWN
  dt=1.0/fs
  tt=NSPS*dt
  baud=12000.0/8192.0

  txt=NZ*dt                              !Transmission length (s)
  bandwidth_ratio=2500.0/(fs/2.0)
  sig=sqrt(bandwidth_ratio) * 10.0**(0.05*snrdb)
  if(snrdb.gt.90.0) sig=1.0
  txt=NN*NSPS0/12000.0

  call genwspr(msg,msgsent,itone)       !Encode the message, get itone

  write(*,1000) f0,xdt,txt,snrdb,fspread,delay,nfiles,msgsent
1000 format('f0:',f9.3,'   DT:',f6.2,'   txt:',f6.1,'   SNR:',f6.1,    &
          '   fspread:',f6.1,'   delay:',f6.1,'  nfiles:',i3,2x,a22)
!  write(*,*) "Channel symbols: "
!  write(*,'(162i2)') itone

  h=1.0
  phi=0.0
  c0=0.
  k=-1 + nint(xdt/dt)
  do j=1,NN
     dphi=-twopi*(f0+h*(itone(j)-1.5)*baud)*dt
     do i=1,NSPS
        k=k+1
        phi=mod(phi+dphi,twopi)
        if(k.ge.0 .and. k.lt.NMAX/NDOWN) c0(k)=cmplx(cos(phi),sin(phi))
     enddo
  enddo
  call sgran()
  do ifile=1,nfiles
     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)
        endif
        c=c*sig
        if(snrdb.lt.90) then
           do i=0,NMAX/NDOWN-1                   !Add gaussian noise at specified SNR
              xnoise=gran()
              ynoise=gran()
              c(i)=c(i) + cmplx(xnoise,ynoise)
           enddo
        endif
        write(fname14,1100) ifile
1100    format('000000_',i4.4,'.c2')
        open(10,file=fname14,status='unknown',access='stream')
        fMHz=10.1387d0
        nmin=2
        write(10) fname14,nmin,fMHz,c      !Save to *.c2 file
        close(10)
        write(*,1108) ifile,xdt,f0,snrdb,fname14
1108    format(i4,f7.2,f8.2,f7.1,2x,a14)
     else
        freq=1500.0+f0
        call wspr_wav(baud,xdt,h,freq,itone,snrdb,iwave)
        hwav=default_header(12000,NMAX)
        write(fname15,1102) ifile
1102    format('000000_',i4.4,'.wav')
        open(10,file=fname15,status='unknown',access='stream')
        write(10) hwav,iwave                !Save to *.wav file
        close(10)
        write(*,1110) ifile,xdt,f0,snrdb,fname15
1110    format(i4,f7.2,f8.2,f7.1,2x,a15)
     endif
  enddo

999 end program wsprsim