code:
		Quick Start for DXpedition Mode
		-------------------------------

These notes are intended for operators already familiar with WSJT-X
and FT8 mode.  QSOs between the Dxpedition ("Fox") and other stations
("Hounds") are completed with as little as one transmission per Hound,
as in the following examples:

----------------------------------------------------------------------------
    Fox (300-600 Hz)                 Hounds
----------------------------------------------------------------------------
1.  CQ KH1DX AJ10
2.                                   KH1DX K1ABC FN42, KH1DX W9XYZ EN37, ...
3.  K1ABC KH1DX -13
4.                                   KH1DX K1ABC R-11
5.  K1ABC RR73; W9XYZ <KH1DX> -17
6.                                   KH1DX W9XYZ R-16
7.  W9XYZ RR73; G4AAA <KH1DX> -09
8.  ...
----------------------------------------------------------------------------

Everybody sets dial frequency to an agreed number and uses CAT control
with Split Operation (either *Rig* or *Fake It*).  Fox transmits up to
5 signals simultaneously, at audio frequencies 300, 360, ... 540
Hz.  Hounds make initial calls (e.g., line 2 above) anywhere in the
range 1000 - 4000 Hz.  They send "R+rpt" 350 Hz above the frequency
where Fox called them.

INSTRUCTIONS FOR FOX
--------------------

1. Start WSJT-X in FT8 mode.  Select *Fox* on the *Settings ->
Advanced tab*.  On the main window, check *Tx even/1st*, *Auto Seq*,
and *Hold Tx Freq*; uncheck *Call 1st*.  Set *Tx 300 Hz* and select
Tab 3.

2. In Fox mode the left window (called "Band Activity" in normal FT8
mode) is labeled "Stations calling DXpedition <MyCall>".  It will be
filled with a sorted list of calling Hounds.  You can sort by Call,
Grid, S/N, Distance, or Random order by using the comboBox at top
right of Tab 3.  You can limit the displayed Hound callsigns to those
no stronger than *Max dB*.  Fox might use this feature to discourage
Hounds from engaging in a QRO arms race.

3. *N Slots* sets the number of simultaneous Fox signals to be used.
Fox carries out as many as *N Slots* QSOs simultaneously.

4. *Repeats* sets the maximum number of repeat transmissions of the
same message.  A QSO is aborted when this number would be exceeded.

5. The *CQ* comboBox on Tab 3 offers a selection of directed CQ
messages.  *Reset* clears the QSO queue.

6. The Fox operator's main task is to select Hounds to be called and
worked.  The text box on Tab 3 holds the "QSO queue": a list of Hound
calls to be worked.  Hit Enter to select the top callsign from the
sorted list of callers (left window), or double-click on any
particular call.  Either actiion moves that Hound into the "QSO
queue".

7. The right window displays decodes of signals below 1000 Hz.
Normally these should include only Hound messages containing "R+rpt"
and Fox's own transmissions.

8. To get things started, toggle *Enable Tx* to red.  If a Hound call
is available in the QSO queue, that station will be called.  If the
QSO queue is empty, Fox calls CQ.

9. If you're using Nslots = 2 or higher, your signal no longer has
a constant envelope.  To avoid producing intermod sidebands you need
to ensure linearity in your Tx system.  One way to get things about right
is to use the WSJT-X *Tune* button to generate a pure tone.  Reduce the
Tx audio level until your power output decreases by 10% or so.  Use this
level for your Fox transmissions.

NOTE: If you are generating Nslots signals, the average power in each one
will be 1/Nslots^2 of its normal value for single-signal transmissions.

Nslots  Relative dB
-------------------
  1        0
  2       -6
  3      -9.5
  4      -12
  5      -14

The following features are not yet implemented for Fox:

1. Enforce all required settings
2. Tx message timeout
3. Manual abort of selected QSO
4. All Tx and Rx messages to all.txt
5. Additional sort criteria for Hound calls
6. Selectable timeout for keeping Hounds in the sorted list
7. Display number of active callers
8. Display QSO rate

INSTRUCTIONS FOR HOUND
----------------------

1. Start WSJT-X in FT8 mode.  Select *Hound* On the *Settings ->
Advanced* tab.  On the main window check *Auto Seq* and uncheck *Tx
even/1st*, *Call 1st*, and *Hold Tx Freq*.  Set *Tx nnnn Hz* to some
frequency between 1000 and 4000 Hz, and select *Tab 1*.  Enter Fox's
callsign and locator in DX Call and DX Grid, select Tx1, and start
*Monitor*.

2. When you have copied Fox, hit *Enable Tx* to call him.  You may
keep calling until he answers.  You may wish to move your TxFreq
around, hoping to find a clear calling frequency.

3. When you are called by Fox with a signal report, your next
transmission will automatically be sent as Tx3 ("R+rpt").  When Fox
receives that message he responds with "RR73", and your QSO is
complete!

The following features are not yet implemented for Hound:

1. Force all required settings
2. React properly to directed CQs from Fox
3. Disable Tx2, 4, 5, 6
4. For Tx1, enforce TxFreq >= 1000 Hz


code:
subroutine baseline(s,nfa,nfb,sbase)

! Fit baseline to spectrum (for FT8)
! Input:  s(npts)         Linear scale in power
! Output: sbase(npts)    Baseline

  implicit real*8 (a-h,o-z)
  real*4 s(1920)
  real*4 sbase(1920)
  real*4 base
  real*8 x(1000),y(1000),a(5)
  data nseg/10/,npct/10/

  df=12000.0/3840.0                    !3.125 Hz
  ia=max(1,nint(nfa/df))
  ib=nint(nfb/df)
  do i=ia,ib
     s(i)=10.0*log10(s(i))            !Convert to dB scale
  enddo

  nterms=5
  nlen=(ib-ia+1)/nseg                 !Length of test segment
  i0=(ib-ia+1)/2                      !Midpoint
  k=0
  do n=1,nseg                         !Loop over all segments
     ja=ia + (n-1)*nlen
     jb=ja+nlen-1
     call pctile(s(ja),nlen,npct,base) !Find lowest npct of points
     do i=ja,jb
        if(s(i).le.base) then
           if (k.lt.1000) k=k+1       !Save all "lower envelope" points
           x(k)=i-i0
           y(k)=s(i)
        endif
     enddo
  enddo
  kz=k
  a=0.
  call polyfit(x,y,y,kz,nterms,0,a,chisqr)  !Fit a low-order polynomial
  do i=ia,ib
     t=i-i0
     sbase(i)=a(1)+t*(a(2)+t*(a(3)+t*(a(4)+t*(a(5))))) + 0.65
!     write(51,3051) i*df,s(i),sbase(i)
!3051 format(3f12.3)
  enddo

  return
end subroutine baseline

license: mit


code:
subroutine chkcrc12a(decoded,nbadcrc)

  use crc
  integer*1 decoded(87)
  integer*1, target:: i1Dec8BitBytes(11)
  character*87 cbits

! Write decoded bits into cbits: 75-bit message plus 12-bit CRC
  write(cbits,1000) decoded
1000 format(87i1)
  read(cbits,1001) i1Dec8BitBytes
1001 format(11b8)
  read(cbits,1002) ncrc12                         !Received CRC12
1002 format(75x,b12)

  i1Dec8BitBytes(10)=iand(i1Dec8BitBytes(10),128+64+32)
  i1Dec8BitBytes(11)=0
  icrc12=crc12(c_loc(i1Dec8BitBytes),11)          !CRC12 computed from 75 msg bits
  ! icrc12=xor(icrc12, 41) ! TODO: jsherer - could change the crc here

  nbadcrc=1
  if(ncrc12.eq.icrc12) nbadcrc=0

  return
end subroutine chkcrc12a


code:
subroutine compress(c)

  parameter (NMAX=15*12000)             !Samples in iwave (180,000)
  complex c(0:NMAX-1)
  real xr(0:NMAX-1),xi(0:NMAX-1)

  xr=real(c)
  call wavestats(xr,NMAX,rms,pk,pwr_pk,pwr_ave)
  xr=xr/rms
  xi=aimag(c)/rms

  do i=0,NMAX-1
     c(i)=rms*cmplx(h1(xr(i)),h1(xi(i)))
  enddo

!  par=pwr_pk/pwr_ave
!  write(*,1010) 5,rms,pk,pwr_pk,pwr_ave,par
!1010 format(i3,2f10.3,3f10.2)
!  call wavestats(xi,NMAX,rms,pk,pwr_pk,pwr_ave)
!  par=pwr_pk/pwr_ave
!  write(*,1010) 5,rms,pk,pwr_pk,pwr_ave,par

  return
end subroutine compress

subroutine wavestats(x,kz,rms,pk,pwr_pk,pwr_ave)

  real x(kz)

  sumsq=dot_product(x,x)
  rms=sqrt(sumsq/kz)
  pk=max(maxval(x),-minval(x))
  pwr_pk=pk*pk
  pwr_ave=sumsq/kz

  return
end subroutine wavestats


code:
subroutine encode174(message,codeword)
! Encode an 87-bit message and return a 174-bit codeword.
! The generator matrix has dimensions (87,87).
! The code is a (174,87) regular ldpc code with column weight 3.
! The code was generated using the PEG algorithm.
! After creating the codeword, the columns are re-ordered according to
! "colorder" to make the codeword compatible with the parity-check matrix
!

include "ldpc_174_87_params.f90"

integer*1 codeword(N)
integer*1 gen(M,K)
integer*1 itmp(N)
integer*1 message(K)
integer*1 pchecks(M)
logical first
data first/.true./

save first,gen

if( first ) then ! fill the generator matrix
  gen=0
  do i=1,M
    do j=1,11
      read(g(i)( (j-1)*2+1:(j-1)*2+2 ),"(Z2)") istr
        do jj=1, 8
          icol=(j-1)*8+jj
          if( icol .le. 87 ) then
            if( btest(istr,8-jj) ) gen(i,icol)=1
          endif
        enddo
    enddo
  enddo
first=.false.
endif

do i=1,M
  nsum=0
  do j=1,K
    nsum=nsum+message(j)*gen(i,j)
  enddo
  pchecks(i)=mod(nsum,2)
enddo
itmp(1:M)=pchecks
itmp(M+1:N)=message(1:K)
codeword(colorder+1)=itmp(1:N)

return
end subroutine encode174


code:
subroutine extractmessage174(decoded,msgreceived,ncrcflag)
  use iso_c_binding, only: c_loc,c_size_t
  use crc
  use packjt
  character*68 alphabet
  character*22 msgreceived
  character*87 cbits
  integer*1 decoded(87)
  integer*1, target::  i1Dec8BitBytes(11)
  integer*4 i4Dec6BitWords(12)

  alphabet='0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz-+/?.'

! Write decoded bits into cbits: 75-bit message plus 12-bit CRC
  write(cbits,1000) decoded
1000 format(87i1)
  read(cbits,1001) i1Dec8BitBytes
1001 format(11b8)
  read(cbits,1002) ncrc12                         !Received CRC12
1002 format(75x,b12)

  i1Dec8BitBytes(10)=iand(i1Dec8BitBytes(10),128+64+32)
  i1Dec8BitBytes(11)=0
  icrc12=crc12(c_loc(i1Dec8BitBytes),11)          !CRC12 computed from 75 msg bits
!  icrc12=xor(icrc12, 41) ! TODO: jsherer - could change the crc here

  if(ncrc12.eq.icrc12 .or. sum(decoded(57:87)).eq.0) then  !### Kludge ###
! CRC12 checks out --- unpack 72-bit message
    do ibyte=1,12
      itmp=0
      do ibit=1,6
        itmp=ishft(itmp,1)+iand(1,decoded((ibyte-1)*6+ibit))
      enddo
      i4Dec6BitWords(ibyte)=itmp
    enddo

    !call unpackmsg(i4Dec6BitWords,msgreceived,.false.,'      ')

    msgreceived='                      '
    do ibyte=1,12
      itmp=i4Dec6BitWords(ibyte)
      msgreceived(ibyte:ibyte) = alphabet(itmp+1:itmp+1)
    enddo

    ncrcflag=1
  else
    msgreceived=' '
    ncrcflag=-1
  endif
  return
  end subroutine extractmessage174


code:
subroutine filt8(f0,nslots,width,wave)

  parameter (NFFT=180000,NH=NFFT/2)
  real wave(NFFT)
  real x(NFFT)
  real s1(0:NH)
  real s2(0:NH)
  complex cx(0:NH)
  equivalence (x,cx)

  x=wave
  call four2a(x,NFFT,1,-1,0)              !r2c
  df=12000.0/NFFT
  fa=f0 - 0.5*6.25
  fb=f0 + 7.5*6.25 + (nslots-1)*60.0
  ia2=nint(fa/df)
  ib1=nint(fb/df)
  ia1=nint(ia2-width/df)
  ib2=nint(ib1+width/df)
  pi=4.0*atan(1.0)
  do i=ia1,ia2
     fil=(1.0 + cos(pi*df*(i-ia2)/width))/2.0
     cx(i)=fil*cx(i)
  enddo
  do i=ib1,ib2
     fil=(1.0 + cos(pi*df*(i-ib1)/width))/2.0
     cx(i)=fil*cx(i)
  enddo
  cx(0:ia1-1)=0.
  cx(ib2+1:)=0.

  call four2a(cx,nfft,1,1,-1)                  !c2r
  wave=x/nfft

!###
  if(nslots.ne.99) return
  x=wave
  call four2a(x,NFFT,1,-1,0)              !r2c
  do i=0,NH
     s1(i)=real(cx(i))**2 + aimag(cx(i))**2
  enddo
  nadd=20
  call smo(s1,NH+1,s2,nadd)
  do i=0,NH
     freq=i*df
     write(29,3101) freq,db(s2(i)) - 72.0
3101 format(2f12.3)
  enddo
!###

  return
end subroutine filt8


code:
subroutine foxfilt(nslots,nfreq,width,wave)

  parameter (NN=79,ND=58,KK=87,NSPS=4*1920)
  parameter (NWAVE=NN*NSPS,NFFT=614400,NH=NFFT/2)
  real wave(NWAVE)
  real x(NFFT)
  complex cx(0:NH)
  equivalence (x,cx)

  x(1:NWAVE)=wave
  x(NWAVE+1:)=0.
  call four2a(x,NFFT,1,-1,0)              !r2c
  df=48000.0/NFFT
  fa=nfreq - 0.5*6.25
  fb=nfreq + 7.5*6.25 + (nslots-1)*60.0
  ia2=nint(fa/df)
  ib1=nint(fb/df)
  ia1=nint(ia2-width/df)
  ib2=nint(ib1+width/df)
  pi=4.0*atan(1.0)
  do i=ia1,ia2
     fil=(1.0 + cos(pi*df*(i-ia2)/width))/2.0
     cx(i)=fil*cx(i)
  enddo
  do i=ib1,ib2
     fil=(1.0 + cos(pi*df*(i-ib1)/width))/2.0
     cx(i)=fil*cx(i)
  enddo
  cx(0:ia1-1)=0.
  cx(ib2+1:)=0.

  call four2a(cx,nfft,1,1,-1)                  !c2r
  wave=x(1:NWAVE)/nfft

  return
end subroutine foxfilt


code:
subroutine foxgen()

  ! Called from MainWindow::foxTxSequencer() to generate the Tx waveform in
  ! FT8 Fox mode.  The Tx message can contain up to 5 "slots", each carrying
  ! its own FT8 signal.

  ! Encoded messages can be of the form "HoundCall FoxCall rpt" (a standard FT8
  ! message with i3bit=0) or "HoundCall_1 RR73; HoundCall_2 <FoxCall> rpt",
  ! a new message type with i3bit=1.  The waveform is generated with
  ! fsample=48000 Hz; it is compressed to reduce the PEP-to-average power ratio,
  ! with (currently disabled) filtering afterware to reduce spectral growth.

  ! Input message information is provided in character array cmsg(5), in
  ! common/foxcom/.  The generated wave(NWAVE) is passed back in the same
  ! common block.

  use crc
  parameter (NN=79,ND=58,KK=87,NSPS=4*1920)
  parameter (NWAVE=NN*NSPS,NFFT=614400,NH=NFFT/2)
  character*40 cmsg
  character*22 msg,msgsent
  character*6 mygrid
  character*87 cbits
  character*88 cb88
  logical bcontest
  integer itone(NN)
  integer icos7(0:6)
  integer*1 msgbits(KK),codeword(3*ND),msgbits2
  integer*1, target:: i1Msg8BitBytes(11)
  integer*1, target:: mycall
  real x(NFFT)
  real*8 dt,twopi,f0,fstep,dfreq,phi,dphi
  complex cx(0:NH)
  common/foxcom/wave(NWAVE),nslots,nfreq,i3bit(5),cmsg(5),mycall(12)
  common/foxcom2/itone2(NN),msgbits2(KK)
  equivalence (x,cx),(y,cy)
  data icos7/4,2,5,6,1,3,0/                   !Costas 7x7 tone pattern

  bcontest=.false.
  fstep=60.d0
  dfreq=6.25d0
  dt=1.d0/48000.d0
  twopi=8.d0*atan(1.d0)
  mygrid='      '
  irpt=0
  nplot=0
  wave=0.

  do n=1,nslots
     i3b=i3bit(n)
     if(i3b.eq.0) then
        msg=cmsg(n)(1:22)                     !Standard FT8 message
     else
        i1=index(cmsg(n),' ')                 !Special Fox message
        i2=index(cmsg(n),';')
        i3=index(cmsg(n),'<')
        i4=index(cmsg(n),'>')
        msg=cmsg(n)(1:i1)//cmsg(n)(i2+1:i3-2)//'                   '
        read(cmsg(n)(i4+2:i4+4),*) irpt
     endif
     call genft8(msg,mygrid,bcontest,0,msgsent,msgbits,itone)
!     print*,'Foxgen:',n,cmsg(n),msgsent

     if(i3b.eq.1) then
        icrc10=crc10(c_loc(mycall),12)
        nrpt=irpt+30
        write(cbits,1001) msgbits(1:56),icrc10,nrpt,i3b,0
1001    format(56b1.1,b10.10,b6.6,b3.3,b12.12)
        read(cbits,1002) msgbits
1002    format(87i1)

        cb88=cbits//'0'
        read(cb88,1003) i1Msg8BitBytes(1:11)
1003    format(11b8)
        icrc12=crc12(c_loc(i1Msg8BitBytes),11)
!        icrc12=xor(icrc12, 41) ! TODO: jsherer - could change the crc here

        write(cbits,1001) msgbits(1:56),icrc10,nrpt,i3b,icrc12
        read(cbits,1002) msgbits

        call encode174(msgbits,codeword)      !Encode the test message

! Message structure: S7 D29 S7 D29 S7
        itone(1:7)=icos7
        itone(36+1:36+7)=icos7
        itone(NN-6:NN)=icos7
        k=7
        do j=1,ND
           i=3*j -2
           k=k+1
           if(j.eq.30) k=k+7
           itone(k)=codeword(i)*4 + codeword(i+1)*2 + codeword(i+2)
        enddo
     endif

! Make copies of itone() and msgbits() for ft8sim
     itone2=itone
     msgbits2=msgbits
     f0=nfreq + fstep*(n-1)
     phi=0.d0
     k=0
     do j=1,NN
        f=f0 + dfreq*itone(j)
        dphi=twopi*f*dt
        do ii=1,NSPS
           k=k+1
           phi=phi+dphi
           xphi=phi
           wave(k)=wave(k)+sin(xphi)
        enddo
     enddo
  enddo
  kz=k

  peak1=maxval(abs(wave))
  wave=wave/peak1
!  call plotspec(1,wave)          !Plot the spectrum

! Apply compression
!  rms=sqrt(dot_product(wave,wave)/kz)
!  wave=wave/rms
!  do i=1,NWAVE
!     wave(i)=h1(wave(i))
!  enddo
!  peak2=maxval(abs(wave))
!  wave=wave/peak2

!  call plotspec(2,wave)          !Plot the spectrum

  width=50.0
  call foxfilt(nslots,nfreq,width,wave)
  peak3=maxval(abs(wave))
  wave=wave/peak3

!  nadd=1000
!  j=0
!  do i=1,NWAVE,nadd
!     sx=dot_product(wave(i:i+nadd-1),wave(i:i+nadd-1))
!     j=j+1
!     write(30,3001) j,sx/nadd
!3001 format(i8,f12.6)
!  enddo

!  call plotspec(3,wave)          !Plot the spectrum

  return
end subroutine foxgen

! include 'plotspec.f90'


code:
subroutine foxgen_wrap(msg40,msgbits,itone)

  parameter (NN=79,ND=58,KK=87,NSPS=4*1920)
  parameter (NWAVE=NN*NSPS)

  character*40 msg40,cmsg
  character*12 mycall12
  integer*1 msgbits(KK),msgbits2
  integer itone(NN)
  common/foxcom/wave(NWAVE),nslots,nfreq,i3bit(5),cmsg(5),mycall12
  common/foxcom2/itone2(NN),msgbits2(KK)

  nslots=1
  nfreq=300
  i1=index(msg40,'<')
  i2=index(msg40,'>')
  mycall12=msg40(i1+1:i2-1)//'    '
  cmsg(1)=msg40
  i3bit(1)=1
  call foxgen()
  msgbits=msgbits2
  itone=itone2

  return
end subroutine foxgen_wrap


code:
subroutine ft8_downsample(dd,newdat,f0,c1)

! Downconvert to complex data sampled at 200 Hz ==> 32 samples/symbol

  parameter (NMAX=15*12000,NSPS=1920)
  parameter (NFFT1=192000,NFFT2=3200)      !192000/60 = 3200

  logical newdat
  complex c1(0:NFFT2-1)
  complex cx(0:NFFT1/2)
  real dd(NMAX),x(NFFT1)
  equivalence (x,cx)
  save cx

  if(newdat) then
! Data in dd have changed, recompute the long FFT
     x(1:NMAX)=dd
     x(NMAX+1:NFFT1)=0.                       !Zero-pad the x array
     call four2a(cx,NFFT1,1,-1,0)             !r2c FFT to freq domain
     newdat=.false.
  endif

  df=12000.0/NFFT1
  baud=12000.0/NSPS
  i0=nint(f0/df)
  ft=f0+8.0*baud
  it=min(nint(ft/df),NFFT1/2)
  fb=f0-1.0*baud
  ib=max(1,nint(fb/df))
  k=0
  c1=0.
  do i=ib,it
   c1(k)=cx(i)
   k=k+1
  enddo
  c1=cshift(c1,i0-ib)
  call four2a(c1,NFFT2,1,1,1)            !c2c FFT back to time domain
  fac=1.0/sqrt(float(NFFT1)*NFFT2)
  c1=fac*c1

  return
end subroutine ft8_downsample


code:
! LDPC (174,87) code
parameter (KK=87)                     !Information bits (75 + CRC12)
parameter (ND=58)                     !Data symbols
parameter (NS=21)                     !Sync symbols (3 @ Costas 7x7)
parameter (NN=NS+ND)                  !Total channel symbols (79)
parameter (NSPS=1920)                 !Samples per symbol at 12000 S/s
parameter (NZ=NSPS*NN)                !Samples in full 15 s waveform (151,680)
parameter (NMAX=15*12000)             !Samples in iwave (180,000)
parameter (NFFT1=2*NSPS, NH1=NFFT1/2) !Length of FFTs for symbol spectra
parameter (NSTEP=NSPS/4)              !Rough time-sync step size
parameter (NHSYM=NMAX/NSTEP-3)        !Number of symbol spectra (1/4-sym steps)
parameter (NDOWN=60)                  !Downsample factor


code:
  parameter (MAXTEST=75,NTEST=68)
  character*40 testmsg(MAXTEST)
  character*40 testmsgchk(MAXTEST)
  ! Test msgs should include the extremes for the different types
  ! See pfx.f90
  ! Type 1 P & A
  ! Type 1 1A & E5
  data testmsg(1:NTEST)/         &
       "CQ WB9XYZ EN34",         &
       "CQ DX WB9XYZ EN34",      &
       "QRZ WB9XYZ EN34",        &
       "KA1ABC WB9XYZ EN34",     &
       "KA1ABC WB9XYZ RO",       &
       "KA1ABC WB9XYZ -21",      &
       "KA1ABC WB9XYZ R-19",     &
       "KA1ABC WB9XYZ RRR",      &
       "KA1ABC WB9XYZ 73",       &
       "KA1ABC WB9XYZ",          &
       "CQ 000 WB9XYZ EN34",     &
       "CQ 999 WB9XYZ EN34",     &
       "CQ EU WB9XYZ EN34",      &
       "CQ WY WB9XYZ EN34",      &
       "1A/KA1ABC WB9XYZ",       &
       "E5/KA1ABC WB9XYZ",       &
       "KA1ABC 1A/WB9XYZ",       &
       "KA1ABC E5/WB9XYZ",       &
       "KA1ABC/P WB9XYZ",        &
       "KA1ABC/A WB9XYZ",        &
       "KA1ABC WB9XYZ/P",        &
       "KA1ABC WB9XYZ/A",        &
       "CQ KA1ABC/P",            &
       "CQ WB9XYZ/A",            &
       "QRZ KA1ABC/P",           &
       "QRZ WB9XYZ/A",           &
       "DE KA1ABC/P",            &
       "DE WB9XYZ/A",            &
       "CQ 1A/KA1ABC",           &
       "CQ E5/KA1ABC",           &
       "DE 1A/KA1ABC",           &
       "DE E5/KA1ABC",           &
       "QRZ 1A/KA1ABC",          &
       "QRZ E5/KA1ABC",          &
       "CQ WB9XYZ/1A",           &
       "CQ WB9XYZ/E5",           &
       "QRZ WB9XYZ/1A",          &
       "QRZ WB9XYZ/E5",          &
       "DE WB9XYZ/1A",           &
       "DE WB9XYZ/E5",           &
       "CQ A000/KA1ABC FM07",    &
       "CQ ZZZZ/KA1ABC FM07",    &
       "QRZ W4/KA1ABC FM07",     &
       "DE W4/KA1ABC FM07",      &
       "CQ W4/KA1ABC -22",       &
       "DE W4/KA1ABC -22",       &
       "QRZ W4/KA1ABC -22",      &
       "CQ W4/KA1ABC R-22",      &
       "DE W4/KA1ABC R-22",      &
       "QRZ W4/KA1ABC R-22",     &
       "DE W4/KA1ABC 73",        &
       "CQ KA1ABC FM07",         &
       "QRZ KA1ABC FM07",        &
       "DE KA1ABC/VE6 FM07",     &
       "CQ KA1ABC/VE6 -22",      &
       "DE KA1ABC/VE6 -22",      &
       "QRZ KA1ABC/VE6 -22",     &
       "CQ KA1ABC/VE6 R-22",     &
       "DE KA1ABC/VE6 R-22",     &
       "QRZ KA1ABC/VE6 R-22",    &
       "DE KA1ABC 73",           &
       "HELLO WORLD",            &
       "ZL4/KA1ABC 73",          &
       "KA1ABC XL/WB9XYZ",       &
       "KA1ABC WB9XYZ/W4",       &
       "DE KA1ABC/QRP 2W",       &
       "KA1ABC/1 WB9XYZ/1",      &
       "123456789ABCDEFGH"/
  data testmsgchk(1:NTEST)/      &
       "CQ WB9XYZ EN34",         &
       "CQ DX WB9XYZ EN34",      &
       "QRZ WB9XYZ EN34",        &
       "KA1ABC WB9XYZ EN34",     &
       "KA1ABC WB9XYZ RO",       &
       "KA1ABC WB9XYZ -21",      &
       "KA1ABC WB9XYZ R-19",     &
       "KA1ABC WB9XYZ RRR",      &
       "KA1ABC WB9XYZ 73",       &
       "KA1ABC WB9XYZ",          &
       "CQ 000 WB9XYZ EN34",     &
       "CQ 999 WB9XYZ EN34",     &
       "CQ EU WB9XYZ EN34",      &
       "CQ WY WB9XYZ EN34",      &
       "1A/KA1ABC WB9XYZ",       &
       "E5/KA1ABC WB9XYZ",       &
       "KA1ABC 1A/WB9XYZ",       &
       "KA1ABC E5/WB9XYZ",       &
       "KA1ABC/P WB9XYZ",        &
       "KA1ABC/A WB9XYZ",        &
       "KA1ABC WB9XYZ/P",        &
       "KA1ABC WB9XYZ/A",        &
       "CQ KA1ABC/P",            &
       "CQ WB9XYZ/A",            &
       "QRZ KA1ABC/P",           &
       "QRZ WB9XYZ/A",           &
       "DE KA1ABC/P",            &
       "DE WB9XYZ/A",            &
       "CQ 1A/KA1ABC",           &
       "CQ E5/KA1ABC",           &
       "DE 1A/KA1ABC",           &
       "DE E5/KA1ABC",           &
       "QRZ 1A/KA1ABC",          &
       "QRZ E5/KA1ABC",          &
       "CQ WB9XYZ/1A",           &
       "CQ WB9XYZ/E5",           &
       "QRZ WB9XYZ/1A",          &
       "QRZ WB9XYZ/E5",          &
       "DE WB9XYZ/1A",           &
       "DE WB9XYZ/E5",           &
       "CQ A000/KA1ABC FM07",    &
       "CQ ZZZZ/KA1ABC FM07",    &
       "QRZ W4/KA1ABC FM07",     &
       "DE W4/KA1ABC FM07",      &
       "CQ W4/KA1ABC -22",       &
       "DE W4/KA1ABC -22",       &
       "QRZ W4/KA1ABC -22",      &
       "CQ W4/KA1ABC R-22",      &
       "DE W4/KA1ABC R-22",      &
       "QRZ W4/KA1ABC R-22",     &
       "DE W4/KA1ABC 73",        &
       "CQ KA1ABC FM07",         &
       "QRZ KA1ABC FM07",        &
       "DE KA1ABC/VE6 FM07",     &
       "CQ KA1ABC/VE6 -22",      &
       "DE KA1ABC/VE6 -22",      &
       "QRZ KA1ABC/VE6 -22",     &
       "CQ KA1ABC/VE6 R-22",     &
       "DE KA1ABC/VE6 R-22",     &
       "QRZ KA1ABC/VE6 R-22",    &
       "DE KA1ABC 73",           &
       "HELLO WORLD",            &
       "ZL4/KA1ABC 73",          &
       "KA1ABC XL/WB9",          &
       "KA1ABC WB9XYZ",          &
       "DE KA1ABC/QRP",          &
       "KA1ABC/1 WB9X",          &
       "123456789ABCD"/


code:
subroutine ft8apset(mycall12,mygrid6,hiscall12,hisgrid6,bcontest,apsym,iaptype)
  parameter(NAPM=4,KK=87)
  character*12 mycall12,hiscall12
  character*22 msg,msgsent
  character*6 mycall,hiscall
  character*6 mygrid6,hisgrid6
  character*4 hisgrid
  logical bcontest
  integer apsym(KK)
  integer*1 msgbits(KK)
  integer itone(KK)

!  mycall=mycall12(1:6)
!  hiscall=hiscall12(1:6)
!  hisgrid=hisgrid6(1:4)
!  if(len_trim(hiscall).eq.0) then
!    iaptype=1
!    hiscall="K9AN"
!  else
!    iaptype=2
!  endif
!  hisgrid=hisgrid6(1:4)
!!  if(len_trim(hisgrid).eq.0) hisgrid="EN50"
!  if(index(hisgrid," ").eq.0) hisgrid="EN50"
  msg='tZQpZP-slh4+'                            ! HELLO WORLD
  i3bit=0                                       ! ### TEMPORARY ??? ###
  call genft8(msg,mygrid6,bcontest,i3bit,msgsent,msgbits,itone)
  apsym=2*msgbits-1

  return
end subroutine ft8apset


code:
subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,lapcqonly,   &
     napwid,lsubtract,nagain,iaptype,mycall12,mygrid6,hiscall12,bcontest,    &
     sync0,f1,xdt,xbase,apsym,nharderrors,dmin,nbadcrc,ipass,iera,msg37,xsnr)

  use crc
  use timer_module, only: timer
  include 'ft8_params.f90'
  parameter(NP2=2812)
  character*37 msg37
  character message*22,msgsent*22,origmsg*22
  character*12 mycall12,hiscall12
  character*6 mycall6,mygrid6,hiscall6,c1,c2
  character*87 cbits
  logical bcontest
  real a(5)
  real s1(0:7,ND),s2(0:7,NN),s1sort(8*ND)
  real ps(0:7),psl(0:7)
  real bmeta(3*ND),bmetb(3*ND),bmetap(3*ND)
  real llr(3*ND),llra(3*ND),llr0(3*ND),llr1(3*ND),llrap(3*ND)           !Soft symbols
  real dd0(15*12000)
  integer*1 decoded(KK),decoded0(KK),apmask(3*ND),cw(3*ND)
  integer*1 msgbits(KK)
  integer apsym(KK)
  integer mcq(28),mde(28),mrrr(16),m73(16),mrr73(16)
  integer itone(NN)
  integer indxs1(8*ND)
  integer icos7(0:6),ip(1)
  integer nappasses(0:5)  !Number of decoding passes to use for each QSO state
  integer naptypes(0:5,4) ! (nQSOProgress, decoding pass)  maximum of 4 passes for now
  integer*1, target:: i1hiscall(12)
  complex cd0(3200)
  complex ctwk(32)
  complex csymb(32)
  logical first,newdat,lsubtract,lapon,lapcqonly,nagain
  equivalence (s1,s1sort)
  data icos7/4,2,5,6,1,3,0/
  data mcq/1,1,1,1,1,0,1,0,0,0,0,0,1,0,0,0,0,0,1,1,0,0,0,1,1,0,0,1/
  data mrrr/0,1,1,1,1,1,1,0,1,1,0,0,1,1,1,1/
  data m73/0,1,1,1,1,1,1,0,1,1,0,1,0,0,0,0/
  data mde/1,1,1,1,1,1,1,1,0,1,1,0,0,1,0,0,0,0,0,1,1,1,0,1,0,0,0,1/
  data mrr73/0,0,0,0,0,0,1,0,0,0,0,1,0,1,0,1/
  data first/.true./
  save nappasses,naptypes

  if(first) then
     mcq=2*mcq-1
     mde=2*mde-1
     mrrr=2*mrrr-1
     m73=2*m73-1
     mrr73=2*mrr73-1
     nappasses(0)=2
     nappasses(1)=2
     nappasses(2)=2
     nappasses(3)=4
     nappasses(4)=4
     nappasses(5)=3

! iaptype
!------------------------
!   1        CQ     ???    ???
!   2        MyCall ???    ???
!   3        MyCall DxCall ???
!   4        MyCall DxCall RRR
!   5        MyCall DxCall 73
!   6        MyCall DxCall RR73
!   7        ???    DxCall ???

     naptypes(0,1:4)=(/1,2,0,0/)
     naptypes(1,1:4)=(/2,3,0,0/)
     naptypes(2,1:4)=(/2,3,0,0/)
     naptypes(3,1:4)=(/3,4,5,6/)
     naptypes(4,1:4)=(/3,4,5,6/)
     naptypes(5,1:4)=(/3,1,2,0/)
     first=.false.
  endif

  max_iterations=30
  nharderrors=-1
  fs2=12000.0/NDOWN
  dt2=1.0/fs2
  twopi=8.0*atan(1.0)
  delfbest=0.
  ibest=0

  call timer('ft8_down',0)
  call ft8_downsample(dd0,newdat,f1,cd0)   !Mix f1 to baseband and downsample
  call timer('ft8_down',1)

  i0=nint((xdt+0.5)*fs2)                   !Initial guess for start of signal
  smax=0.0
  do idt=i0-8,i0+8                         !Search over +/- one quarter symbol
     call sync8d(cd0,idt,ctwk,0,sync)
     if(sync.gt.smax) then
        smax=sync
        ibest=idt
     endif
  enddo
  xdt2=ibest*dt2                           !Improved estimate for DT

! Now peak up in frequency
  i0=nint(xdt2*fs2)
  smax=0.0
  do ifr=-5,5                              !Search over +/- 2.5 Hz
    delf=ifr*0.5
    dphi=twopi*delf*dt2
    phi=0.0
    do i=1,32
      ctwk(i)=cmplx(cos(phi),sin(phi))
      phi=mod(phi+dphi,twopi)
    enddo
   call sync8d(cd0,i0,ctwk,1,sync)
    if( sync .gt. smax ) then
      smax=sync
      delfbest=delf
    endif
  enddo
  a=0.0
  a(1)=-delfbest
  call twkfreq1(cd0,NP2,fs2,a,cd0)
  xdt=xdt2
  f1=f1+delfbest                           !Improved estimate of DF

  call sync8d(cd0,i0,ctwk,2,sync)

  j=0
  do k=1,NN
    i1=ibest+(k-1)*32
    csymb=cmplx(0.0,0.0)
    if( i1.ge.1 .and. i1+31 .le. NP2 ) csymb=cd0(i1:i1+31)
    call four2a(csymb,32,1,-1,1)
    s2(0:7,k)=abs(csymb(1:8))/1e3
  enddo

! sync quality check
  is1=0
  is2=0
  is3=0
  do k=1,7
    ip=maxloc(s2(:,k))
    if(icos7(k-1).eq.(ip(1)-1)) is1=is1+1
    ip=maxloc(s2(:,k+36))
    if(icos7(k-1).eq.(ip(1)-1)) is2=is2+1
    ip=maxloc(s2(:,k+72))
    if(icos7(k-1).eq.(ip(1)-1)) is3=is3+1
  enddo
! hard sync sum - max is 21
  nsync=is1+is2+is3
  if(nsync .le. 6) then ! bail out
    call timer('badnsync', 0)
    nbadcrc=1
    call timer('badnsync', 1)
    return
  endif

  j=0
  do k=1,NN
     if(k.le.7) cycle
     if(k.ge.37 .and. k.le.43) cycle
     if(k.gt.72) cycle
     j=j+1
     s1(0:7,j)=s2(0:7,k)
  enddo

  call indexx(s1sort,8*ND,indxs1)
  xmeds1=s1sort(indxs1(nint(0.5*8*ND)))
  s1=s1/xmeds1

  do j=1,ND
     i4=3*j-2
     i2=3*j-1
     i1=3*j
! Max amplitude
     ps=s1(0:7,j)
     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))
     bmeta(i4)=r4
     bmeta(i2)=r2
     bmeta(i1)=r1
     bmetap(i4)=r4
     bmetap(i2)=r2
     bmetap(i1)=r1
! Max log metric
     psl=log(ps)
     r1=max(psl(1),psl(3),psl(5),psl(7))-max(psl(0),psl(2),psl(4),psl(6))
     r2=max(psl(2),psl(3),psl(6),psl(7))-max(psl(0),psl(1),psl(4),psl(5))
     r4=max(psl(4),psl(5),psl(6),psl(7))-max(psl(0),psl(1),psl(2),psl(3))
     bmetb(i4)=r4
     bmetb(i2)=r2
     bmetb(i1)=r1

! Metric for Cauchy noise
!     r1=log(ps(1)**3+ps(3)**3+ps(5)**3+ps(7)**3)- &
!        log(ps(0)**3+ps(2)**3+ps(4)**3+ps(6)**3)
!     r2=log(ps(2)**3+ps(3)**3+ps(6)**3+ps(7)**3)- &
!        log(ps(0)**3+ps(1)**3+ps(4)**3+ps(5)**3)
!     r4=log(ps(4)**3+ps(5)**3+ps(6)**3+ps(7)**3)- &
!        log(ps(0)**3+ps(1)**3+ps(2)**3+ps(3)**3)
! Metric for AWGN, no fading
!     bscale=2.5
!     b0=bessi0(bscale*ps(0))
!     b1=bessi0(bscale*ps(1))
!     b2=bessi0(bscale*ps(2))
!     b3=bessi0(bscale*ps(3))
!     b4=bessi0(bscale*ps(4))
!     b5=bessi0(bscale*ps(5))
!     b6=bessi0(bscale*ps(6))
!     b7=bessi0(bscale*ps(7))
!     r1=log(b1+b3+b5+b7)-log(b0+b2+b4+b6)
!     r2=log(b2+b3+b6+b7)-log(b0+b1+b4+b5)
!     r4=log(b4+b5+b6+b7)-log(b0+b1+b2+b3)

     if(nQSOProgress .eq. 0 .or. nQSOProgress .eq. 5) then
! When bits 88:115 are set as ap bits, bit 115 lives in symbol 39 along
! with no-ap bits 116 and 117. Take care of metrics for bits 116 and 117.
        if(j.eq.39) then  ! take care of bits that live in symbol 39
           if(apsym(28).lt.0) then
              bmetap(i2)=max(ps(2),ps(3))-max(ps(0),ps(1))
              bmetap(i1)=max(ps(1),ps(3))-max(ps(0),ps(2))
           else
              bmetap(i2)=max(ps(6),ps(7))-max(ps(4),ps(5))
              bmetap(i1)=max(ps(5),ps(7))-max(ps(4),ps(6))
           endif
        endif
     endif

! When bits 116:143 are set as ap bits, bit 115 lives in symbol 39 along
! with ap bits 116 and 117. Take care of metric for bit 115.
!        if(j.eq.39) then  ! take care of bit 115
!           iii=2*(apsym(29)+1)/2 + (apsym(30)+1)/2  ! known values of bits 116 & 117
!           if(iii.eq.0) bmetap(i4)=ps(4)-ps(0)
!           if(iii.eq.1) bmetap(i4)=ps(5)-ps(1)
!           if(iii.eq.2) bmetap(i4)=ps(6)-ps(2)
!           if(iii.eq.3) bmetap(i4)=ps(7)-ps(3)
!        endif

! bit 144 lives in symbol 48 and will be 1 if it is set as an ap bit.
! take care of metrics for bits 142 and 143
     if(j.eq.48) then  ! bit 144 is always 1
       bmetap(i4)=max(ps(5),ps(7))-max(ps(1),ps(3))
       bmetap(i2)=max(ps(3),ps(7))-max(ps(1),ps(5))
     endif

! bit 154 lives in symbol 52 and will be 0 if it is set as an ap bit
! take care of metrics for bits 155 and 156
     if(j.eq.52) then  ! bit 154 will be 0 if it is set as an ap bit.
        bmetap(i2)=max(ps(2),ps(3))-max(ps(0),ps(1))
        bmetap(i1)=max(ps(1),ps(3))-max(ps(0),ps(2))
     endif

  enddo

  call normalizebmet(bmeta,3*ND)
  call normalizebmet(bmetb,3*ND)
  call normalizebmet(bmetap,3*ND)

  scalefac=2.83
  llr0=scalefac*bmeta
  llr1=scalefac*bmetb
  llra=scalefac*bmetap  ! llr's for use with ap
  apmag=scalefac*(maxval(abs(bmetap))*1.01)

! pass #
!------------------------------
!   1        regular decoding
!   2        erase 24
!   3        erase 48
!   4        ap pass 1
!   5        ap pass 2
!   6        ap pass 3
!   7        ap pass 4, etc.

  if(lapon) then
     if(.not.lapcqonly) then
        npasses=4+nappasses(nQSOProgress)
     else
        npasses=5
     endif
  else
     npasses=4
  endif

  do ipass=1,npasses

     llr=llr0
     if(ipass.eq.2) llr=llr1
     if(ipass.eq.3) llr(1:24)=0.
     if(ipass.eq.4) llr(1:48)=0.
     if(ipass.le.4) then
        apmask=0
        llrap=llr
        iaptype=0
     endif

     if(ipass .gt. 4) then
        if(.not.lapcqonly) then
           iaptype=naptypes(nQSOProgress,ipass-4)
        else
           iaptype=1
        endif
        if(iaptype.ge.3 .and. (abs(f1-nfqso).gt.napwid .and. abs(f1-nftx).gt.napwid) ) cycle
        if(iaptype.eq.1 .or. iaptype.eq.2 ) then ! AP,???,???
           apmask=0
           apmask(88:115)=1    ! first 28 bits are AP
           apmask(144)=1       ! not free text
           llrap=llr
           if(iaptype.eq.1) llrap(88:115)=apmag*mcq
           if(iaptype.eq.2) llrap(88:115)=apmag*apsym(1:28)
           llrap(116:117)=llra(116:117)
           llrap(142:143)=llra(142:143)
           llrap(144)=-apmag
        endif
        if(iaptype.eq.3) then   ! mycall, dxcall, ???
           apmask=0
           apmask(88:115)=1   ! mycall
           apmask(116:143)=1  ! hiscall
           apmask(144)=1      ! not free text
           llrap=llr
           llrap(88:143)=apmag*apsym(1:56)
           llrap(144)=-apmag
        endif
        if(iaptype.eq.4 .or. iaptype.eq.5 .or. iaptype.eq.6) then
           apmask=0
           apmask(88:115)=1   ! mycall
           apmask(116:143)=1  ! hiscall
           apmask(144:159)=1  ! RRR or 73 or RR73
           llrap=llr
           llrap(88:143)=apmag*apsym(1:56)
           if(iaptype.eq.4) llrap(144:159)=apmag*mrrr
           if(iaptype.eq.5) llrap(144:159)=apmag*m73
           if(iaptype.eq.6) llrap(144:159)=apmag*mrr73
        endif
        if(iaptype.eq.7) then   ! ???, dxcall, ???
           apmask=0
           apmask(116:143)=1  ! hiscall
           apmask(144)=1      ! not free text
           llrap=llr
           llrap(115)=llra(115)
           llrap(116:143)=apmag*apsym(29:56)
           llrap(144)=-apmag
        endif
     endif

     cw=0
     call timer('bpd174  ',0)
     call bpdecode174(llrap,apmask,max_iterations,decoded,cw,nharderrors,  &
          niterations)
     call timer('bpd174  ',1)
     dmin=0.0
     if(ndepth.eq.3 .and. nharderrors.lt.0) then
        ndeep=3
        if(abs(nfqso-f1).le.napwid .or. abs(nftx-f1).le.napwid) then
          if((ipass.eq.3 .or. ipass.eq.4) .and. .not.nagain) then
            ndeep=3
          else
            ndeep=4
          endif
        endif
        if(nagain) ndeep=5
        call timer('osd174  ',0)
        call osd174(llrap,apmask,ndeep,decoded,cw,nharderrors,dmin)
        call timer('osd174  ',1)
     endif
     nbadcrc=1
     message='                      '
     xsnr=-99.0
     if(count(cw.eq.0).eq.174) cycle           !Reject the all-zero codeword
     if(nharderrors.ge.0 .and. nharderrors+dmin.lt.60.0 .and. &
        .not.(sync.lt.2.0 .and. nharderrors.gt.35)      .and. &
        .not.(ipass.gt.2 .and. nharderrors.gt.39)       .and. &
        .not.(ipass.eq.4 .and. nharderrors.gt.30)             &
       ) then
        call chkcrc12a(decoded,nbadcrc)
     else
        nharderrors=-1
        cycle
     endif
     i3bit=4*decoded(73) + 2*decoded(74) + decoded(75)
     if(nbadcrc.eq.0) then
        decoded0=decoded
        call extractmessage174(decoded,origmsg,ncrcflag)
        decoded=decoded0

        message(1:12)=origmsg(1:12)
        call genft8(message,mygrid6,bcontest,i3bit,msgsent,msgbits,itone)
        if(lsubtract) call subtractft8(dd0,itone,f1,xdt2)
        xsig=0.0
        xnoi=0.0
        do i=1,79
           xsig=xsig+s2(itone(i),i)**2
           ios=mod(itone(i)+4,7)
           xnoi=xnoi+s2(ios,i)**2
        enddo
        xsnr=0.001
        if(xnoi.gt.0 .and. xnoi.lt.xsig) xsnr=xsig/xnoi-1.0
        xsnr=10.0*log10(xsnr)-27.0
        xsnr2=db(xsig/xbase - 1.0) - 32.0
        if(.not.nagain) xsnr=xsnr2
        if(xsnr .lt. -24.0) xsnr=-24.0

        msg37=origmsg//'               '

        msg37(22:22) = char(48 + i3bit)

        return
     endif
  enddo

  return
end subroutine ft8b

subroutine normalizebmet(bmet,n)
  real bmet(n)

  bmetav=sum(bmet)/real(n)
  bmet2av=sum(bmet*bmet)/real(n)
  var=bmet2av-bmetav*bmetav
  if( var .gt. 0.0 ) then
     bmetsig=sqrt(var)
  else
     bmetsig=sqrt(bmet2av)
  endif
  bmet=bmet/bmetsig
  return
end subroutine normalizebmet

function bessi0(x)
! From Numerical Recipes
   real bessi0,x
   double precision p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y
   save p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9
   data p1,p2,p3,p4,p5,p6,p7/1.0d0,3.5156229d0,3.0899424d0,1.2067492d0, &
      0.2659732d0,0.360768d-1,0.45813d-2/
   data q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,0.1328592d-1,           &
      0.225319d-2,-0.157565d-2,0.916281d-2,-0.2057706d-1,               &
      0.2635537d-1,-0.1647633d-1,0.392377d-2/

   if (abs(x).lt.3.75) then
      y=(x/3.75)**2
      bessi0=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))
   else
      ax=abs(x)
      y=3.75/ax
      bessi0=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4         &
           +y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))))
   endif
   return
end function bessi0


code:
program ft8code

! Provides examples of message packing, LDPC(144,87) encoding, bit and
! symbol ordering, and other details of the FT8 protocol.

  use packjt
  use crc
  include 'ft8_params.f90'               !Set various constants
  include 'ft8_testmsg.f90'
  parameter (NWAVE=NN*NSPS)

  character*40 msg,msgchk
  character*37 msg37
  character*6 c1,c2
  character*9 comment
  character*22 msgsent,message
  character*6 mygrid6
  character bad*1,msgtype*10
  character*87 cbits
  logical bcontest
  integer itone(NN)
  integer dgen(12)
  integer*1 msgbits(KK),decoded(KK),decoded0(KK)
  data mygrid6/'EM48  '/

! Get command-line argument(s)
  nargs=iargc()
  if(nargs.ne.1 .and. nargs.ne.3) then
     print*
     print*,'Program ft8code:  Provides examples of message packing, ',       &
          'LDPC(174,87) encoding,'
     print*,'bit and symbol ordering, and other details of the FT8 protocol.'
     print*
     print*,'Usage: ft8code [-c grid] "message"  # Results for specified message'
     print*,'       ft8code -t                   # Examples of all message types'
     go to 999
  endif

  bcontest=.false.
  call getarg(1,msg)                    !Message to be transmitted
  if(len(trim(msg)).eq.2 .and. msg(1:2).eq.'-t') then
     testmsg(NTEST+1)='KA1ABC RR73; WB9XYZ <KH1/KH7Z> -11'
     nmsg=NTEST+1
  else if(len(trim(msg)).eq.2 .and. msg(1:2).eq.'-c') then
     bcontest=.true.
     call getarg(2,mygrid6)
     call getarg(3,msg)
     msgchk=msg
     nmsg=1
  else
     msgchk=msg
     call fmtmsg(msgchk,iz)          !To upper case; collapse multiple blanks
     nmsg=1
  endif

  write(*,1010)
1010 format("    Message                Decoded              Err? Type"/76("-"))

  do imsg=1,nmsg
     if(nmsg.gt.1) msg=testmsg(imsg)
     call fmtmsg(msg,iz)               !To upper case, collapse multiple blanks
     msgchk=msg

! Generate msgsent, msgbits, and itone
     call packmsg(msg(1:22),dgen,itype,bcontest)
     msgtype=""
     if(itype.eq.1) msgtype="Std Msg"
     if(itype.eq.2) msgtype="Type 1 pfx"
     if(itype.eq.3) msgtype="Type 1 sfx"
     if(itype.eq.4) msgtype="Type 2 pfx"
     if(itype.eq.5) msgtype="Type 2 sfx"
     if(itype.eq.6) msgtype="Free text"
     i3bit=0
     call genft8(msg(1:22),mygrid6,bcontest,i3bit,msgsent,msgbits,itone)

     decoded=msgbits
     i3bit=4*decoded(73) + 2*decoded(74) + decoded(75)
     iFreeText=decoded(57)
     decoded0=decoded
     if(i3bit.eq.1) decoded(57:)=0
     call extractmessage174(decoded,message,ncrcflag)
     decoded=decoded0

    if(bcontest) call fix_contest_msg(mygrid6,message)
    bad=" "
    comment='         '
    if(itype.ne.6 .and. message.ne.msgchk) bad="*"
    if(itype.eq.6 .and. message(1:13).ne.msgchk(1:13)) bad="*"
    if(itype.eq.6 .and. len(trim(msgchk)).gt.13) comment='truncated'
    write(*,1020) imsg,msgchk,message,bad,i3bit,itype,msgtype,comment
1020    format(i2,'.',1x,a22,1x,a22,1x,a1,2i2,1x,a10,1x,a9)

  enddo

  if(nmsg.eq.1) then
     write(*,1030) msgbits(1:56)
1030 format(/'Call1: ',28i1,'    Call2: ',28i1)
     write(*,1032) msgbits(57:72),msgbits(73:75),msgbits(76:87)
1032 format('Grid:  ',16i1,'   3Bit: ',3i1,'    CRC12: ',12i1)
     write(*,1034) itone
1034 format(/'Channel symbols:'/79i1)
  endif

999 end program ft8code


code:
program ft8d

! Decode FT8 data read from *.wav files.

  include 'ft8_params.f90'
  character*12 arg
  character infile*80,datetime*13,message*22
  real s(NH1,NHSYM)
  real candidate(3,100)
  integer ihdr(11)
  integer*2 iwave(NMAX)                 !Generated full-length waveform
  real dd(NMAX)

  nargs=iargc()
  if(nargs.lt.3) then
     print*,'Usage:   ft8d MaxIt Norder file1 [file2 ...]'
     print*,'Example  ft8d   40     2   *.wav'
     go to 999
  endif
  call getarg(1,arg)
  read(arg,*) max_iterations
  call getarg(2,arg)
  read(arg,*) norder
  nfiles=nargs-2

  twopi=8.0*atan(1.0)
  fs=12000.0                             !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 (baud)
  txt=NZ*dt                              !Transmission length (s)
  nfa=100.0
  nfb=3000.0
  nfqso=1500.0

  do ifile=1,nfiles
     call getarg(ifile+2,infile)
     open(10,file=infile,status='old',access='stream')
     read(10,end=999) ihdr,iwave
     close(10)
     j2=index(infile,'.wav')
     read(infile(j2-6:j2-1),*) nutc
     datetime=infile(j2-13:j2-1)
     call sync8(iwave,nfa,nfb,nfqso,s,candidate,ncand)
     syncmin=2.0
     dd=iwave
     do icand=1,ncand
       sync=candidate(3,icand)
       if( sync.lt.syncmin) cycle
       f1=candidate(1,icand)
       xdt=candidate(2,icand)
       nsnr=min(99,nint(10.0*log10(sync)-25.5))
       call ft8b(dd,nfqso,f1,xdt,nharderrors,dmin,nbadcrc,message,xsnr)
       nsnr=xsnr
       xdt=xdt-0.6
       write(*,1110) datetime,0,nsnr,xdt,f1,message,nharderrors,dmin
1110   format(a13,2i4,f6.2,f7.1,' ~ ',a22,i6,f7.1)
     enddo
  enddo   ! ifile loop

999 end program ft8d


code:
program ft8sim

! Generate simulated data for a 15-second HF/6m mode using 8-FSK.
! Output is saved to a *.wav file.

  use wavhdr
  include 'ft8_params.f90'               !Set various constants
  parameter (NWAVE=NN*NSPS)
  type(hdr) h                            !Header for .wav file
  character arg*12,fname*17
  character msg40*40,msg*22,msgsent*22,msg0*22
  character*6 mygrid6
  logical bcontest
  complex c0(0:NMAX-1)
  complex c(0:NMAX-1)
  real wave(NMAX)
  integer itone(NN)
  integer*1 msgbits(KK)
  integer*2 iwave(NMAX)                  !Generated full-length waveform
  data mygrid6/'EM48  '/

! Get command-line argument(s)
  nargs=iargc()
  if(nargs.ne.8) then
     print*,'Usage:    ft8sim "message"         nsig|f0  DT fdop del width nfiles snr'
     print*,'Examples: ft8sim "K1ABC W9XYZ EN37" 1500.0 0.0  0.1 1.0   0     10   -18'
     print*,'          ft8sim "K1ABC W9XYZ EN37"   10   0.0  0.1 1.0  25     10   -18'
     print*,'          ft8sim "K1ABC W9XYZ EN37"   25   0.0  0.1 1.0  25     10   -18'
     print*,'          ft8sim "K1ABC RR73; W9XYZ <KH1/KH7Z> -11" 300 0 0 0 25 1 -10'
     print*,'Make nfiles negative to invoke 72-bit contest mode.'
     go to 999
  endif
  call getarg(1,msg40)                   !Message to be transmitted
  call getarg(2,arg)
  read(arg,*) f0                         !Frequency (only used for single-signal)
  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,*) width                      !Filter transition width (Hz)
  call getarg(7,arg)
  read(arg,*) nfiles                     !Number of files
  call getarg(8,arg)
  read(arg,*) snrdb                      !SNR_2500
  nsig=1
  if(f0.lt.100.0) then
     nsig=f0
     f0=1500
  endif

  bcontest=nfiles.lt.0
  nfiles=abs(nfiles)
  twopi=8.0*atan(1.0)
  fs=12000.0                             !Sample rate (Hz)
  dt=1.0/fs                              !Sample interval (s)
  tt=NSPS*dt                             !Duration of symbols (s)
  baud=1.0/tt                            !Keying rate (baud)
  bw=8*baud                              !Occupied bandwidth (Hz)
  txt=NZ*dt                              !Transmission length (s)
  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
  txt=NN*NSPS/12000.0

! Source-encode, then get itone()
  if(index(msg40,';').le.0) then
     i3bit=0
     msg=msg40(1:22)
     call genft8(msg,mygrid6,bcontest,i3bit,msgsent,msgbits,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)
  else
     call foxgen_wrap(msg40,msgbits,itone)
     write(*,1001) f0,xdt,txt,snrdb,bw,msg40
1001 format('f0:',f9.3,'   DT:',f6.2,'   TxT:',f6.1,'   SNR:',f6.1,    &
          '  BW:',f4.1,2x,a40)
  endif

  write(*,1030) msgbits(1:56)
1030 format(/'Call1: ',28i1,'    Call2: ',28i1)
  write(*,1032) msgbits(57:72),msgbits(73:75),msgbits(76:87)
1032 format('Grid:  ',16i1,'   3Bit: ',3i1,'    CRC12: ',12i1)
  write(*,1034) itone
1034 format(/'Channel symbols:'/79i1/)

  msg0=msg
  do ifile=1,nfiles
     c=0.
     do isig=1,nsig
        c0=0.
        if(nsig.eq.2) then
           if(index(msg,'R-').gt.0) f0=500
           i1=index(msg,' ')
           msg(i1+4:i1+4)=char(ichar('A')+isig-1)
           if(isig.eq.2) then
              f0=f0+100
           endif
           call genft8(msg,mygrid6,bcontest,i3bit,msgsent,msgbits,itone)
        endif
        if(nsig.eq.25) then
           f0=(isig+2)*100.0
        else if(nsig.eq.50) then
           msg=msg0
           f0=1000.0 + (isig-1)*60.0
           i1=index(msg,' ')
           i2=index(msg(i1+1:),' ') + i1
           msg(i1+2:i1+2)=char(ichar('0')+mod(isig-1,10))
           msg(i1+3:i1+3)=char(ichar('A')+mod(isig-1,26))
           msg(i1+4:i1+4)=char(ichar('A')+mod(isig-1,26))
           msg(i1+5:i1+5)=char(ichar('A')+mod(isig-1,26))
           write(msg(i2+3:i2+4),'(i2.2)') isig-1
           if(ifile.ge.2 .and. isig.eq.ifile-1) then
              write(msg(i2+1:i2+4),1002) -isig
1002          format('R',i3.2)
              f0=600.0 + mod(isig-1,5)*60.0
           endif
           call genft8(msg,mygrid6,bcontest,i3bit,msgsent,msgbits,itone)
        endif
        k=-1 + nint((xdt+0.5+0.01*gran())/dt)
!        k=-1 + nint((xdt+0.5)/dt)
        ia=k+1
        phi=0.0
        do j=1,NN                             !Generate complex waveform
           dphi=twopi*(f0+itone(j)*baud)*dt
           do i=1,NSPS
              k=k+1
              phi=mod(phi+dphi,twopi)
              if(k.ge.0 .and. k.lt.NMAX) c0(k)=cmplx(cos(phi),sin(phi))
           enddo
        enddo
        if(fspread.ne.0.0 .or. delay.ne.0.0) call watterson(c0,NMAX,fs,delay,fspread)
        c=c+sig*c0
     enddo
     ib=k
     wave=real(c)
     peak=maxval(abs(wave(ia:ib)))
     rms=sqrt(dot_product(wave(ia:ib),wave(ia:ib))/NWAVE)
     nslots=1
     if(width.gt.0.0) call filt8(f0,nslots,width,wave)

     if(snrdb.lt.90) then
        do i=1,NMAX                   !Add gaussian noise at specified SNR
           xnoise=gran()
!           wave(i)=wave(i) + xnoise
!           if(i.ge.ia .and. i.le.ib) write(30,3001) i,wave(i)/peak
!3001       format(i8,f12.6)
           wave(i)=wave(i) + xnoise
        enddo
     endif

     fac=32767.0
     rms=100.0
     if(snrdb.ge.90.0) iwave(1:NMAX)=nint(fac*wave)
     if(snrdb.lt.90.0) iwave(1:NMAX)=nint(rms*wave)

     h=default_header(12000,NMAX)
     write(fname,1102) ifile
1102 format('000000_',i6.6,'.wav')
     open(10,file=fname,status='unknown',access='stream')
     write(10) h,iwave                !Save to *.wav file
     close(10)
     write(*,1110) ifile,xdt,f0,snrdb,fname
1110 format(i4,f7.2,f8.2,f7.1,2x,a17)
  enddo

999 end program ft8sim


code:
subroutine genft8(msg,mygrid,bcontest,i3bit,msgsent,msgbits,itone)

! Encode an FT8 message, producing array itone().

  use crc
  use packjt
  include 'ft8_params.f90'
  character*68 alphabet
  character*22 msg,msgsent
  character*6 mygrid
  character*87 cbits
  logical bcontest,checksumok
  integer*4 i4Msg6BitWords(12)                !72-bit message as 6-bit words
  integer*1 msgbits(KK),codeword(3*ND)
  integer*1, target:: i1Msg8BitBytes(11)
  integer itone(NN)
  integer icos7(0:6)
  data icos7/4,2,5,6,1,3,0/                   !Costas 7x7 tone pattern

  alphabet='0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz-+/?.'

  itype=6
  do i=1,12
    v=index(alphabet, msg(i:i))
    if(v.eq.0) exit
    i4Msg6BitWords(i)=v - 1
  enddo
  msgsent='                      '
  msgsent(1:12)=msg(1:12)

  ! call packmsg(msg,i4Msg6BitWords,itype,bcontest) !Pack into 12 6-bit bytes
  ! call unpackmsg(i4Msg6BitWords,msgsent,bcontest,mygrid) !Unpack to get msgsent

  write(cbits,1000) i4Msg6BitWords,32*i3bit
1000 format(12b6.6,b8.8)
  read(cbits,1001) i1Msg8BitBytes(1:10)
1001 format(10b8)
  i1Msg8BitBytes(10)=iand(i1Msg8BitBytes(10),128+64+32)
  i1Msg8BitBytes(11)=0
  icrc12=crc12(c_loc(i1Msg8BitBytes),11)
!  icrc12=xor(icrc12, 41) ! TODO: jsherer - could change the crc here

! For reference, here's how to check the CRC
!  i1Msg8BitBytes(10)=icrc12/256
!  i1Msg8BitBytes(11)=iand (icrc12,255)
!  checksumok = crc12_check(c_loc (i1Msg8BitBytes), 11)
!  if( checksumok ) write(*,*) 'Good checksum'

  write(cbits,1003) i4Msg6BitWords,i3bit,icrc12
1003 format(12b6.6,b3.3,b12.12)
  read(cbits,1004) msgbits
1004 format(87i1)

  call encode174(msgbits,codeword)      !Encode the test message

! Message structure: S7 D29 S7 D29 S7
  itone(1:7)=icos7
  itone(36+1:36+7)=icos7
  itone(NN-6:NN)=icos7
  k=7
  do j=1,ND
     i=3*j -2
     k=k+1
     if(j.eq.30) k=k+7
     itone(k)=codeword(i)*4 + codeword(i+1)*2 + codeword(i+2)
  enddo

  return
end subroutine genft8


code:
subroutine genft8refsig(itone,cref,f0)
  complex cref(79*1920)
  integer itone(79)
  real*8 twopi,phi,dphi,dt,xnsps
  data twopi/0.d0/
  save twopi
  if( twopi .lt. 0.1 ) twopi=8.d0*atan(1.d0)

  xnsps=1920.d0
  dt=1.d0/12000.d0
  phi=0.d0
  k=1
  do i=1,79
    dphi=twopi*(f0*dt+itone(i)/xnsps)
    do is=1,1920
      cref(k)=cmplx(cos(phi),sin(phi))
      phi=mod(phi+dphi,twopi)
      k=k+1
    enddo
  enddo
  return
end subroutine genft8refsig


code:
real function h1(x)

!  sigma=1.0/sqrt(2.0)
  sigma=1.0
  xlim=sigma/sqrt(6.0)
  ax=abs(x)
  sgnx=1.0
  if(x.lt.0) sgnx=-1.0
  if(ax.le.xlim) then
     h1=x
  else
     z=exp(1.0/6.0 - (ax/sigma)**2)
     h1=sgnx*sqrt(6.0)*sigma*(2.0/3.0 - 0.5*z)
  endif

  return
end function h1


code:
program ldpcsim174
! End to end test of the (174,75)/crc12 encoder and decoder.
use crc
use packjt

character*22 msg,msgsent,msgreceived
character*8 arg
character*6 grid
integer*1, allocatable ::  codeword(:), decoded(:), message(:)
integer*1, target:: i1Msg8BitBytes(11)
integer*1 msgbits(87)
integer*1 apmask(174), cw(174)
integer*2 checksum
integer*4 i4Msg6BitWords(13)
integer colorder(174)
integer nerrtot(174),nerrdec(174),nmpcbad(87)
logical checksumok,fsk,bpsk
real*8, allocatable ::  rxdata(:)
real, allocatable :: llr(:)

data colorder/            &
   0,  1,  2,  3, 30,  4,  5,  6,  7,  8,  9, 10, 11, 32, 12, 40, 13, 14, 15, 16,&
  17, 18, 37, 45, 29, 19, 20, 21, 41, 22, 42, 31, 33, 34, 44, 35, 47, 51, 50, 43,&
  36, 52, 63, 46, 25, 55, 27, 24, 23, 53, 39, 49, 59, 38, 48, 61, 60, 57, 28, 62,&
  56, 58, 65, 66, 26, 70, 64, 69, 68, 67, 74, 71, 54, 76, 72, 75, 78, 77, 80, 79,&
  73, 83, 84, 81, 82, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99,&
 100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,&
 120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,&
 140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,&
 160,161,162,163,164,165,166,167,168,169,170,171,172,173/

nerrtot=0
nerrdec=0
nmpcbad=0  ! Used to collect the number of errors in the message+crc part of the codeword

nargs=iargc()
if(nargs.ne.4) then
   print*,'Usage: ldpcsim  niter  ndepth  #trials   s '
   print*,'eg:    ldpcsim    10     2      1000    0.84'
   print*,'belief propagation iterations: niter, ordered-statistics depth: ndepth'
   print*,'If s is negative, then value is ignored and sigma is calculated from SNR.'
   return
endif
call getarg(1,arg)
read(arg,*) max_iterations
call getarg(2,arg)
read(arg,*) ndepth
call getarg(3,arg)
read(arg,*) ntrials
call getarg(4,arg)
read(arg,*) s

fsk=.false.
bpsk=.true.

! don't count crc bits as data bits
N=174
K=87
! scale Eb/No for a (174,87) code
rate=real(K)/real(N)

write(*,*) "rate: ",rate
write(*,*) "niter= ",max_iterations," s= ",s

allocate ( codeword(N), decoded(K), message(K) )
allocate ( rxdata(N), llr(N) )

  msg="K1JT K9AN EN50"
!  msg="G4WJS K9AN EN50"
  call packmsg(msg,i4Msg6BitWords,itype,.false.) !Pack into 12 6-bit bytes
  call unpackmsg(i4Msg6BitWords,msgsent,.false.,grid) !Unpack to get msgsent
  write(*,*) "message sent ",msgsent

  i4=0
  ik=0
  im=0
  do i=1,12
    nn=i4Msg6BitWords(i)
    do j=1, 6
      ik=ik+1
      i4=i4+i4+iand(1,ishft(nn,j-6))
      i4=iand(i4,255)
      if(ik.eq.8) then
        im=im+1
!           if(i4.gt.127) i4=i4-256
        i1Msg8BitBytes(im)=i4
        ik=0
      endif
    enddo
  enddo

  i1Msg8BitBytes(10:11)=0
  checksum = crc12 (c_loc (i1Msg8BitBytes), 11)
!  checksum = xor(checksum, 41) ! TODO: jsherer - could change the crc here
! For reference, the next 3 lines show how to check the CRC
  i1Msg8BitBytes(10)=checksum/256
  i1Msg8BitBytes(11)=iand (checksum,255)
  checksumok = crc12_check(c_loc (i1Msg8BitBytes), 11)
  if( checksumok ) write(*,*) 'Good checksum'

! K=87, For now:
! msgbits(1:72) JT message bits
! msgbits(73:75) 3 free message bits (set to 0)
! msgbits(76:87) CRC12
  mbit=0
  do i=1, 9
    i1=i1Msg8BitBytes(i)
    do ibit=1,8
      mbit=mbit+1
      msgbits(mbit)=iand(1,ishft(i1,ibit-8))
    enddo
  enddo
  msgbits(73:75)=0  ! the three extra message bits go here
  i1=i1Msg8BitBytes(10) ! First 4 bits of crc12 are LSB of this byte
  do ibit=1,4
    msgbits(75+ibit)=iand(1,ishft(i1,ibit-4))
  enddo
  i1=i1Msg8BitBytes(11) ! Now shift in last 8 bits of the CRC
  do ibit=1,8
    msgbits(79+ibit)=iand(1,ishft(i1,ibit-8))
  enddo

  write(*,*) 'message'
  write(*,'(11(8i1,1x))') msgbits

  call encode174(msgbits,codeword)
  call init_random_seed()
!  call sgran()

  write(*,*) 'codeword'
  write(*,'(22(8i1,1x))') codeword

write(*,*) "Es/N0   SNR2500   ngood  nundetected nbadcrc   sigma"
do idb = 20,-10,-1
!do idb = -3,-3,-1
  db=idb/2.0-1.0
  sigma=1/sqrt( 2*(10**(db/10.0)) )
  ngood=0
  nue=0
  nbadcrc=0
  nberr=0
  do itrial=1, ntrials
! Create a realization of a noisy received word
    do i=1,N
      if( bpsk ) then
        rxdata(i) = 2.0*codeword(i)-1.0 + sigma*gran()
      elseif( fsk ) then
        if( codeword(i) .eq. 1 ) then
          r1=(1.0 + sigma*gran())**2 + (sigma*gran())**2
          r2=(sigma*gran())**2 + (sigma*gran())**2
        elseif( codeword(i) .eq. 0 ) then
          r2=(1.0 + sigma*gran())**2 + (sigma*gran())**2
          r1=(sigma*gran())**2 + (sigma*gran())**2
        endif
!        rxdata(i)=0.35*(sqrt(r1)-sqrt(r2))
!        rxdata(i)=0.35*(exp(r1)-exp(r2))
        rxdata(i)=0.12*(log(r1)-log(r2))
      endif
    enddo
    nerr=0
    do i=1,N
      if( rxdata(i)*(2*codeword(i)-1.0) .lt. 0 ) nerr=nerr+1
    enddo
    if(nerr.ge.1) nerrtot(nerr)=nerrtot(nerr)+1
    nberr=nberr+nerr

! Correct signal normalization is important for this decoder.
    rxav=sum(rxdata)/N
    rx2av=sum(rxdata*rxdata)/N
    rxsig=sqrt(rx2av-rxav*rxav)
    rxdata=rxdata/rxsig
! To match the metric to the channel, s should be set to the noise standard deviation.
! For now, set s to the value that optimizes decode probability near threshold.
! The s parameter can be tuned to trade a few tenth's dB of threshold for an order of
! magnitude in UER
    if( s .lt. 0 ) then
      ss=sigma
    else
      ss=s
    endif

    llr=2.0*rxdata/(ss*ss)
    nap=0 ! number of AP bits
    llr(colorder(174-87+1:174-87+nap)+1)=5*(2.0*msgbits(1:nap)-1.0)
    apmask=0
    apmask(colorder(174-87+1:174-87+nap)+1)=1

! max_iterations is max number of belief propagation iterations
    call bpdecode174(llr, apmask, max_iterations, decoded, cw, nharderrors,niterations)
    if( ndepth .ge. 0 .and. nharderrors .lt. 0 ) call osd174(llr, apmask, ndepth, decoded, cw,  nharderrors, dmin)
! If the decoder finds a valid codeword, nharderrors will be .ge. 0.
    if( nharderrors .ge. 0 ) then
      call extractmessage174(decoded,msgreceived,ncrcflag)
      if( ncrcflag .ne. 1 ) then
        nbadcrc=nbadcrc+1
      endif

      nueflag=0
      nerrmpc=0
      do i=1,K   ! find number of errors in message+crc part of codeword
        if( msgbits(i) .ne. decoded(i) ) then
          nueflag=1
          nerrmpc=nerrmpc+1
        endif
      enddo
      if(nerrmpc.ge.1) nmpcbad(nerrmpc)=nmpcbad(nerrmpc)+1
      if( ncrcflag .eq. 1 ) then
        if( nueflag .eq. 0 ) then
          ngood=ngood+1
          if(nerr.ge.1) nerrdec(nerr)=nerrdec(nerr)+1
        else if( nueflag .eq. 1 ) then
          nue=nue+1;
        endif
      endif
    endif
  enddo
  baud=12000/1920
  snr2500=db+10.0*log10((baud/2500.0))
  pberr=real(nberr)/(real(ntrials*N))
  write(*,"(f4.1,4x,f5.1,1x,i8,1x,i8,1x,i8,8x,f5.2,8x,e10.3)") db,snr2500,ngood,nue,nbadcrc,ss,pberr

enddo

open(unit=23,file='nerrhisto.dat',status='unknown')
do i=1,174
  write(23,'(i4,2x,i10,i10,f10.2)') i,nerrdec(i),nerrtot(i),real(nerrdec(i))/real(nerrtot(i)+1e-10)
enddo
close(23)
open(unit=25,file='nmpcbad.dat',status='unknown')
do i=1,87
  write(25,'(i4,2x,i10)') i,nmpcbad(i)
enddo
close(25)

end program ldpcsim174


code:
subroutine osd174(llr,apmask,ndeep,decoded,cw,nhardmin,dmin)
!
! An ordered-statistics decoder for the (174,87) code.
!
include "ldpc_174_87_params.f90"

integer*1 apmask(N),apmaskr(N)
integer*1 gen(K,N)
integer*1 genmrb(K,N),g2(N,K)
integer*1 temp(K),m0(K),me(K),mi(K),misub(K),e2sub(N-K),e2(N-K),ui(N-K)
integer*1 r2pat(N-K)
integer indices(N),nxor(N)
integer*1 cw(N),ce(N),c0(N),hdec(N)
integer*1 decoded(K)
integer indx(N)
real llr(N),rx(N),absrx(N)
logical first,reset
data first/.true./
save first,gen

if( first ) then ! fill the generator matrix
  gen=0
  do i=1,M
    do j=1,22
      read(g(i)(j:j),"(Z1)") istr
        do jj=1, 4
          irow=(j-1)*4+jj
          if( btest(istr,4-jj) ) gen(irow,i)=1
        enddo
    enddo
  enddo
  do irow=1,K
    gen(irow,M+irow)=1
  enddo
first=.false.
endif

! Re-order received vector to place systematic msg bits at the end.
rx=llr(colorder+1)
apmaskr=apmask(colorder+1)

! Hard decisions on the received word.
hdec=0
where(rx .ge. 0) hdec=1

! Use magnitude of received symbols as a measure of reliability.
absrx=abs(rx)
call indexx(absrx,N,indx)

! Re-order the columns of the generator matrix in order of decreasing reliability.
do i=1,N
  genmrb(1:K,i)=gen(1:K,indx(N+1-i))
  indices(i)=indx(N+1-i)
enddo

! Do gaussian elimination to create a generator matrix with the most reliable
! received bits in positions 1:K in order of decreasing reliability (more or less).
do id=1,K ! diagonal element indices
  do icol=id,K+20  ! The 20 is ad hoc - beware
    iflag=0
    if( genmrb(id,icol) .eq. 1 ) then
      iflag=1
      if( icol .ne. id ) then ! reorder column
        temp(1:K)=genmrb(1:K,id)
        genmrb(1:K,id)=genmrb(1:K,icol)
        genmrb(1:K,icol)=temp(1:K)
        itmp=indices(id)
        indices(id)=indices(icol)
        indices(icol)=itmp
      endif
      do ii=1,K
        if( ii .ne. id .and. genmrb(ii,id) .eq. 1 ) then
          genmrb(ii,1:N)=ieor(genmrb(ii,1:N),genmrb(id,1:N))
        endif
      enddo
      exit
    endif
  enddo
enddo

g2=transpose(genmrb)

! The hard decisions for the K MRB bits define the order 0 message, m0.
! Encode m0 using the modified generator matrix to find the "order 0" codeword.
! Flip various combinations of bits in m0 and re-encode to generate a list of
! codewords. Return the member of the list that has the smallest Euclidean
! distance to the received word.

hdec=hdec(indices)   ! hard decisions from received symbols
m0=hdec(1:K)         ! zero'th order message
absrx=absrx(indices)
rx=rx(indices)
apmaskr=apmaskr(indices)

call mrbencode(m0,c0,g2,N,K)
nxor=ieor(c0,hdec)
nhardmin=sum(nxor)
dmin=sum(nxor*absrx)

cw=c0
ntotal=0
nrejected=0

if(ndeep.eq.0) goto 998  ! norder=0
if(ndeep.gt.5) ndeep=5
if( ndeep.eq. 1) then
   nord=1
   npre1=0
   npre2=0
   nt=40
   ntheta=12
elseif(ndeep.eq.2) then
   nord=1
   npre1=1
   npre2=0
   nt=40
   ntheta=12
elseif(ndeep.eq.3) then
   nord=1
   npre1=1
   npre2=1
   nt=40
   ntheta=12
   ntau=14
elseif(ndeep.eq.4) then
   nord=2
   npre1=1
   npre2=0
   nt=40
   ntheta=12
   ntau=19
elseif(ndeep.eq.5) then
   nord=2
   npre1=1
   npre2=1
   nt=40
   ntheta=12
   ntau=19
endif

do iorder=1,nord
   misub(1:K-iorder)=0
   misub(K-iorder+1:K)=1
   iflag=K-iorder+1
   do while(iflag .ge.0)
      if(iorder.eq.nord .and. npre1.eq.0) then
         iend=iflag
      else
         iend=1
      endif
      do n1=iflag,iend,-1
         mi=misub
         mi(n1)=1
         if(any(iand(apmaskr(1:K),mi).eq.1)) cycle
         ntotal=ntotal+1
         me=ieor(m0,mi)
         if(n1.eq.iflag) then
            call mrbencode(me,ce,g2,N,K)
            e2sub=ieor(ce(K+1:N),hdec(K+1:N))
            e2=e2sub
            nd1Kpt=sum(e2sub(1:nt))+1
            d1=sum(ieor(me(1:K),hdec(1:K))*absrx(1:K))
         else
            e2=ieor(e2sub,g2(K+1:N,n1))
            nd1Kpt=sum(e2(1:nt))+2
         endif
         if(nd1Kpt .le. ntheta) then
            call mrbencode(me,ce,g2,N,K)
            nxor=ieor(ce,hdec)
            if(n1.eq.iflag) then
               dd=d1+sum(e2sub*absrx(K+1:N))
            else
               dd=d1+ieor(ce(n1),hdec(n1))*absrx(n1)+sum(e2*absrx(K+1:N))
            endif
            if( dd .lt. dmin ) then
               dmin=dd
               cw=ce
               nhardmin=sum(nxor)
               nd1Kptbest=nd1Kpt
            endif
         else
            nrejected=nrejected+1
         endif
      enddo
! Get the next test error pattern, iflag will go negative
! when the last pattern with weight iorder has been generated.
      call nextpat(misub,k,iorder,iflag)
   enddo
enddo

if(npre2.eq.1) then
   reset=.true.
   ntotal=0
   do i1=K,1,-1
      do i2=i1-1,1,-1
         ntotal=ntotal+1
         mi(1:ntau)=ieor(g2(K+1:K+ntau,i1),g2(K+1:K+ntau,i2))
         call boxit(reset,mi(1:ntau),ntau,ntotal,i1,i2)
      enddo
   enddo

   ncount2=0
   ntotal2=0
   reset=.true.
! Now run through again and do the second pre-processing rule
   misub(1:K-nord)=0
   misub(K-nord+1:K)=1
   iflag=K-nord+1
   do while(iflag .ge.0)
      me=ieor(m0,misub)
      call mrbencode(me,ce,g2,N,K)
      e2sub=ieor(ce(K+1:N),hdec(K+1:N))
      do i2=0,ntau
         ntotal2=ntotal2+1
         ui=0
         if(i2.gt.0) ui(i2)=1
         r2pat=ieor(e2sub,ui)
778      continue
            call fetchit(reset,r2pat(1:ntau),ntau,in1,in2)
            if(in1.gt.0.and.in2.gt.0) then
               ncount2=ncount2+1
               mi=misub
               mi(in1)=1
               mi(in2)=1
               if(sum(mi).lt.nord+npre1+npre2.or.any(iand(apmaskr(1:K),mi).eq.1)) cycle
               me=ieor(m0,mi)
               call mrbencode(me,ce,g2,N,K)
               nxor=ieor(ce,hdec)
               dd=sum(nxor*absrx)
               if( dd .lt. dmin ) then
                  dmin=dd
                  cw=ce
                  nhardmin=sum(nxor)
               endif
               goto 778
             endif
      enddo
      call nextpat(misub,K,nord,iflag)
   enddo
endif

998 continue
! Re-order the codeword to place message bits at the end.
cw(indices)=cw
hdec(indices)=hdec
decoded=cw(M+1:N)
cw(colorder+1)=cw ! put the codeword back into received-word order
return
end subroutine osd174

subroutine mrbencode(me,codeword,g2,N,K)
integer*1 me(K),codeword(N),g2(N,K)
! fast encoding for low-weight test patterns
  codeword=0
  do i=1,K
    if( me(i) .eq. 1 ) then
      codeword=ieor(codeword,g2(1:N,i))
    endif
  enddo
return
end subroutine mrbencode

subroutine nextpat(mi,k,iorder,iflag)
  integer*1 mi(k),ms(k)
! generate the next test error pattern
  ind=-1
  do i=1,k-1
     if( mi(i).eq.0 .and. mi(i+1).eq.1) ind=i
  enddo
  if( ind .lt. 0 ) then ! no more patterns of this order
    iflag=ind
    return
  endif
  ms=0
  ms(1:ind-1)=mi(1:ind-1)
  ms(ind)=1
  ms(ind+1)=0
  if( ind+1 .lt. k ) then
     nz=iorder-sum(ms)
     ms(k-nz+1:k)=1
  endif
  mi=ms
  do i=1,k  ! iflag will point to the lowest-index 1 in mi
    if(mi(i).eq.1) then
      iflag=i
      exit
    endif
  enddo
  return
end subroutine nextpat

subroutine boxit(reset,e2,ntau,npindex,i1,i2)
  integer*1 e2(1:ntau)
  integer   indexes(4000,2),fp(0:525000),np(4000)
  logical reset
  common/boxes/indexes,fp,np

  if(reset) then
    patterns=-1
    fp=-1
    np=-1
    sc=-1
    indexes=-1
    reset=.false.
  endif

  indexes(npindex,1)=i1
  indexes(npindex,2)=i2
  ipat=0
  do i=1,ntau
    if(e2(i).eq.1) then
      ipat=ipat+ishft(1,ntau-i)
    endif
  enddo

  ip=fp(ipat)   ! see what's currently stored in fp(ipat)
  if(ip.eq.-1) then
    fp(ipat)=npindex
  else
     do while (np(ip).ne.-1)
      ip=np(ip)
     enddo
     np(ip)=npindex
  endif
  return
end subroutine boxit

subroutine fetchit(reset,e2,ntau,i1,i2)
  integer   indexes(4000,2),fp(0:525000),np(4000)
  integer   lastpat
  integer*1 e2(ntau)
  logical reset
  common/boxes/indexes,fp,np
  save lastpat,inext

  if(reset) then
    lastpat=-1
    reset=.false.
  endif

  ipat=0
  do i=1,ntau
    if(e2(i).eq.1) then
      ipat=ipat+ishft(1,ntau-i)
    endif
  enddo
  index=fp(ipat)

  if(lastpat.ne.ipat .and. index.gt.0) then ! return first set of indices
     i1=indexes(index,1)
     i2=indexes(index,2)
     inext=np(index)
  elseif(lastpat.eq.ipat .and. inext.gt.0) then
     i1=indexes(inext,1)
     i2=indexes(inext,2)
     inext=np(inext)
  else
     i1=-1
     i2=-1
     inext=-1
  endif
  lastpat=ipat
  return
end subroutine fetchit


code:
subroutine subtractft8(dd,itone,f0,dt)

! Subtract an ft8 signal
!
! Measured signal  : dd(t)    = a(t)cos(2*pi*f0*t+theta(t))
! Reference signal : cref(t)  = exp( j*(2*pi*f0*t+phi(t)) )
! Complex amp      : cfilt(t) = LPF[ dd(t)*CONJG(cref(t)) ]
! Subtract         : dd(t)    = dd(t) - 2*REAL{cref*cfilt}

  use timer_module, only: timer

  parameter (NMAX=15*12000,NFRAME=1920*79)
  parameter (NFFT=NMAX,NFILT=1400)
  real*4  dd(NMAX), window(-NFILT/2:NFILT/2)
  complex cref,camp,cfilt,cw
  integer itone(79)
  logical first
  data first/.true./
  common/heap8/cref(NFRAME),camp(NMAX),cfilt(NMAX),cw(NMAX)
  save first

  nstart=dt*12000+1
  call genft8refsig(itone,cref,f0)
  camp=0.
  do i=1,nframe
    id=nstart-1+i
    if(id.ge.1.and.id.le.NMAX) camp(i)=dd(id)*conjg(cref(i))
  enddo

  if(first) then
! Create and normalize the filter
     pi=4.0*atan(1.0)
     fac=1.0/float(nfft)
     sum=0.0
     do j=-NFILT/2,NFILT/2
        window(j)=cos(pi*j/NFILT)**2
        sum=sum+window(j)
     enddo
     cw=0.
     cw(1:NFILT+1)=window/sum
     cw=cshift(cw,NFILT/2+1)
     call four2a(cw,nfft,1,-1,1)
     cw=cw*fac
     first=.false.
  endif

  cfilt=0.0
  cfilt(1:nframe)=camp(1:nframe)
  call four2a(cfilt,nfft,1,-1,1)
  cfilt(1:nfft)=cfilt(1:nfft)*cw(1:nfft)
  call four2a(cfilt,nfft,1,1,1)

! Subtract the reconstructed signal
  do i=1,nframe
     j=nstart+i-1
     if(j.ge.1 .and. j.le.NMAX) dd(j)=dd(j)-2*REAL(cfilt(i)*cref(i))
  enddo

  return
end subroutine subtractft8


code:
subroutine sync8(dd,nfa,nfb,syncmin,nfqso,s,candidate,ncand,sbase)

  include 'ft8_params.f90'
! Search over +/- 2.5s relative to 0.5s TX start time.
  parameter (JZ=62)
  complex cx(0:NH1)
  real s(NH1,NHSYM)
  real savg(NH1)
  real sbase(NH1)
  real x(NFFT1)
  real sync2d(NH1,-JZ:JZ)
  real red(NH1)
  real candidate0(3,200)
  real candidate(3,200)
  real dd(NMAX)
  integer jpeak(NH1)
  integer indx(NH1)
  integer ii(1)
  integer icos7(0:6)
  data icos7/4,2,5,6,1,3,0/                   !Costas 7x7 tone pattern
  equivalence (x,cx)

! Compute symbol spectra, stepping by NSTEP steps.
  savg=0.
  tstep=NSTEP/12000.0
  df=12000.0/NFFT1                            !3.125 Hz
  fac=1.0/300.0
  do j=1,NHSYM
     ia=(j-1)*NSTEP + 1
     ib=ia+NSPS-1
     x(1:NSPS)=fac*dd(ia:ib)
     x(NSPS+1:)=0.
     call four2a(x,NFFT1,1,-1,0)              !r2c FFT
     do i=1,NH1
        s(i,j)=real(cx(i))**2 + aimag(cx(i))**2
     enddo
     savg=savg + s(1:NH1,j)                   !Average spectrum
  enddo
  call baseline(savg,nfa,nfb,sbase)
!  savg=savg/NHSYM
!  do i=1,NH1
!     write(51,3051) i*df,savg(i),db(savg(i))
!3051 format(f10.3,e12.3,f12.3)
!  enddo

  ia=max(1,nint(nfa/df))
  ib=nint(nfb/df)
  nssy=NSPS/NSTEP   ! # steps per symbol
  nfos=NFFT1/NSPS   ! # frequency bin oversampling factor
  jstrt=0.5/tstep

  do i=ia,ib
     do j=-JZ,+JZ
        ta=0.
        tb=0.
        tc=0.
        t0a=0.
        t0b=0.
        t0c=0.
        do n=0,6
           k=j+jstrt+nssy*n
           if(k.ge.1.and.k.le.NHSYM) then
              ta=ta + s(i+nfos*icos7(n),k)
              t0a=t0a + sum(s(i:i+nfos*6:nfos,k))
           endif
           tb=tb + s(i+nfos*icos7(n),k+nssy*36)
           t0b=t0b + sum(s(i:i+nfos*6:nfos,k+nssy*36))
           if(k+nssy*72.le.NHSYM) then
              tc=tc + s(i+nfos*icos7(n),k+nssy*72)
              t0c=t0c + sum(s(i:i+nfos*6:nfos,k+nssy*72))
           endif
        enddo
        t=ta+tb+tc
        t0=t0a+t0b+t0c
        t0=(t0-t)/6.0
        sync_abc=t/t0

        t=tb+tc
        t0=t0b+t0c
        t0=(t0-t)/6.0
        sync_bc=t/t0
        sync2d(i,j)=max(sync_abc,sync_bc)
     enddo
  enddo

  red=0.
  do i=ia,ib
     ii=maxloc(sync2d(i,-JZ:JZ)) - 1 - JZ
     j0=ii(1)
     jpeak(i)=j0
     red(i)=sync2d(i,j0)
!     write(52,3052) i*df,red(i),db(red(i))
!3052 format(3f12.3)
  enddo
  iz=ib-ia+1
  call indexx(red(ia:ib),iz,indx)
  ibase=indx(nint(0.40*iz)) - 1 + ia
  base=red(ibase)
  red=red/base

  candidate0=0.
  k=0
  do i=1,200
     ji=iz+1-1
     if(ji.eq.0) exit
     n=ia + indx(ji) - 1
     if(red(n).lt.syncmin) exit
     if(k.lt.200) k=k+1
     candidate0(1,k)=n*df
     candidate0(2,k)=(jpeak(n)-1)*tstep
     candidate0(3,k)=red(n)
  enddo
  ncand=k

! Put nfqso at top of list, and save only the best of near-dupe freqs.
  do i=1,ncand
     if(abs(candidate0(1,i)-nfqso).lt.10.0) candidate0(1,i)=-candidate0(1,i)
     if(i.ge.2) then
        do j=1,i-1
           fdiff=abs(candidate0(1,i))-abs(candidate0(1,j))
           if(abs(fdiff).lt.4.0) then
              if(candidate0(3,i).ge.candidate0(3,j)) candidate0(3,j)=0.
              if(candidate0(3,i).lt.candidate0(3,j)) candidate0(3,i)=0.
           endif
        enddo
!        write(*,3001) i,candidate0(1,i-1),candidate0(1,i),candidate0(3,i-1),  &
!             candidate0(3,i)
!3001    format(i2,4f8.1)
     endif
  enddo

  fac=20.0/maxval(s)
  s=fac*s

! Sort by sync
!  call indexx(candidate0(3,1:ncand),ncand,indx)
! Sort by frequency
  call indexx(candidate0(1,1:ncand),ncand,indx)
  k=1
!  do i=ncand,1,-1
  do i=1,ncand
     j=indx(i)
!     if( candidate0(3,j) .ge. syncmin .and. candidate0(2,j).ge.-1.5 ) then
     if( candidate0(3,j) .ge. syncmin ) then
       candidate(1,k)=abs(candidate0(1,j))
       candidate(2,k)=candidate0(2,j)
       candidate(3,k)=candidate0(3,j)
       k=k+1
     endif
  enddo
  ncand=k-1
  return
end subroutine sync8


code:
subroutine sync8d(cd0,i0,ctwk,itwk,sync)

! Compute sync power for a complex, downsampled FT8 signal.

  parameter(NP2=2812,NDOWN=60)
  complex cd0(3125)
  complex csync(0:6,32)
  complex csync2(32)
  complex ctwk(32)
  complex z1,z2,z3
  logical first
  integer icos7(0:6)
  data icos7/4,2,5,6,1,3,0/
  data first/.true./
  save first,twopi,fs2,dt2,taus,baud,csync

  p(z1)=real(z1)**2 + aimag(z1)**2          !Statement function for power

! Set some constants and compute the csync array.
  if( first ) then
    twopi=8.0*atan(1.0)
    fs2=12000.0/NDOWN                       !Sample rate after downsampling
    dt2=1/fs2                               !Corresponding sample interval
    taus=32*dt2                             !Symbol duration
    baud=1.0/taus                           !Keying rate
    do i=0,6
      phi=0.0
      dphi=twopi*icos7(i)*baud*dt2
      do j=1,32
        csync(i,j)=cmplx(cos(phi),sin(phi)) !Waveform for 7x7 Costas array
        phi=mod(phi+dphi,twopi)
      enddo
    enddo
    first=.false.
  endif

  sync=0
  do i=0,6                              !Sum over 7 Costas frequencies and
     i1=i0+i*32                         !three Costas arrays
     i2=i1+36*32
     i3=i1+72*32
     csync2=csync(i,1:32)
     if(itwk.eq.1) csync2=ctwk*csync2      !Tweak the frequency
     z1=0.
     z2=0.
     z3=0.
     if(i1.ge.1 .and. i1+31.le.NP2) z1=sum(cd0(i1:i1+31)*conjg(csync2))
     if(i2.ge.1 .and. i2+31.le.NP2) z2=sum(cd0(i2:i2+31)*conjg(csync2))
     if(i3.ge.1 .and. i3+31.le.NP2) z3=sum(cd0(i3:i3+31)*conjg(csync2))
     sync = sync + p(z1) + p(z2) + p(z3)
  enddo

  return
end subroutine sync8d


code:
subroutine twkfreq1(ca,npts,fsample,a,cb)

  complex ca(npts)
  complex cb(npts)
  complex w,wstep
  real a(5)
  data twopi/6.283185307/

! Mix the complex signal
  w=1.0
  wstep=1.0
  x0=0.5*(npts+1)
  s=2.0/npts
  do i=1,npts
     x=s*(i-x0)
     p2=1.5*x*x - 0.5
     p3=2.5*(x**3) - 1.5*x
     p4=4.375*(x**4) - 3.75*(x**2) + 0.375
     dphi=(a(1) + x*a(2) + p2*a(3) + p3*a(4) + p4*a(5)) * (twopi/fsample)
     wstep=cmplx(cos(dphi),sin(dphi))
     w=w*wstep
     cb(i)=w*ca(i)
  enddo

  return
end subroutine twkfreq1


code:
subroutine watterson(c,npts,fs,delay,fspread)

  complex c(0:npts-1)
  complex c2(0:npts-1)
  complex cs1(0:npts-1)
  complex cs2(0:npts-1)

  nonzero=0
  df=fs/npts
  if(fspread.gt.0.0) then
     do i=0,npts-1
        xx=gran()
        yy=gran()
        cs1(i)=0.707*cmplx(xx,yy)
        xx=gran()
        yy=gran()
        cs2(i)=0.707*cmplx(xx,yy)
     enddo
     call four2a(cs1,npts,1,-1,1)     !To freq domain
     call four2a(cs2,npts,1,-1,1)
     do i=0,npts-1
        f=i*df
        if(i.gt.npts/2) f=(i-npts)*df
        x=(f/(0.5*fspread))**2
        a=0.
        if(x.le.50.0) then
           a=exp(-x)
        endif
        cs1(i)=a*cs1(i)
        cs2(i)=a*cs2(i)
        if(abs(f).lt.10.0) then
           p1=real(cs1(i))**2 + aimag(cs1(i))**2
           p2=real(cs2(i))**2 + aimag(cs2(i))**2
           if(p1.gt.0.0) nonzero=nonzero+1
!           write(62,3101) f,p1,p2,db(p1+1.e-12)-60,db(p2+1.e-12)-60
!3101       format(f10.3,2f12.3,2f10.3)
        endif
     enddo
     call four2a(cs1,npts,1,1,1)     !Back to time domain
     call four2a(cs2,npts,1,1,1)
     cs1(0:npts-1)=cs1(0:npts-1)/npts
     cs2(0:npts-1)=cs2(0:npts-1)/npts
  endif

  nshift=nint(0.001*delay*fs)
  c2(0:npts-1)=cshift(c(0:npts-1),nshift)
  sq=0.
  do i=0,npts-1
     if(nonzero.gt.1) then
        c(i)=0.5*(cs1(i)*c(i) + cs2(i)*c2(i))
     else
        c(i)=0.5*(c(i) + c2(i))
     endif
     sq=sq + real(c(i))**2 + aimag(c(i))**2
!     write(61,3001) i/12000.0,c(i)
!3001 format(3f12.6)
  enddo
  rms=sqrt(sq/npts)
  c=c/rms

  return
end subroutine watterson


code:
module ft8_decode

  parameter (MAXFOX=1000)
  character*12 c2fox(MAXFOX)
  character*4  g2fox(MAXFOX)
  integer nsnrfox(MAXFOX)
  integer nfreqfox(MAXFOX)
  integer n30fox(MAXFOX)
  integer n30z
  integer nfox

  type :: ft8_decoder
     procedure(ft8_decode_callback), pointer :: callback
   contains
     procedure :: decode
  end type ft8_decoder

  abstract interface
     subroutine ft8_decode_callback (this,sync,snr,dt,freq,decoded,nap,qual)
       import ft8_decoder
       implicit none
       class(ft8_decoder), intent(inout) :: this
       real, intent(in) :: sync
       integer, intent(in) :: snr
       real, intent(in) :: dt
       real, intent(in) :: freq
       character(len=37), intent(in) :: decoded
       integer, intent(in) :: nap
       real, intent(in) :: qual
     end subroutine ft8_decode_callback
  end interface

contains

  subroutine decode(this,callback,iwave,nQSOProgress,nfqso,nftx,newdat,  &
       nutc,nfa,nfb,nexp_decode,ndepth,nagain,lft8apon,lapcqonly,napwid, &
       mycall12,mygrid6,hiscall12,hisgrid6)
!    use wavhdr
    use timer_module, only: timer
    include 'ft8/ft8_params.f90'
!    type(hdr) h

    class(ft8_decoder), intent(inout) :: this
    procedure(ft8_decode_callback) :: callback
    real s(NH1,NHSYM)
    real sbase(NH1)
    real candidate(3,200)
    real dd(15*12000)
    logical, intent(in) :: lft8apon,lapcqonly,nagain
    logical newdat,lsubtract,ldupe,bcontest
    character*12 mycall12, hiscall12
    character*6 mygrid6,hisgrid6
    integer*2 iwave(15*12000)
    integer apsym(KK)
    character datetime*13,message*22,msg37*37
    character*22 allmessages(100)
    integer allsnrs(100)
    save s,dd

    bcontest=iand(nexp_decode,128).ne.0
    this%callback => callback
    write(datetime,1001) nutc        !### TEMPORARY ###
1001 format("000000_",i6.6)

    call ft8apset(mycall12,mygrid6,hiscall12,hisgrid6,bcontest,apsym,iaptype)

    dd=iwave
    ndecodes=0
    allmessages='                      '
    allsnrs=0
    ifa=nfa
    ifb=nfb
    if(nagain) then
       ifa=nfqso-10
       ifb=nfqso+10
    endif

! For now:
! ndepth=1: no subtraction, 1 pass, belief propagation only
! ndepth=2: subtraction, 3 passes, belief propagation only
! ndepth=3: subtraction, 3 passes, bp+osd
    if(ndepth.eq.1) npass=1
    if(ndepth.ge.2) npass=3
    do ipass=1,npass
      newdat=.true.  ! Is this a problem? I hijacked newdat.
      syncmin=1.5
      if(ipass.eq.1) then
        lsubtract=.true.
        if(ndepth.eq.1) lsubtract=.false.
      elseif(ipass.eq.2) then
        n2=ndecodes
        if(ndecodes.eq.0) cycle
        lsubtract=.true.
      elseif(ipass.eq.3) then
        if((ndecodes-n2).eq.0) cycle
        lsubtract=.false.
      endif

      call timer('sync8   ',0)
      call sync8(dd,ifa,ifb,syncmin,nfqso,s,candidate,ncand,sbase)
      call timer('sync8   ',1)
      do icand=1,ncand
        sync=candidate(3,icand)
        f1=candidate(1,icand)
        xdt=candidate(2,icand)
        xbase=10.0**(0.1*(sbase(nint(f1/3.125))-40.0))
        nsnr0=min(99,nint(10.0*log10(sync) - 25.5))    !### empirical ###
        call timer('ft8b    ',0)
        call ft8b(dd,newdat,nQSOProgress,nfqso,nftx,ndepth,lft8apon,       &
             lapcqonly,napwid,lsubtract,nagain,iaptype,mycall12,mygrid6,   &
             hiscall12,bcontest,sync,f1,xdt,xbase,apsym,nharderrors,dmin,  &
             nbadcrc,iappass,iera,msg37,xsnr)
        message=msg37(1:22)   !###
        nsnr=nint(xsnr)
        xdt=xdt-0.5
        hd=nharderrors+dmin
        call timer('ft8b    ',1)
        if(nbadcrc.eq.0) then
!           call jtmsg(message,iflag)
           if(bcontest) then
              call fix_contest_msg(mygrid6,message)
              msg37(1:22)=message
           endif
!           if(iand(iflag,31).ne.0) message(22:22)='?'
           ldupe=.false.
           do id=1,ndecodes
              if(message.eq.allmessages(id).and.nsnr.le.allsnrs(id)) ldupe=.true.
           enddo
           if(.not.ldupe) then
              ndecodes=ndecodes+1
              allmessages(ndecodes)=message
              allsnrs(ndecodes)=nsnr
           endif
!           write(81,1004) nutc,ncand,icand,ipass,iaptype,iappass,        &
!                nharderrors,dmin,hd,min(sync,999.0),nint(xsnr),          &
!                xdt,nint(f1),message
!1004          format(i6.6,2i4,3i2,i3,3f6.1,i4,f6.2,i5,2x,a22)
!           flush(81)
           if(.not.ldupe .and. associated(this%callback)) then
              qual=1.0-(nharderrors+dmin)/60.0 ! scale qual to [0.0,1.0]
              call this%callback(sync,nsnr,xdt,f1,msg37,iaptype,qual)
           endif
        endif
      enddo
  enddo
  return
  end subroutine decode

end module ft8_decode


code:
# Makefile for Windows in JTSDK-PY environment

# Re-direct stdout and stderr:     cmd.exe              bash
#                              make > junk 2>&1      make &> junk

CC = gcc
FC = gfortran

FFLAGS = -O2 -DWIN32 -fbounds-check -fno-second-underscore -Wall \
       -Wno-conversion -Wno-character-truncation
CFLAGS = -I. -DWIN32 -DWin32 -DBIGSYM -DHAVE_STRUCT_TIMESPEC

# Default rules
%.o: %.c
	${CC} ${CFLAGS} -c $<
%.o: %.f
	${FC} ${FFLAGS} -c $<
%.o: %.F
	${FC} ${FFLAGS} -c $<
%.o: %.f90
	${FC} ${FFLAGS} -c $<
%.o: %.F90
	${FC} ${FFLAGS} -c $<

all:    libftrsd.a

OBJS1 = extract2.o ftrsd2.o init_rs_int.o encode_rs_int.o decode_rs_int.o
libftrsd.a: $(OBJS1)
	ar cr libftrsd.a $(OBJS1)
	ranlib libftrsd.a
	cp libftrsd.a ..

# Build rsdtest
OBJS2 = rsdtest.o
rsdtest: $(OBJS2) ../libjt.a
	$(FC) -o rsdtest $(OBJS2) libftrsd.a ../libjt.a ../libpthreadGC2.a

ftrsd: ftrsd.o encode_rs_int.o decode_rs_int.o init_rs_int.o
	gcc -g -o $@ $^

encode_rs_int.o: encode_rs.c
	gcc -DBIGSYM=1 $(CFLAGS) -c -o $@ $^

decode_rs_int.o: decode_rs.c
	gcc -DBIGSYM=1 $(CFLAGS) -c -o $@ $^

init_rs_int.o: init_rs.c
	gcc -DBIGSYM=1 $(CFLAGS) -c -o $@ $^

.PHONY : clean

clean:
	rm -rf *.o libjt.a rsdtest ftrsd

path: lib/ftrsd/Makefile


code:
srcdir = .
prefix = /usr/local
exec_prefix=${prefix}
CC=gcc

CFLAGS=-I/usr/local/include -Wall -O2

all: encode_rs_int.o decode_rs_int.o init_rs_int.o sfrsd2.o sfrsd.o sfrsd

encode_rs_int.o: encode_rs.c
	gcc -DBIGSYM=1 $(CFLAGS) -c -o $@ $^

decode_rs_int.o: decode_rs.c
	gcc -DBIGSYM=1 $(CFLAGS) -c -o $@ $^

init_rs_int.o: init_rs.c
	gcc -DBIGSYM=1 $(CFLAGS) -c -o $@ $^

sfrsd2.o: sfrsd2.c
	gcc -DBIGSYM=1 $(CFLAGS) -c -o $@ $^

sfrsd.o: sfrsd.c
	gcc -DBIGSYM=1 $(CFLAGS) -c -o $@ $^

sfrsd: sfrsd.o encode_rs_int.o decode_rs_int.o init_rs_int.o sfrsd2.o
	gcc -g -o $@ $^

clean:
	rm -f *.o *.a sfrsd


code:
# Makefile for Windows in JTSDK-PY environment

# Re-direct stdout and stderr:     cmd.exe              bash
#                              make > junk 2>&1      make &> junk

CC = gcc
FC = gfortran

FFLAGS = -O2 -DWIN32 -fbounds-check -fno-second-underscore -Wall \
       -Wno-conversion -Wno-character-truncation
CFLAGS = -I. -DWIN32 -DWin32 -DBIGSYM -DHAVE_STRUCT_TIMESPEC

# Default rules
%.o: %.c
	${CC} ${CFLAGS} -c $<
%.o: %.f
	${FC} ${FFLAGS} -c $<
%.o: %.F
	${FC} ${FFLAGS} -c $<
%.o: %.f90
	${FC} ${FFLAGS} -c $<
%.o: %.F90
	${FC} ${FFLAGS} -c $<

all:    rsdtest

# Build rsdtest
OBJS2 = rsdtest.o extract2.o demod64b.o sfrsd3.o
rsdtest: $(OBJS2) ../libjt.a
	$(FC) -o rsdtest $(OBJS2) ../libjt.a ../libpthreadGC2.a

sfrsd: sfrsd.o encode_rs_int.o decode_rs_int.o init_rs_int.o
	gcc -g -o $@ $^

.PHONY : clean

clean:
	rm -rf *.o libjt.a rsdtest sfrsd


code:
/* Reed-Solomon decoder
 * Copyright 2002 Phil Karn, KA9Q
 * May be used under the terms of the GNU General Public License (GPL)
 * Modified by Steve Franke, K9AN, for use in a soft-symbol RS decoder
 */

#ifdef DEBUG
#include <stdio.h>
#endif

#include <stdlib.h>
#include <string.h>

#define	min(a,b)	((a) < (b) ? (a) : (b))

#ifdef FIXED
#include "fixed.h"
#elif defined(BIGSYM)
#include "int.h"
#else
#include "char.h"
#endif

int DECODE_RS(
#ifndef FIXED
              void *p,
#endif
              DTYPE *data, int *eras_pos, int no_eras, int calc_syn){

#ifndef FIXED
    struct rs *rs = (struct rs *)p;
#endif
    int deg_lambda, el, deg_omega;
    int i, j, r,k;
    DTYPE u,q,tmp,num1,num2,den,discr_r;
    DTYPE lambda[NROOTS+1];	// Err+Eras Locator poly
    static DTYPE s[51];					 // and syndrome poly
    DTYPE b[NROOTS+1], t[NROOTS+1], omega[NROOTS+1];
    DTYPE root[NROOTS], reg[NROOTS+1], loc[NROOTS];
    int syn_error, count;

    if( calc_syn ) {
        /* form the syndromes; i.e., evaluate data(x) at roots of g(x) */
        for(i=0;i<NROOTS;i++)
            s[i] = data[0];

        for(j=1;j<NN;j++){
            for(i=0;i<NROOTS;i++){
                if(s[i] == 0){
                    s[i] = data[j];
                } else {
                    s[i] = data[j] ^ ALPHA_TO[MODNN(INDEX_OF[s[i]] + (FCR+i)*PRIM)];
                }
            }
        }

        /* Convert syndromes to index form, checking for nonzero condition */
        syn_error = 0;
        for(i=0;i<NROOTS;i++){
            syn_error |= s[i];
            s[i] = INDEX_OF[s[i]];
        }

        if (!syn_error) {
            /* if syndrome is zero, data[] is a codeword and there are no
             * errors to correct. So return data[] unmodified
             */
            count = 0;
            goto finish;
        }

    }

    memset(&lambda[1],0,NROOTS*sizeof(lambda[0]));
    lambda[0] = 1;

    if (no_eras > 0) {
        /* Init lambda to be the erasure locator polynomial */
        lambda[1] = ALPHA_TO[MODNN(PRIM*(NN-1-eras_pos[0]))];
        for (i = 1; i < no_eras; i++) {
            u = MODNN(PRIM*(NN-1-eras_pos[i]));
            for (j = i+1; j > 0; j--) {
                tmp = INDEX_OF[lambda[j - 1]];
                if(tmp != A0)
                    lambda[j] ^= ALPHA_TO[MODNN(u + tmp)];
            }
        }

#if DEBUG >= 1
        /* Test code that verifies the erasure locator polynomial just constructed
         Needed only for decoder debugging. */

        /* find roots of the erasure location polynomial */
        for(i=1;i<=no_eras;i++)
            reg[i] = INDEX_OF[lambda[i]];

        count = 0;
        for (i = 1,k=IPRIM-1; i <= NN; i++,k = MODNN(k+IPRIM)) {
            q = 1;
            for (j = 1; j <= no_eras; j++)
                if (reg[j] != A0) {
                    reg[j] = MODNN(reg[j] + j);
                    q ^= ALPHA_TO[reg[j]];
                }
            if (q != 0)
                continue;
            /* store root and error location number indices */
            root[count] = i;
            loc[count] = k;
            count++;
        }
        if (count != no_eras) {
            printf("count = %d no_eras = %d\n lambda(x) is WRONG\n",count,no_eras);
            count = -1;
            goto finish;
        }
#if DEBUG >= 2
        printf("\n Erasure positions as determined by roots of Eras Loc Poly:\n");
        for (i = 0; i < count; i++)
            printf("%d ", loc[i]);
        printf("\n");
#endif
#endif
    }
    for(i=0;i<NROOTS+1;i++)
        b[i] = INDEX_OF[lambda[i]];

    /*
     * Begin Berlekamp-Massey algorithm to determine error+erasure
     * locator polynomial
     */
    r = no_eras;
    el = no_eras;
    while (++r <= NROOTS) {	/* r is the step number */
        /* Compute discrepancy at the r-th step in poly-form */
        discr_r = 0;
        for (i = 0; i < r; i++){
            if ((lambda[i] != 0) && (s[r-i-1] != A0)) {
                discr_r ^= ALPHA_TO[MODNN(INDEX_OF[lambda[i]] + s[r-i-1])];
            }
        }
        discr_r = INDEX_OF[discr_r];	/* Index form */
        if (discr_r == A0) {
            /* 2 lines below: B(x) <-- x*B(x) */
            memmove(&b[1],b,NROOTS*sizeof(b[0]));
            b[0] = A0;
        } else {
            /* 7 lines below: T(x) <-- lambda(x) - discr_r*x*b(x) */
            t[0] = lambda[0];
            for (i = 0 ; i < NROOTS; i++) {
                if(b[i] != A0)
                    t[i+1] = lambda[i+1] ^ ALPHA_TO[MODNN(discr_r + b[i])];
                else
                    t[i+1] = lambda[i+1];
            }
            if (2 * el <= r + no_eras - 1) {
                el = r + no_eras - el;
                /*
                 * 2 lines below: B(x) <-- inv(discr_r) *
                 * lambda(x)
                 */
                for (i = 0; i <= NROOTS; i++)
                    b[i] = (lambda[i] == 0) ? A0 : MODNN(INDEX_OF[lambda[i]] - discr_r + NN);
            } else {
                /* 2 lines below: B(x) <-- x*B(x) */
                memmove(&b[1],b,NROOTS*sizeof(b[0]));
                b[0] = A0;
            }
            memcpy(lambda,t,(NROOTS+1)*sizeof(t[0]));
        }
    }

    /* Convert lambda to index form and compute deg(lambda(x)) */
    deg_lambda = 0;
    for(i=0;i<NROOTS+1;i++){
        lambda[i] = INDEX_OF[lambda[i]];
        if(lambda[i] != A0)
            deg_lambda = i;
    }
    /* Find roots of the error+erasure locator polynomial by Chien search */
    memcpy(&reg[1],&lambda[1],NROOTS*sizeof(reg[0]));
    count = 0;		/* Number of roots of lambda(x) */
    for (i = 1,k=IPRIM-1; i <= NN; i++,k = MODNN(k+IPRIM)) {
        q = 1; /* lambda[0] is always 0 */
        for (j = deg_lambda; j > 0; j--){
            if (reg[j] != A0) {
                reg[j] = MODNN(reg[j] + j);
                q ^= ALPHA_TO[reg[j]];
            }
        }
        if (q != 0)
            continue; /* Not a root */
        /* store root (index-form) and error location number */
#if DEBUG>=2
        printf("count %d root %d loc %d\n",count,i,k);
#endif
        root[count] = i;
        loc[count] = k;
        /* If we've already found max possible roots,
         * abort the search to save time
         */
        if(++count == deg_lambda)
            break;
    }
    if (deg_lambda != count) {
        /*
         * deg(lambda) unequal to number of roots => uncorrectable
         * error detected
         */
        count = -1;
        goto finish;
    }
    /*
     * Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo
     * x**NROOTS). in index form. Also find deg(omega).
     */
    deg_omega = 0;
    for (i = 0; i < NROOTS;i++){
        tmp = 0;
        j = (deg_lambda < i) ? deg_lambda : i;
        for(;j >= 0; j--){
            if ((s[i - j] != A0) && (lambda[j] != A0))
                tmp ^= ALPHA_TO[MODNN(s[i - j] + lambda[j])];
        }
        if(tmp != 0)
            deg_omega = i;
        omega[i] = INDEX_OF[tmp];
    }
    omega[NROOTS] = A0;

    /*
     * Compute error values in poly-form. num1 = omega(inv(X(l))), num2 =
     * inv(X(l))**(FCR-1) and den = lambda_pr(inv(X(l))) all in poly-form
     */
    for (j = count-1; j >=0; j--) {
        num1 = 0;
        for (i = deg_omega; i >= 0; i--) {
            if (omega[i] != A0)
                num1  ^= ALPHA_TO[MODNN(omega[i] + i * root[j])];
        }
        num2 = ALPHA_TO[MODNN(root[j] * (FCR - 1) + NN)];
        den = 0;

        /* lambda[i+1] for i even is the formal derivative lambda_pr of lambda[i] */
        for (i = min(deg_lambda,NROOTS-1) & ~1; i >= 0; i -=2) {
            if(lambda[i+1] != A0)
                den ^= ALPHA_TO[MODNN(lambda[i+1] + i * root[j])];
        }
        if (den == 0) {
#if DEBUG >= 1
            printf("\n ERROR: denominator = 0\n");
#endif
            count = -1;
            goto finish;
        }
        /* Apply error to data */
        if (num1 != 0) {
            data[loc[j]] ^= ALPHA_TO[MODNN(INDEX_OF[num1] + INDEX_OF[num2] + NN - INDEX_OF[den])];
        }
    }
finish:
    if(eras_pos != NULL){
        for(i=0;i<count;i++)
            eras_pos[i] = loc[i];
    }
    return count;
}


code:
/* Reed-Solomon encoder
 * Copyright 2002, Phil Karn, KA9Q
 * May be used under the terms of the GNU General Public License (GPL)
 */
#include <string.h>

#ifdef FIXED
#include "fixed.h"
#elif defined(BIGSYM)
#include "int.h"
#else
#include "char.h"
#endif

void ENCODE_RS(
#ifndef FIXED
void *p,
#endif
DTYPE *data, DTYPE *bb){
#ifndef FIXED
  struct rs *rs = (struct rs *)p;
#endif
  int i, j;
  DTYPE feedback;

  memset(bb,0,NROOTS*sizeof(DTYPE));

  for(i=0;i<NN-NROOTS;i++){
    feedback = INDEX_OF[data[i] ^ bb[0]];
    if(feedback != A0){      /* feedback term is non-zero */
#ifdef UNNORMALIZED
      /* This line is unnecessary when GENPOLY[NROOTS] is unity, as it must
       * always be for the polynomials constructed by init_rs()
       */
      feedback = MODNN(NN - GENPOLY[NROOTS] + feedback);
#endif
      for(j=1;j<NROOTS;j++)
	bb[j] ^= ALPHA_TO[MODNN(feedback + GENPOLY[NROOTS-j])];
    }
    /* Shift */
    memmove(&bb[0],&bb[1],sizeof(DTYPE)*(NROOTS-1));
    if(feedback != A0)
      bb[NROOTS-1] = ALPHA_TO[MODNN(feedback + GENPOLY[0])];
    else
      bb[NROOTS-1] = 0;
  }
}


code:
subroutine extract2(s3,nadd,ntrials,param,msg)

  real s3(64,63)
  real tmp(4032)
  integer np1(0:7,0:7),np2(0:7,0:7),np0(0:7,0:7)
  integer ns(0:7,0:7)
  integer perr(0:7,0:7),pmr2(0:7,0:7)
  character msg*22
  integer dat4(12)
  integer mrsym(0:62),mr2sym(0:62),mrprob(0:62),mr2prob(0:62)
  integer correct(0:62)
  integer param(0:7)
  integer indx(0:62)
  logical first
  real*8 tt
  common/extcom/ntdecode
  data ndone/0/,ngood/0/,nbad/0/,first/.true./
  save

  if(first) then
     np1=0
     np2=0
     np0=0
     ns=0
     first=.false.
  endif

  nfail=0
1 call demod64a(s3,nadd,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow)

!  if(ntest.lt.50 .or. nlow.gt.20) then
!     ncount=-999                         !Flag bad data
!     go to 900
!  endif

  call chkhist(mrsym,nhist,ipk)
  if(nhist.ge.20) then
     nfail=nfail+1
     call pctile(s3,tmp,4032,50,base)     ! ### or, use ave from demod64a ?
     do j=1,63
        s3(ipk,j)=base
     enddo
     go to 1
  endif

  call graycode(mrsym,63,-1)
  call interleave63(mrsym,-1)
  call interleave63(mrprob,-1)

  call graycode(mr2sym,63,-1)
  call interleave63(mr2sym,-1)
  call interleave63(mr2prob,-1)

  nverbose=0
  call ftrsd2(mrsym,mrprob,mr2sym,mr2prob,ntrials,nverbose,correct,   &
       param,indx,tt,ntry)
  ncandidates=param(0)
  nhard=param(1)
  nsoft=param(2)
  nera=param(3)
  ngmd=param(4)
  ndone=ndone+1

  do i=1,12
     dat4(i)=correct(12-i)
  enddo

  msg='                      '
! if(nhard.ge.0) then
  if(nhard.ge.0 .and. nhard.le.42 .and. nsoft.le.32 .and.              &
       (nhard+nsoft).le.73) then
     call unpackmsg(dat4,msg) !Unpack the user message
     if(msg.eq.'VK7MO K1JT FN20       ') then
        ngood=ngood+1
        do k=0,62
           j=indx(k)
           i=(62-j)
           p1=mrprob(i)/255.0 + 1.e-10
           p2=mr2prob(i)/255.0
           ii=k/8
           jj=7.9999*p2/p1
           ns(ii,jj)=ns(ii,jj)+1
           if(correct(j).eq.mrsym(i)) then
              np1(ii,jj)=np1(ii,jj)+1
           else
              np0(ii,jj)=np0(ii,jj)+1
           endif
           if(correct(j).eq.mr2sym(i)) np2(ii,jj)=np2(ii,jj)+1
        enddo
     else
        nbad=nbad+1
     endif
  endif
  fgood=float(ngood)/ndone
  fbad=float(nbad)/ndone

  nboth=nhard+nsoft
  if(nhard.lt.0) then
     nsoft=99
     nera=99
     nboth=99
  endif
  write(*,1010) ndone,fgood,fbad,ncandidates,nhard,nsoft,nera,nboth,   &
       ntry,tt,msg
  write(32,1010) ndone,fgood,fbad,ncandidates,nhard,nsoft,nera,nboth,  &
       ntry,tt,msg
1010 format(i4,2f7.3,i7,4i3,i8,f8.1,1x,a22)
  flush(32)

  if(msg.eq.'VK7MO K1JT FN20       ') then
     write(33,1012) ndone,ngood,nbad,ntry,log10(float(1+ntry))
1012 format(4i8,f8.3)
     flush(33)
  endif

  rewind 40
  write(40,1080)
1080 format('Totals:')
  write(40,1090) ns
1090 format(8i7)
  write(40,1091)
1091 format(/'error:')
  write(40,1090) np0
  write(40,1092)
1092 format(/'sym = mrsym:')
  write(40,1090) np1
  write(40,1093)
1093 format(/'sym = mr2sym:')
  write(40,1090) np2

  write(40,1095)
1095 format(/'Probability of error:')
  perr=nint(100.0*float(np0)/(ns+0.001))
  write(40,1096) perr
1096 format(8i7)

  write(40,1097)
1097 format(/'P(mr2 correct | mr not correct) :')
  pmr2=nint(100.0*float(np2)/(np0+0.001))
  write(40,1096) pmr2
  flush(40)

  return
end subroutine extract2