code:
module jt9_decode

  type :: jt9_decoder
     procedure(jt9_decode_callback), pointer :: callback
   contains
     procedure :: decode
  end type jt9_decoder

  abstract interface
     subroutine jt9_decode_callback (this, sync, snr, dt, freq, drift, &
          decoded)
       import jt9_decoder
       implicit none
       class(jt9_decoder), intent(inout) :: this
       real, intent(in) :: sync
       integer, intent(in) :: snr
       real, intent(in) :: dt
       real, intent(in) :: freq
       integer, intent(in) :: drift
       character(len=22), intent(in) :: decoded
     end subroutine jt9_decode_callback
  end interface

contains

  subroutine decode(this,callback,ss,id2,nfqso,newdat,npts8,nfa,    &
       nfsplit,nfb,ntol,nzhsym,nagain,ndepth,nmode,nsubmode,nexp_decode)
    use timer_module, only: timer

    include 'constants.f90'
    class(jt9_decoder), intent(inout) :: this
    procedure(jt9_decode_callback) :: callback
    real ss(184,NSMAX)
    logical, intent(in) :: newdat, nagain
    character*22 msg
    real*4 ccfred(NSMAX)
    real*4 red2(NSMAX)
    logical ccfok(NSMAX)
    logical done(NSMAX)
    integer*2 id2(NTMAX*12000)
    integer*1 i1SoftSymbols(207)
    common/decstats/ntry65a,ntry65b,n65a,n65b,num9,numfano
    save ccfred,red2

    if(nexp_decode.eq.-99) stop     !Silence compiler warning
    this%callback => callback
    if(nmode.eq.9 .and. nsubmode.ge.1) then
       call decode9w(nfqso,ntol,nsubmode,ss,id2,sync,nsnr,xdt,freq,msg)
       if (associated(this%callback)) then
          ndrift=0
          call this%callback(sync,nsnr,xdt,freq,ndrift,msg)
       end if
       go to 999
    endif

    nsynced=0
    ndecoded=0
    nsps=6912                                   !Params for JT9-1
    df3=1500.0/2048.0

    tstep=0.5*nsps/12000.0                      !Half-symbol step (seconds)
    done=.false.

    nf0=0
    nf1=nfa
    if(nmode.eq.65+9) nf1=nfsplit
    ia=max(1,nint((nf1-nf0)/df3))
    ib=min(NSMAX,nint((nfb-nf0)/df3))
    lag1=-int(2.5/tstep + 0.9999)
    lag2=int(5.0/tstep + 0.9999)
    if(newdat) then
       call timer('sync9   ',0)
       call sync9(ss,nzhsym,lag1,lag2,ia,ib,ccfred,red2,ipk)
       call timer('sync9   ',1)
    endif

    nsps8=nsps/8
    df8=1500.0/nsps8
    dblim=db(864.0/nsps8) - 26.2

    ia1=1                         !quel compiler gripe
    ib1=1                         !quel compiler gripe
    do nqd=1,0,-1
       limit=5000
       ccflim=3.0
       red2lim=1.6
       schklim=2.2
       if(iand(ndepth,7).eq.2) then
          limit=10000
          ccflim=2.7
       endif
       if(iand(ndepth,7).eq.3 .or. nqd.eq.1) then
          limit=30000
          ccflim=2.5
          schklim=2.0
       endif
       if(nagain) then
          limit=100000
          ccflim=2.4
          schklim=1.8
       endif
       ccfok=.false.

       if(nqd.eq.1) then
          nfa1=nfqso-ntol
          nfb1=nfqso+ntol
          ia=max(1,nint((nfa1-nf0)/df3))
          ib=min(NSMAX,nint((nfb1-nf0)/df3))
          ccfok(ia:ib)=(ccfred(ia:ib).gt.(ccflim-2.0)) .and.               &
               (red2(ia:ib).gt.(red2lim-1.0))
          ia1=ia
          ib1=ib
       else
          nfa1=nf1
          nfb1=nfb
          ia=max(1,nint((nfa1-nf0)/df3))
          ib=min(NSMAX,nint((nfb1-nf0)/df3))
          do i=ia,ib
             ccfok(i)=ccfred(i).gt.ccflim .and. red2(i).gt.red2lim
          enddo
          ccfok(ia1:ib1)=.false.
       endif

       fgood=0.
       do i=ia,ib
          if(done(i) .or. (.not.ccfok(i))) cycle
          f=(i-1)*df3
          if(nqd.eq.1 .or.                                                   &
               (ccfred(i).ge.ccflim .and. abs(f-fgood).gt.10.0*df8)) then

             call timer('softsym ',0)
             fpk=nf0 + df3*(i-1)
             call softsym(id2,npts8,nsps8,newdat,fpk,syncpk,snrdb,xdt,    &
                  freq,drift,a3,schk,i1SoftSymbols)
             call timer('softsym ',1)

             sync=(syncpk+1)/4.0
             if(nqd.eq.1 .and. ((sync.lt.0.5) .or. (schk.lt.1.0))) cycle
             if(nqd.ne.1 .and. ((sync.lt.1.0) .or. (schk.lt.1.5))) cycle

             call timer('jt9fano ',0)
             call jt9fano(i1SoftSymbols,limit,nlim,msg)
             call timer('jt9fano ',1)

             if(sync.lt.0.0 .or. snrdb.lt.dblim-2.0) sync=0.0
             nsync=int(sync)
             if(nsync.gt.10) nsync=10
             nsnr=nint(snrdb)
             ndrift=nint(drift/df3)
             num9=num9+1

             if(msg.ne.'                      ') then
                numfano=numfano+1
                if (associated(this%callback)) then
                   call this%callback(sync,nsnr,xdt,freq,ndrift,msg)
                end if
                iaa=max(1,i-1)
                ibb=min(NSMAX,i+22)
                fgood=f
                nsynced=1
                ndecoded=1
                ccfok(iaa:ibb)=.false.
                done(iaa:ibb)=.true.
             endif
          endif
       enddo
       if(nagain) exit
    enddo

999 return
  end subroutine decode
end module jt9_decode


code:
subroutine jt9a()
  use, intrinsic :: iso_c_binding, only: c_f_pointer
  use prog_args
  use timer_module, only: timer
  use timer_impl, only: init_timer !, limtrace

  include 'jt9com.f90'

! These routines connect the shared memory region to the decoder.
  interface
     function address_jt9()
       use, intrinsic :: iso_c_binding, only: c_ptr
       type(c_ptr) :: address_jt9
     end function address_jt9
  end interface

  integer*1 attach_jt9
!  integer*1 lock_jt9,unlock_jt9
  integer size_jt9
! Multiple instances:
  character*80 mykey
  type(dec_data), pointer :: shared_data
  type(params_block) :: local_params
  logical fileExists

! Multiple instances:
  i0 = len(trim(shm_key))

  call init_timer (trim(data_dir)//'/timer.out')
!  open(23,file=trim(data_dir)//'/CALL3.TXT',status='unknown')

!  limtrace=-1                            !Disable all calls to timer()

! Multiple instances: set the shared memory key before attaching
  mykey=trim(repeat(shm_key,1))
  i0 = len(mykey)
  i0=setkey_jt9(trim(mykey))
  i1=attach_jt9()
  msdelay=10

10 inquire(file=trim(temp_dir)//'/.lock',exist=fileExists)
  if(fileExists) then
     call sleep_msec(msdelay)
     go to 10
  endif

  inquire(file=trim(temp_dir)//'/.quit',exist=fileExists)
  if(fileExists) then
     i1=detach_jt9()
     go to 999
  endif
  if(i1.eq.999999) stop                  !Silence compiler warning

  nbytes=size_jt9()
  if(nbytes.le.0) then
     print*,'jt9a: Shared memory mem_jt9 does not exist.'
     print*,"Must start 'jt9 -s <thekey>' from within WSJT-X."
     go to 999
  endif
  call c_f_pointer(address_jt9(),shared_data)
  local_params=shared_data%params !save a copy because wsjtx carries on accessing
  call flush(6)
  call timer('decoder ',0)
  call multimode_decoder(shared_data%ss,shared_data%id2,local_params,12000)
  call timer('decoder ',1)

100 inquire(file=trim(temp_dir)//'/.lock',exist=fileExists)
  if(fileExists) go to 10
  call sleep_msec(msdelay)
  go to 100

999 call timer('decoder ',101)

  return
end subroutine jt9a

license: mit


code:
program jt9code

! Generate simulated data for testing of WSJT-X

  character*22 msg,msgchk,msg0,msg1,decoded,cok*3,bad*1,msgtype*10,expected
  integer*4 i4tone(85)                     !Channel symbols (values 0-8)
  include 'testmsg.f90'
  include 'jt9sync.f90'

  nargs=iargc()
  if(nargs.ne.1) then
     print*,'Usage: jt9code "message"'
     print*,'       jt9code -t'
     go to 999
  endif

  call getarg(1,msg)
  nmsg=1
  if(msg(1:2).eq."-t") then
     if (NTEST+5 > MAXTEST) then
        write(*,*) "NTEST exceed MAXTEST"
     endif
     testmsg(NTEST+1)="KA1ABC WB9XYZ EN34 OOO"
     testmsg(NTEST+2)="KA1ABC WB9XYZ OOO"
     testmsg(NTEST+3)="RO"
     testmsg(NTEST+4)="RRR"
     testmsg(NTEST+5)="73"
     testmsgchk(NTEST+1)="KA1ABC WB9XYZ EN34 OOO"
     testmsgchk(NTEST+2)="KA1ABC WB9XYZ OOO"
     testmsgchk(NTEST+3)="RO"
     testmsgchk(NTEST+4)="RRR"
     testmsgchk(NTEST+5)="73"
     nmsg=NTEST+5
  endif

  write(*,1010)
1010 format("    Message                Decoded              Err? Type          Expected"/   &
            76("-"))
  do imsg=1,nmsg
     if(nmsg.gt.1) then
        msg=testmsg(imsg)
        msgchk=testmsgchk(imsg)
     endif
     call fmtmsg(msg,iz)                !To upper case, collapse multiple blanks
     msg0=msg
     ichk=0
     call chkmsg(msg,cok,nspecial,flip)     !See if it includes "OOO" report
     msg1=msg                               !Message without "OOO"

     if(nspecial.gt.0) then                  !or is a shorthand message
        if(nspecial.eq.2) decoded="RO"
        if(nspecial.eq.3) decoded="RRR"
        if(nspecial.eq.4) decoded="73"
        itype=-1
        msgtype="Shorthand"
        go to 10
     endif

     call gen9(msg,ichk,decoded,i4tone,itype)   !Encode message into tone #s

     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"

     if(cok.eq."OOO") decoded(20:22)=cok
     call fmtmsg(decoded,iz)

10   bad=" "
     expected = 'EXACT'
     if (msg0.ne.msgchk) expected = 'TRUNCATED'
     if (nmsg.eq.1) expected = 'UNKNOWN'
     if(decoded.ne.msgchk) bad="*"
     write(*,1020) imsg,msg0,decoded,bad,itype,msgtype,expected
1020 format(i2,'.',1x,a22,1x,a22,1x,a1,i3,":",a10,2x,a22)
  enddo

  if(nmsg.eq.1) write(*,1030) i4tone
1030 format(/'Channel symbols'/(30i2))

999 end program jt9code


code:
  use, intrinsic :: iso_c_binding, only: c_int, c_short, c_float, c_char, c_bool

  include 'constants.f90'

  !
  ! these structures must be kept in sync with ../commons.h
  !
  type, bind(C) :: params_block
     integer(c_int) :: nutc
     logical(c_bool) :: ndiskdat
     integer(c_int) :: ntr
     integer(c_int) :: nQSOProgress ! See MainWindow::m_QSOProgress for values
     integer(c_int) :: nfqso
     integer(c_int) :: nftx
     logical(c_bool) :: newdat
     integer(c_int) :: npts8
     integer(c_int) :: nfa
     integer(c_int) :: nfsplit
     integer(c_int) :: nfb
     integer(c_int) :: ntol
     integer(c_int) :: kin
     integer(c_int) :: nzhsym
     integer(c_int) :: nsubmode
     logical(c_bool) :: nagain
     integer(c_int) :: ndepth
     logical(c_bool) :: lft8apon
     logical(c_bool) :: lapcqonly
     logical(c_bool) :: ljt65apon
     integer(c_int) :: napwid
     integer(c_int) :: ntxmode
     integer(c_int) :: nmode
     integer(c_int) :: minw
     logical(c_bool) :: nclearave
     integer(c_int) :: minsync
     real(c_float) :: emedelay
     real(c_float) :: dttol
     integer(c_int) :: nlist
     integer(c_int) :: listutc(10)
     integer(c_int) :: n2pass
     integer(c_int) :: nranera
     integer(c_int) :: naggressive
     logical(c_bool) :: nrobust
     integer(c_int) :: nexp_decode
     character(kind=c_char) :: datetime(20)
     character(kind=c_char) :: mycall(12)
     character(kind=c_char) :: mygrid(6)
     character(kind=c_char) :: hiscall(12)
     character(kind=c_char) :: hisgrid(6)
  end type params_block

  type, bind(C) :: dec_data
     real(c_float) :: ss(184,NSMAX)
     real(c_float) :: savg(NSMAX)
     real(c_float) :: sred(5760)
     integer(c_short) :: id2(NMAX)
     type(params_block) :: params
  end type dec_data


code:
program jt9sim

! Generate simulated data for testing of WSJT-X

  use wavhdr
  use packjt
  parameter (NTMAX=120)
  parameter (NMAX=NTMAX*12000)
  type(hdr) h
  integer*2 iwave(NMAX)                  !Generated waveform (no noise)
  real*4 dat(NMAX)
  real*8 f0,f,dt,twopi,phi,dphi,baud,fspan,fsample,freq
  character msg*22,msg0*22,message*22,msgsent*22,arg*8,fname*11

  integer*4 i4tone(85)             !Channel symbols (values 0-8)
  integer*4 i4DataSymNoGray(69)    !Data Symbols, values 0-7
  integer*1 i1ScrambledBits(207)   !Unpacked bits, scrambled order
  integer*1 i1Bits(207)            !Encoded information-carrying bits
  integer*1 i1SoftSymbols(207)
  include 'jt9sync.f90'

  nargs=iargc()
  if(nargs.ne.6) then
     print*,'Usage:   jt9sim "message"     fspan nsigs minutes SNR nfiles'
     print*,'Example: jt9sim "CQ K1ABC FN42" 200  20      2    -28    1'
     print*,' '
     print*,'Enter message = "" to use entries in msgs.txt.'
     print*,'Enter SNR = 0 to generate a range of SNRs.'
     print*,'Enter SNR = 99 to generate a noiseless signal at frequency fspan'
     go to 999
  endif

  call getarg(1,msg0)
  call fmtmsg(msg0,iz)
  message=msg0                       !Transmitted message
  call getarg(2,arg)
  read(arg,*) fspan                  !Total freq range (Hz)
  call getarg(3,arg)
  read(arg,*) nsigs                  !Number of signals in each file
  call getarg(4,arg)
  read(arg,*) minutes                !Length of file (1 2 5 10 30 minutes)
  call getarg(5,arg)
  read(arg,*) snrdb                  !S/N in dB (2500 hz reference BW)
  call getarg(6,arg)
  read(arg,*) nfiles                 !Number of files

  rmsdb=25.
  rms=10.0**(0.05*rmsdb)
  fsample=12000.d0                   !Sample rate (Hz)
  dt=1.d0/fsample                    !Sample interval (s)
  twopi=8.d0*atan(1.d0)
  npts=12000*(60*minutes-6)
  nsps=0
  if(minutes.eq.1)  nsps=6912
  if(minutes.eq.2)  nsps=15360
  if(minutes.eq.5)  nsps=40960
  if(minutes.eq.10) nsps=82944
  if(minutes.eq.30) nsps=252000
  if(nsps.eq.0) stop 'Bad value for minutes.'

  f0=1400.d0                         !Center frequency (Hz)
!  f0=3000.d0                         !Center frequency (Hz)

!  f0=1500.0
!  if(minutes.eq.5)  f0=1100.
!  if(minutes.eq.10) f0=1050.
!  if(minutes.eq.30) f0=1025.

  call init_random_seed()       ! seed Fortran RANDOM_NUMBER generator
  call sgran()                  ! see C rand generator (used in gran)

  h=default_header(12000,npts)
  k=0                                !Silence compiler warning

  if(msg0(1:3).eq.'sin') read(msg0(4:),*) sinfreq

  if(message.eq."") open(12,file='msgs.txt',status='old')

  write(*,1000)
1000 format('File  N    freq      S/N  Message'/    &
            '---------------------------------------------------')

  do ifile=1,nfiles                            !Loop over all files
     nmin=(ifile-1)*minutes
     ihr=nmin/60
     imin=mod(nmin,60)
     write(fname,1002) ihr,imin                !Create output filename
1002 format('000000_',2i2.2)
     open(10,file=fname//'.wav',access='stream',status='unknown')

     if(snrdb.lt.90) then
        do i=1,npts
           dat(i)=gran()
        enddo
     else
        dat(1:npts)=0.
     endif

     if(msg0.ne.'                      ') then
        call gen9(message,0,msgsent,i4tone,itype) !Encode message into tone #s
     endif

     rewind 12
     do isig=1,nsigs                            !Loop over multiple signals

        if(msg0.eq.'                      ') then
           read(12,1004) message                !Use pre-generated message texts
1004       format(a22)
           call gen9(message,0,msgsent,i4tone,itype)
        endif

        f=f0
        if(nsigs.gt.1) f=f0 - 0.5d0*fspan + fspan*(isig-1.d0)/(nsigs-1.d0)
        snrdbx=snrdb
!        snrdbx=snrdb + (ifile-1)*4.0
        sig=10.0**(0.05*snrdbx)
        if(snrdb.gt.90.0) sig=1.0
        write(*,1020) ifile,isig,f,snrdbx,msgsent
1020    format(i3,i4,f10.3,f7.1,2x,a22)

        phi=0.
        baud=12000.d0/nsps
        k=12000                             !Start audio at t = 1.0 s
!        f1=0.0001 * (ifile-1)
        f1=0.
!        print*,ifile-1,f1
        dphi2=0.
        ddphi2=twopi*f1*dt/60.0
        do isym=1,85
           freq=f + i4tone(isym)*baud
           if(msg0(1:3).eq.'sin') freq=sinfreq
           dphi=twopi*freq*dt + dphi2
           do i=1,nsps
              phi=phi + dphi
              dphi2=dphi2 + ddphi2
              if(phi.lt.-twopi) phi=phi+twopi
              if(phi.gt.twopi) phi=phi-twopi
              xphi=phi
              k=k+1
              dat(k)=dat(k) + sig*sin(xphi)  !Use lookup table for i*2 sin(x) ?
           enddo
        enddo
     enddo

     fac=32767.0/nsigs
     if(snrdb.ge.90.0) iwave(1:npts)=nint(fac*dat(1:npts))
     if(snrdb.lt.90.0) iwave(1:npts)=nint(rms*dat(1:npts))

     write(10) h,iwave(1:npts)
     close(10)

! We're done!  Now decode the data symbols from i4tone, as a test.
     if(msg0.ne.'                      ') then
        j=0
        do i=1,85
           if(isync(i).eq.1) cycle
           j=j+1
           i4DataSymNoGray(j)=igray(i4tone(i)-1,-1)
        enddo
        call unpackbits(i4DataSymNoGray,69,3,i1ScrambledBits)
        call interleave9(i1ScrambledBits,-1,i1Bits)

        do i=1,206
           i4=-10
           if(i1Bits(i).eq.1) i4=10
           i1SoftSymbols(i)=i4
        enddo
        limit=1000
        call jt9fano(i1SoftSymbols,limit,nlim,msg)
        if(msg.ne.msg0) print*,'Decode error: ',msg0,' ',msg
     endif
  enddo

999 end program jt9sim


code:
program jt9w

  parameter (NSMAX=6827,NZMAX=60*12000)
  real ss(184,NSMAX)
  real ref(NSMAX)
  integer*2 id2(NZMAX)
  character*12 arg
  character*22 decoded

  call getarg(1,arg)
  read(arg,*) iutc

  open(20,file='refspec.dat',status='old')
  do i=1,NSMAX
     read(20,*) j,freq,ref(i)
  enddo

  do ifile=1,999
     read(60,end=999) nutc,nfqso,ntol,ndepth,nmode,nsubmode,ss,id2
     if(nutc.ne.iutc) cycle
     ntol=100
     call decode9w(nutc,nfqso,ntol,nsubmode,ss,id2,sync,nsnr,xdt,freq,decoded)
     write(*,1100) nutc,nsnr,xdt1-1.0,nint(freq),decoded
1100 format(i4.4,i4,f5.1,i5,1x,'@',1x,a22)
     exit
  enddo

999 end program jt9w


code:
subroutine jtmsg(msg,iflag)

! Attempts to identify false decodes in JT-style messages.

! Returns iflag with sum of bits as follows:
! ------------------------------------------
!    1 Grid/Report invalid
!    2 Second callsign invalid
!    4 First callsign invalid
!    8 Very unlikely free text
!   16 Questionable free text
!    0 Message is probably OK
! ------------------------------------------

  character*22 msg,t
  character*13 w1,w2,w3,w
  character*6 bc1,bc2,bc3
  character*1 c
  logical c1ok,c2ok,c3ok,isdigit,isletter,isgrid4

! Statement functions
  isdigit(c)=(ichar(c).ge.ichar('0')) .and. (ichar(c).le.ichar('9'))
  isletter(c)=(ichar(c).ge.ichar('A')) .and. (ichar(c).le.ichar('Z'))
  isgrid4(w)=(len_trim(w).eq.4 .and.                                        &
       ichar(w(1:1)).ge.ichar('A') .and. ichar(w(1:1)).le.ichar('R') .and.  &
       ichar(w(2:2)).ge.ichar('A') .and. ichar(w(2:2)).le.ichar('R') .and.  &
       ichar(w(3:3)).ge.ichar('0') .and. ichar(w(3:3)).le.ichar('9') .and.  &
       ichar(w(4:4)).ge.ichar('0') .and. ichar(w(4:4)).le.ichar('9'))

  t=trim(msg)                         !Temporary copy of msg
  nt=len_trim(t)

! Check for standard messages
! Insert underscore in "CQ AA " to "CQ ZZ ", "CQ nnn " to make them one word.
  if(t(1:3).eq.'CQ ' .and. isletter(t(4:4)) .and.                             &
       isletter(t(5:5)) .and. t(6:6).eq.' ') t(3:3)='_'
  if(t(1:3).eq.'CQ ' .and. isdigit(t(4:4)) .and.                              &
       isdigit(t(5:5)) .and. isdigit(t(6:6)) .and. t(7:7).eq.' ') t(3:3)='_'

! Parse first three words
  w1='             '
  w2='             '
  w3='             '
  i1=index(t,' ')
  if(i1.gt.0) w1(1:i1-1)=t(1:i1-1)
  t=t(i1+1:)
  i2=index(t,' ')
  if(i2.gt.0) w2(1:i2-1)=t(1:i2-1)
  t=t(i2+1:)
  i3=index(t,' ')
  if(i3.gt.0) w3(1:i3-1)=t(1:i3-1)

  if(w1(1:3).eq.'CQ ' .or. w1(1:3).eq.'CQ_' .or. w1(1:3).eq.'DE ' .or.   &
       w1(1:4).eq.'QRZ ') then
! CQ/DE/QRZ: Should have one good callsign in w2 and maybe a grid/rpt in w3
     call chkcall(w2,bc2,c2ok)
     iflag=0
     if(.not.c2ok) iflag=iflag+2
     if(len_trim(w3).ne.0 .and. (.not.isgrid4(w3))) iflag=iflag+1
     if(w1(1:3).eq.'DE ' .and. c2ok) iflag=0
     if(iflag.eq.0) return
  endif

! Check for two calls and maybe a grid, rpt, R+rpr, RRR, or 73
  iflag=0
  call chkcall(w1,bc1,c1ok)
  call chkcall(w2,bc2,c2ok)
  if(.not.c1ok) iflag=iflag+4
  if(.not.c2ok) iflag=iflag+2
  if(len_trim(w3).ne.0 .and. (.not.isgrid4(w3)) .and.                  &
          w3(1:1).ne.'+' .and. w3(1:1).ne.'-' .and.                      &
          w3(1:2).ne.'R+' .and. w3(1:2).ne.'R-' .and.                    &
          w3(1:3).ne.'73 ' .and. w3(1:4).ne.'RRR ') iflag=iflag+1
  call chkcall(w3,bc3,c3ok)
! Allow(?) non-standard messages of the form CQ AS OC K1JT
  if(w1(1:3).eq.'CQ_'.and.isletter(w2(1:1)).and.isletter(w2(2:2)).and.   &
     w2(3:3).eq.' '.and.c3ok) iflag=0
  if(iflag.eq.0 .or. nt.gt.13) return

! Check for plausible free text

  nc=0
  np=0
  do i=1,13
     c=msg(i:i)
     if(c.ne.' ') nc=nc+1             !Number of non-blank characters
     if(c.eq.'+') np=np+1             !Number of punctuation characters
     if(c.eq.'-') np=np+1
     if(c.eq.'.') np=np+1
     if(c.eq.'/') np=np+1
     if(c.eq.'?') np=np+1
  enddo
  nb=13-nc                            !Number of blanks
  iflag=16                             !Mark as potentially questionable
  if(nc.ge.12 .or. (nc.ge.11 .and. np.gt.0)) then
     iflag=8                           !Unlikely free text, flag it
  endif

! Save messages containing some common words
  if(msg(1:3).eq.'CQ ') iflag=0
  if(index(msg,'DE ').gt.0) iflag=0
  if(index(msg,'TU ').gt.0) iflag=0
  if(index(msg,' TU').gt.0) iflag=0
  if(index(msg,'73 ').gt.0) iflag=0
  if(index(msg,' 73').gt.0) iflag=0
  if(index(msg,'TNX').gt.0) iflag=0
  if(index(msg,'THX').gt.0) iflag=0
  if(index(msg,'EQSL').gt.0) iflag=0
  if(index(msg,'LOTW').gt.0) iflag=0
  if(index(msg,'DECOD').gt.0) iflag=0
  if(index(msg,'CHK').gt.0) iflag=0
  if(index(msg,'CLK').gt.0) iflag=0
  if(index(msg,'CLOCK').gt.0) iflag=0
  if(index(msg,'LOG').gt.0) iflag=0
  if(index(msg,'QRM').gt.0) iflag=0
  if(index(msg,'QSY').gt.0) iflag=0
  if(index(msg,'TEST').gt.0) iflag=0
  if(index(msg,'CQDX').gt.0) iflag=0
  if(index(msg,'CALL').gt.0) iflag=0
  if(index(msg,'QRZ').gt.0) iflag=0
  if(index(msg,'AUTO').gt.0) iflag=0
  if(index(msg,'PHOTO').gt.0) iflag=0
  if(index(msg,'HYBRID').gt.0) iflag=0

  if(c1ok .and. w1(1:6).eq.bc1) iflag=0
  if(c2ok .and. w2(1:6).eq.bc2) iflag=0

  if(nb.ge.4) iflag=0

  return
end subroutine jtmsg


code:
void ldpc32_table_(int cw[])
{
  // Compute and return the table of 65535 codewords for the (32,16) code.

  // Array y contains the sixteen rows (columns) of the parity-check matrix
  int y[16] = { 0xd452, 0x7ecb, 0xc5d5, 0xf81c,
                0x61d7, 0x0ed8, 0xa3c5, 0x9ef9,
                0xb3bd, 0xe5b6, 0x2fcd, 0xc23a,
                0x5deb, 0xfa0e, 0x35fc, 0x1379 };

  unsigned int c[2];          /* Codeword composed of 16-bit info and 16-bit parity */

  int i,j,k;
  int aux;
  int weight(int vector);

  for(k=0; k<65536; k++) {
    c[0] = k;
    c[1] = 0;
    for (i=0; i<16; i++) {
      aux = 0;
      for (j=0; j<16; j++) {
	aux = aux ^ ((c[0] & y[i]) >> j & 1);
      }
      c[1] = (c[1] << 1) ^ aux;
    }
    cw[k]=65536*c[1] + c[0];
  }
}


code:
program ldpcsim

use, intrinsic :: iso_c_binding
use iso_c_binding, only: c_loc,c_size_t
use hashing
use packjt
parameter(NRECENT=10)
character*12 recent_calls(NRECENT)
character*22 msg,msgsent,msgreceived
character*8 arg
integer*1, allocatable ::  codeword(:), decoded(:), message(:)
integer*1, target:: i1Msg8BitBytes(10)
integer*1 i1hash(4)
integer*1 msgbits(80)
integer*4 i4Msg6BitWords(13)
integer ihash
integer nerrtot(128),nerrdec(128)
real*8, allocatable ::  lratio(:), rxdata(:), rxavgd(:)
real, allocatable :: yy(:), llr(:)
equivalence(ihash,i1hash)

do i=1,NRECENT
  recent_calls(i)='            '
enddo
nerrtot=0
nerrdec=0

nargs=iargc()
if(nargs.ne.4) then
   print*,'Usage: ldpcsim  niter   navg  #trials  s '
   print*,'eg:    ldpcsim    10     1     1000    0.75'
   return
endif
call getarg(1,arg)
read(arg,*) max_iterations
call getarg(2,arg)
read(arg,*) navg
call getarg(3,arg)
read(arg,*) ntrials
call getarg(4,arg)
read(arg,*) s

! don't count hash bits as data bits
N=128
K=72
rate=real(K)/real(N)

write(*,*) "rate: ",rate

write(*,*) "niter= ",max_iterations," navg= ",navg," s= ",s

allocate ( codeword(N), decoded(K), message(K) )
allocate ( lratio(N), rxdata(N), rxavgd(N), yy(N), llr(N) )

msg="K9AN K1JT EN50"
  call packmsg(msg,i4Msg6BitWords,itype,.false.) !Pack into 12 6-bit bytes
  call unpackmsg(i4Msg6BitWords,msgsent,.false.,'      ') !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

  ihash=nhash(c_loc(i1Msg8BitBytes),int(9,c_size_t),146)
  ihash=2*iand(ihash,32767)                   !Generate the 8-bit hash
  i1Msg8BitBytes(10)=i1hash(1)                !Hash code to byte 10
  mbit=0
  do i=1, 10
    i1=i1Msg8BitBytes(i)
    do ibit=1,8
      mbit=mbit+1
      msgbits(mbit)=iand(1,ishft(i1,ibit-8))
    enddo
  enddo
  call encode_msk144(msgbits,codeword)
  call init_random_seed()

write(*,*) "Eb/N0  SNR2500   ngood  nundetected nbadhash  sigma"
do idb = -6, 14
  db=idb/2.0-1.0
  sigma=1/sqrt( 2*rate*(10**(db/10.0)) )
  ngood=0
  nue=0
  nbadhash=0

  do itrial=1, ntrials
    rxavgd=0d0
    do iav=1,navg
      call sgran()
! Create a realization of a noisy received word
      do i=1,N
        rxdata(i) = 2.0*codeword(i)-1.0 + sigma*gran()
      enddo
      rxavgd=rxavgd+rxdata
    enddo
    rxdata=rxavgd
    nerr=0
    do i=1,N
      if( rxdata(i)*(2*codeword(i)-1.0) .lt. 0 ) nerr=nerr+1
    enddo
    nerrtot(nerr)=nerrtot(nerr)+1

! 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)
    lratio=exp(llr)
    yy=rxdata

! max_iterations is max number of belief propagation iterations
!    call ldpc_decode(lratio, decoded, max_iterations, niterations, max_dither, ndither)
!    call amsdecode(yy, max_iterations, decoded, niterations)
!    call bitflipmsk144(rxdata, decoded, niterations)
    call bpdecode144(llr, max_iterations, decoded, niterations)

! If the decoder finds a valid codeword, niterations will be .ge. 0.
    if( niterations .ge. 0 ) then
      call extractmessage144(decoded,msgreceived,nhashflag,recent_calls,nrecent)
      if( nhashflag .ne. 1 ) then
        nbadhash=nbadhash+1
      endif
      nueflag=0

! Check the message plus hash against what was sent.
      do i=1,K
        if( msgbits(i) .ne. decoded(i) ) then
          nueflag=1
        endif
      enddo
      if( nhashflag .eq. 1 .and. nueflag .eq. 0 ) then
        ngood=ngood+1
        nerrdec(nerr)=nerrdec(nerr)+1
      else if( nhashflag .eq. 1 .and. nueflag .eq. 1 ) then
        nue=nue+1;
      endif
    endif
  enddo
  snr2500=db-3.5
  write(*,"(f4.1,4x,f5.1,1x,i8,1x,i8,1x,i8,8x,f5.2)") db,snr2500,ngood,nue,nbadhash,ss

enddo

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

end program ldpcsim


code:
program ldpcsim

use, intrinsic :: iso_c_binding
use hashing
use packjt

character*22 msg,msgsent,msgreceived
character*8 arg
integer*1, allocatable ::  codeword(:), decoded(:), message(:)
real*8, allocatable ::  rxdata(:), rxavgd(:)
real, allocatable :: llr(:)
integer ihash
integer*1 hardbits(32)

nargs=iargc()
if(nargs.ne.4) then
   print*,'Usage: ldpcsim  niter navg   #trials  s '
   print*,'eg:    ldpcsim   10     1     1000    0.75'
   return
endif
call getarg(1,arg)
read(arg,*) max_iterations
call getarg(2,arg)
read(arg,*) navg
call getarg(3,arg)
read(arg,*) ntrials
call getarg(4,arg)
read(arg,*) s

K=16
N=32
!rate=real(K)/real(N)
! don't count hash bits as data bits
rate=4.0/real(N)
write(*,*) "rate: ",rate
write(*,*) "niter= ",max_iterations,"navg= ",navg," s= ",s

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

msg="K1JT K9AN"
call fmtmsg(msg,iz)
call hash(msg,22,ihash)
irpt=14
ihash=iand(ihash,4095)                 !12-bit hash
ig=16*ihash + irpt                     !4-bit report
write(*,*) irpt,ihash,ig

do i=1,16
  message(i)=iand(1,ishft(ig,1-i))
enddo
write(*,'(16i1)') message
call encode_msk40(message,codeword)
write(*,'(32i1)') codeword
call init_random_seed()

write(*,*) "Eb/N0  SNR2500   ngood  nundetected nbadhash"
do idb = 0, 30
  db=idb/2.0
  sigma=1/sqrt( 2*rate*(10**(db/10.0)) )
  ngood=0
  nue=0
  nbadhash=0

  itsum=0
  do itrial=1, ntrials
    rxavgd=0d0
    do iav=1,navg
      call sgran()
! Create a realization of a noisy received word
      do i=1,N
        rxdata(i) = 2.0*codeword(i)-1.0 + sigma*gran()
      enddo
      rxavgd=rxavgd+rxdata
    enddo
    rxdata=rxavgd

! 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
    if( s .le. 0 ) then
      ss=sigma
    else
      ss=s
    endif

    llr=2.0*rxdata/(ss*ss)

    call bpdecode40(llr, max_iterations, decoded, niterations)
! If the decoder finds a valid codeword, niterations will be .ge. 0.
    if( niterations .ge. 0 ) then
      nueflag=0
      nhashflag=0
      imsg=0
      do i=1,16
        imsg=ishft(imsg,1)+iand(1,decoded(17-i))
      enddo
      nrxrpt=iand(imsg,15)
      nrxhash=(imsg-nrxrpt)/16
      if( nrxhash .ne. ihash ) then
        nbadhash=nbadhash+1
        nhashflag=1
      endif

! Check the message plus hash against what was sent.
      do i=1,K
        if( message(i) .ne. decoded(i) ) then
          nueflag=1
        endif
      enddo

      if( nhashflag .eq. 0 .and. nueflag .eq. 0 ) then
        ngood=ngood+1
        itsum=itsum+niterations
      else if( nhashflag .eq. 0 .and. nueflag .eq. 1 ) then
        nue=nue+1;
      endif
    else
      hardbits=0
      where(llr .gt. 0) hardbits=1
!      write(*,'(32i1)') hardbits
!      write(*,'(32i1)') codeword
      isum=0
      do i=1,32
        if( hardbits(i) .ne. codeword(i) ) isum=isum+1
      enddo
!      write(*,*) 'number of errors ',isum
    endif
  enddo
  avits=real(itsum)/real(ngood+0.1)
  snr2500=db-10.0
  write(*,"(f4.1,4x,f5.1,1x,i8,1x,i8,1x,i8,1x,f8.2,1x,f8.1)") db,snr2500,ngood,nue,nbadhash,ss,avits

enddo

end program ldpcsim


code:
subroutine libration(jd,RA,Dec,xl,b)

! Compute optical libration of moon at jd: that is, the sub-observer
! point (xl,b) in selenographic coordinates.  RA and Dec are
! topocentric values.

  implicit real*8 (a-h,o-z)
  parameter (RADS=0.0174532925199433d0)
  parameter (TWOPI=6.28318530717959d0)
  real*8 jd,j2000,mjd,lambda

  j2000=2451545.0d0
  RA2000=RA
  Dec2000=Dec
  year=2000.0d0 + (jd-j2000)/365.25d0
  mjd=jd-2400000.d0
  call sla_PRECES('FK5',year,2000.d0,RA2000,Dec2000)
  call sla_EQECL(RA2000,Dec2000,mjd,lambda,beta)
  day=jd - j2000
  t = day / 36525.d0
  xi = 1.54242 * RADS
  ft = 93.2720993 + 483202.0175273 * t - .0034029 * t * t
  b= ft / 360
  a = 360 * (b - floor(b))
  if (a.lt.0.d0) a = 360 + a;
  f=a/57.2957795131d0
  omega=sla_dranrm(2.182439196d0 - t*33.7570446085d0 + t*t*3.6236526d-5)
  w = lambda - omega
  y = sin(w) * cos(beta) * cos(xi) - sin(beta) * sin(xi)
  x = cos(w) * cos(beta)
  a = sla_dranrm(atan2(y, x))
  xl = a - f
  if(xl.lt.-0.25*TWOPI) xl=xl+TWOPI             !Fix 'round the back' angles
  if(xl.gt.0.25*TWOPI)  xl=xl-TWOPI
  b = asin(-sin(w) * cos(beta) * sin(xi) - sin(beta) * cos(xi))

  return
end subroutine libration


code:
subroutine lorentzian(y,npts,a)

! Input:  y(npts); assume x(i)=i, i=1,npts
! Output: a(5)
!         a(1) = baseline
!         a(2) = amplitude
!         a(3) = x0
!         a(4) = width
!         a(5) = chisqr

  real y(npts)
  real a(5)
  real deltaa(4)

  a=0.
  df=12000.0/8192.0                               !df = 1.465 Hz
  width=0.
  ipk=0
  ymax=-1.e30
  do i=1,npts
     if(y(i).gt.ymax) then
        ymax=y(i)
        ipk=i
     endif
!     write(50,3001) i,i*df,y(i)
!3001 format(i6,2f12.3)
  enddo
!  base=(sum(y(ipk-149:ipk-50)) + sum(y(ipk+51:ipk+150)))/200.0
  base=(sum(y(1:20)) + sum(y(npts-19:npts)))/40.0
  stest=ymax - 0.5*(ymax-base)
  ssum=y(ipk)
  do i=1,50
     if(ipk+i.gt.npts) exit
     if(y(ipk+i).lt.stest) exit
     ssum=ssum + y(ipk+i)
  enddo
  do i=1,50
     if(ipk-i.lt.1) exit
     if(y(ipk-i).lt.stest) exit
     ssum=ssum + y(ipk-i)
  enddo
  ww=ssum/y(ipk)
  width=2
  t=ww*ww - 5.67
  if(t.gt.0.0) width=sqrt(t)
  a(1)=base
  a(2)=ymax-base
  a(3)=ipk
  a(4)=width

! Now find Lorentzian parameters

  deltaa(1)=0.1
  deltaa(2)=0.1
  deltaa(3)=1.0
  deltaa(4)=1.0
  nterms=4

!  Start the iteration
  chisqr=0.
  chisqr0=1.e6
  do iter=1,5
     do j=1,nterms
        chisq1=fchisq0(y,npts,a)
        fn=0.
        delta=deltaa(j)
10      a(j)=a(j)+delta
        chisq2=fchisq0(y,npts,a)
        if(chisq2.eq.chisq1) go to 10
        if(chisq2.gt.chisq1) then
           delta=-delta                      !Reverse direction
           a(j)=a(j)+delta
           tmp=chisq1
           chisq1=chisq2
           chisq2=tmp
        endif
20      fn=fn+1.0
        a(j)=a(j)+delta
        chisq3=fchisq0(y,npts,a)
        if(chisq3.lt.chisq2) then
           chisq1=chisq2
           chisq2=chisq3
           go to 20
        endif

! Find minimum of parabola defined by last three points
        delta=delta*(1./(1.+(chisq1-chisq2)/(chisq3-chisq2))+0.5)
        a(j)=a(j)-delta
        deltaa(j)=deltaa(j)*fn/3.
!          write(*,4000) iter,j,a,chisq2
!4000      format(i1,i2,4f10.4,f11.3)
     enddo
     chisqr=fchisq0(y,npts,a)
!       write(*,4000) 0,0,a,chisqr
     if(chisqr/chisqr0.gt.0.99) exit
     chisqr0=chisqr
  enddo
  a(5)=chisqr

  return
end subroutine lorentzian


code:
subroutine lpf1(dd,jz,dat,jz2)

  parameter (NFFT1=64*11025,NFFT2=32*11025)
  real dd(jz)
  real dat(jz)
  real x(NFFT1)
  complex cx(0:NFFT1/2)
  equivalence (x,cx)
  save x,cx

  fac=1.0/float(NFFT1)
  x(1:jz)=fac*dd(1:jz)
  x(jz+1:NFFT1)=0.0
  call four2a(x,NFFT1,1,-1,0)                    !Forwarxd FFT, r2c
  cx(NFFT2/2:)=0.0

!  df=11025.0/NFFT1
!  do i=1,NFFT1/2
!     sx=real(cx(i))**2 + aimag(cx(i))**2
!     write(50,3000) i*df,sx
!3000 format(f15.6,e12.3)
!  enddo

  call four2a(cx,NFFT2,1,1,-1)                   !Inverse FFT, c2r
  jz2=jz/2
  dat(1:jz2)=x(1:jz2)

  return
end subroutine lpf1


code:
subroutine makepings(pings,npts,width,sig)

  real pings(npts)
  real*8 t
  real t0(14)

  iping0=-999
  dt=1.0/12000.0
  do i=1,14
     t0(i)=i                            !Make pings at t=1, 2, ... 14 s.
  enddo
  amp=sig

  do i=1,npts
     iping=min(max(1,i/12000),14)
     t=(i*dt-t0(iping))/width
     if(t.lt.0.d0 .and. t.lt.10.0) then    !????
        fac=0.
     else
        fac=2.718*t*dexp(-t)
     endif
     pings(i)=fac*amp
  enddo

  return
end subroutine makepings


code:
subroutine mixlpf(x1,nbfo,c0)

  real*4 x1(512)
  real*8 twopi,phi,dphi
  complex c1(512),c2(105+512)
  complex c0(64)
  data phi/0.d0/
  save phi,c2

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

  do i=1,512
     phi=phi+dphi
     if(phi.gt.twopi) phi=phi-twopi
     xphi=phi
     c1(i)=x1(i)*cmplx(cos(xphi),sin(xphi))
  enddo
  c2(106:105+512)=c1

  call fil3c(c2,105+512,c0,n2)
  c2(1:105)=c1(512-104:512)                       !Save 105 trailing samples

  return
end subroutine mixlpf


code:
subroutine MoonDopJPL(nyear,month,nday,uth4,lon4,lat4,RAMoon4,        &
     DecMoon4,LST4,HA4,AzMoon4,ElMoon4,vr4,dist4)

  implicit real*8 (a-h,o-z)
  real*4 uth4                    !UT in hours
  real*4 lon4                    !East longitude, degrees
  real*4 lat4                    !Latitude, degrees
  real*4 RAMoon4                 !Topocentric RA of moon, hours
  real*4 DecMoon4                !Topocentric Dec of Moon, degrees
  real*4 LST4                    !Locat sidereal time, hours
  real*4 HA4                     !Local Hour angle, degrees
  real*4 AzMoon4                 !Topocentric Azimuth of moon, degrees
  real*4 ElMoon4                 !Topocentric Elevation of moon, degrees
  real*4 vr4                     !Radial velocity of moon wrt obs, km/s
  real*4 dist4                   !Echo time, seconds

  twopi=8.d0*atan(1.d0)          !Define some constants
  rad=360.d0/twopi
  clight=2.99792458d5

  call sla_CLDJ(nyear,month,nday,djutc,j)
  djutc=djutc + uth4/24.d0
  dut=-0.460d0

  east_long=lon4/rad
  geodetic_lat=lat4/rad
  height=40.
  nspecial=0

  call ephem(djutc,dut,east_long,geodetic_lat,height,nspecial,    &
       RA,Dec,Az,El,techo,dop,fspread_1GHz,vr)

  RAMoon4=RA
  DecMoon4=Dec
  LST4=0.                 !These two variables not presently used
  HA4=0.
  AzMoon4=Az*rad
  ElMoon4=El*rad
  vr4=vr
  dist4=techo

  return
end subroutine MoonDopJPL


code:
subroutine move(x,y,n)
  real x(n),y(n)
  do i=1,n
     y(i)=x(i)
  enddo
  return
end subroutine move


code:
W1AAA K2BBB EM00
W2CCC K3DDD EM01
W3EEE K4FFF EM02
W5GGG K6HHH EM03
W7III K8JJJ EM04
W9KKK K0LLL EM05
G0MMM F1NNN JN06
G2OOO F3PPP JN07
G4QQQ F5RRR JN08
G6SSS F7TTT JN09
W1XAA K2XBB EM10
W2XCC K3XDD EM11
W3XEE K4XFF EM12
W5XGG K6XHH EM13
W7XII K8XJJ EM14
W9XKK K0XLL EM15
G0XMM F1XNN JN16
G2XOO F3XPP JN17
G4XQQ F5XRR JN18
G6XSS F7XTT JN19
W1YAA K2YBB EM20
W2YCC K3YDD EM21
W3YEE K4YFF EM22
W5YGG K6YHH EM23
W7YII K8YJJ EM24
W9YKK K0YLL EM25
G0YMM F1YNN JN26
G2YOO F3YPP JN27
G4YQQ F5YRR JN28
G6YSS F7YTT JN29
W1ZAA K2ZBB EM30
W2ZCC K3ZDD EM31
W3ZEE K4ZFF EM32
W5ZGG K6ZHH EM33
W7ZII K8ZJJ EM34
W9ZKK K0ZLL EM35
G0ZMM F1ZNN JN36
G2ZOO F3ZPP JN37
G4ZQQ F5ZRR JN38
G6ZSS F7ZTT JN39
W1AXA K2BXB EM40
W2CXC K3DXD EM41
W3EXE K4FXF EM42
W5GXG K6HXH EM43
W7IXI K8JXJ EM44
W9KXK K0LXL EM45
G0MXM F1NXN JN46
G2OXO F3PXP JN47
G4QXQ F5RXR JN48
G6SXS F7TXT JN49
W1AYA K2BYB EM50
W2CYC K3DYD EM51
W3EYE K4FYF EM52
W5GYG K6HYH EM53
W7IYI K8JYJ EM54
W9KYK K0LYL EM55
G0MYM F1NYN JN56
G2OYO F3PYP JN57
G4QYQ F5RYR JN58
G6SYS F7TYT JN59


code:
subroutine msk144_freq_search(cdat,fc,if1,if2,delf,nframes,navmask,cb,    &
     cdat2,xmax,bestf,cs,xccs)

  parameter (NSPM=864,NZ=7*NSPM)
  complex cdat(NZ)
  complex cdat2(NZ)
  complex c(NSPM)                    !Coherently averaged complex data
  complex ct2(2*NSPM)
  complex cs(NSPM)
  complex cb(42)                     !Complex waveform for sync word
  complex cc(0:NSPM-1)
  real xcc(0:NSPM-1)
  real xccs(0:NSPM-1)
  integer navmask(nframes)           !Tells which frames to average

  navg=sum(navmask)
  n=nframes*NSPM
  fac=1.0/(48.0*sqrt(float(navg)))

  do ifr=if1,if2                     !Find freq that maximizes sync
     ferr=ifr*delf
     call tweak1(cdat,n,-(fc+ferr),cdat2)
     c=0
     sumw=0.
     do i=1,nframes
        ib=(i-1)*NSPM+1
        ie=ib+NSPM-1
        if(navmask(i).eq.1) c=c + cdat2(ib:ie)
     enddo

     cc=0
     ct2(1:NSPM)=c
     ct2(NSPM+1:2*NSPM)=c

     do ish=0,NSPM-1
        cc(ish)=dot_product(ct2(1+ish:42+ish)+ct2(337+ish:378+ish),cb(1:42))
     enddo

     xcc=abs(cc)
     xb=maxval(xcc)*fac
     if(xb.gt.xmax) then
        xmax=xb
        bestf=ferr
        cs=c
        xccs=xcc
     endif
  enddo

  return
end subroutine msk144_freq_search


code:
program msk144code

! Provides examples of message packing, bit and symbol ordering,
! LDPC encoding, and other necessary details of the MSK144 protocol.

  use packjt
  character msg*22,msgsent*22,decoded,bad*1,msgtype*13,mygrid*6
  integer*4 i4tone(144)
  logical*1 bcontest
  include 'testmsg.f90'
  data mygrid/'FN20qi'/

  nargs=iargc()
  if(nargs.ne.1) then
     print*,'Usage: msk144code "message"'
     print*,'       msk144code -t'
     print*,' '
     print*,'Examples:'
     print*,'       msk144code "KA1ABC WB9XYZ EN37"'
     print*,'       msk144code "<KA1ABC WB9XYZ> R-03"'
     print*,'       msk144code "KA1ABC WB9XYZ R EN37"'
     go to 999
  endif

  call getarg(1,msg)
  nmsg=1
  if(msg(1:2).eq."-t") then
     testmsg(NTEST+1)="<KA1ABC WB9XYZ> -03"
     testmsg(NTEST+2)="<KA1ABC WB9XYZ> R+03"
     testmsg(NTEST+3)="<KA1ABC WB9XYZ> RRR"
     testmsg(NTEST+4)="<KA1ABC WB9XYZ> 73"
     testmsg(NTEST+5)="KA1ABC WB9XYZ R EN37"
     nmsg=NTEST+5
  endif

  write(*,1010)
1010 format("     Message                 Decoded                Err? Type"/   &
            74("-"))
  do imsg=1,nmsg
     if(nmsg.gt.1) msg=testmsg(imsg)
     call fmtmsg(msg,iz)                !To upper case, collapse multiple blanks
     i1=len(trim(msg))-5
     bcontest=.false.
     if(msg(i1:i1+1).eq.'R ') bcontest=.true.
     ichk=0
     call genmsk144(msg,mygrid,ichk,bcontest,msgsent,i4tone,itype)

     msgtype=""
     if(itype.eq.1) msgtype="Std Msg"
     if(itype.eq.2) msgtype="Type 1 prefix"
     if(itype.eq.3) msgtype="Type 1 suffix"
     if(itype.eq.4) msgtype="Type 2 prefix"
     if(itype.eq.5) msgtype="Type 2 suffix"
     if(itype.eq.6) msgtype="Free text"
     if(itype.eq.7) msgtype="Hashed calls"

     bad=" "
     if(msgsent.ne.msg) bad="*"
     write(*,1020) imsg,msg,msgsent,bad,itype,msgtype
1020 format(i2,'.',2x,a22,2x,a22,3x,a1,i3,": ",a13)
  enddo

  if(nmsg.eq.1) then
     n=144
     if(msg(1:1).eq."<") n=40
     write(*,1030) i4tone(1:n)
1030 format(/'Channel symbols'/(72i1))
  endif

999 end program msk144code


code:
program msk144d2

  ! Test the msk144 decoder for WSJT-X

  use options
  use timer_module, only: timer
  use timer_impl, only: init_timer
  use readwav

  character c
  character*80 line
  character*512 datadir
  character*500 infile
  character*12 mycall,hiscall
  character*6 mygrid
  character(len=500) optarg

  logical :: display_help=.false.
  logical*1 bShMsgs
  logical*1 bcontest
  logical*1 btrain
  logical*1 bswl

  type(wav_header) :: wav

  integer*2 id2(30*12000)
  integer*2 ichunk(7*1024)

  real*8 pcoeffs(5)

  type (option) :: long_options(9) = [ &
       option ('ndepth',.true.,'c','ndepth',''), &
       option ('dxcall',.true.,'d','hiscall',''), &
       option ('evemode',.true.,'e','Must be used with -s.',''), &
       option ('frequency',.true.,'f','rxfreq',''), &
       option ('help',.false.,'h','Display this help message',''), &
       option ('mycall',.true.,'m','mycall',''), &
       option ('nftol',.true.,'n','nftol',''), &
       option ('rxequalize',.false.,'r','Rx Equalize',''), &
       option ('short',.false.,'s','enable Sh','') &
       ]
  t0=0.0
  ndepth=3
  ntol=100
  nrxfreq=1500
  mycall=''
  mygrid='EN50WC'
  hiscall=''
  bShMsgs=.false.
  bcontest=.false.
  btrain=.false.
  bswl=.false.
  datadir='.'
  pcoeffs=0.d0

  do
     call getopt('c:d:ef:hm:n:rs',long_options,c,optarg,narglen,nstat,noffset,nremain,.true.)
     if( nstat .ne. 0 ) then
        exit
     end if
     select case (c)
     case ('c')
        read (optarg(:narglen), *) ndepth
     case ('d')
        read (optarg(:narglen), *) hiscall
     case ('e')
        bswl=.true.
     case ('f')
        read (optarg(:narglen), *) nrxfreq
     case ('h')
        display_help = .true.
     case ('m')
        read (optarg(:narglen), *) mycall
     case ('n')
        read (optarg(:narglen), *) ntol
     case ('r')
        btrain=.true.
     case ('s')
        bShMsgs=.true.
     end select
  end do

  if(display_help .or. nstat.lt.0 .or. nremain.lt.1) then
     print *, ''
     print *, 'Usage: msk144d [OPTIONS] file1 [file2 ...]'
     print *, ''
     print *, '       msk144 decode pre-recorded .WAV file(s)'
     print *, ''
     print *, 'OPTIONS:'
     do i = 1, size (long_options)
        call long_options(i) % print (6)
     end do
     go to 999
  endif

  call init_timer ('timer.out')
  call timer('msk144  ',0)
  ndecoded=0
  do ifile=noffset+1,noffset+nremain
     call get_command_argument(ifile,optarg,narglen)
     infile=optarg(:narglen)
     call timer('read    ',0)
     call wav%read (infile)
     i1=index(infile,'.wav')
     if( i1 .eq. 0 ) i1=index(infile,'.WAV')
     read(infile(i1-6:i1-1),*,err=998) nutc
     inquire(FILE=infile,SIZE=isize)
     npts=min((isize-216)/2,360000)
     read(unit=wav%lun) id2(1:npts)
     close(unit=wav%lun)
     call timer('read    ',1)

     do i=1,npts-7*1024+1,7*512
       ichunk=id2(i:i+7*1024-1)
       tsec=(i-1)/12000.0
       tt=sum(float(abs(id2(i:i+7*512-1))))
       if( tt .ne. 0.0 ) then
         call mskrtd(ichunk,nutc,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall,bShMsgs, &
                     bcontest,btrain,pcoeffs,bswl,datadir,line)
         if( index(line,"&") .ne. 0 .or.   &
              index(line,"^") .ne. 0 .or.   &
              index(line,"!") .ne. 0 .or.   &
              index(line,"@") .ne. 0 ) then
           write(*,*) line
         endif
       endif
     enddo
  enddo

  call timer('msk144  ',1)
  call timer('msk144  ',101)
  go to 999

998 print*,'Cannot read from file:'
  print*,infile

999 continue
end program msk144d2


code:
subroutine msk144decodeframe(c,softbits,msgreceived,nsuccess,recent_calls,nrecent)
!  use timer_module, only: timer

  parameter (NSPM=864)
  character*22 msgreceived
  character*12 recent_calls(nrecent)
  complex cb(42)
  complex cfac,cca,ccb
  complex c(NSPM)
  integer*1 decoded(80)
  integer s8(8),hardbits(144)
  real*8 dt, fs, pi, twopi
  real cbi(42),cbq(42)
  real pp(12)
  real softbits(144)
  real llr(128)
  logical first
  data first/.true./
  data s8/0,1,1,1,0,0,1,0/
  save first,cb,fs,pi,twopi,dt,s8,pp

  if(first) then
! define half-sine pulse and raised-cosine edge window
     pi=4d0*datan(1d0)
     twopi=8d0*datan(1d0)
     fs=12000.0
     dt=1.0/fs

     do i=1,12
       angle=(i-1)*pi/12.0
       pp(i)=sin(angle)
     enddo

! define the sync word waveforms
     s8=2*s8-1
     cbq(1:6)=pp(7:12)*s8(1)
     cbq(7:18)=pp*s8(3)
     cbq(19:30)=pp*s8(5)
     cbq(31:42)=pp*s8(7)
     cbi(1:12)=pp*s8(2)
     cbi(13:24)=pp*s8(4)
     cbi(25:36)=pp*s8(6)
     cbi(37:42)=pp(1:6)*s8(8)
     cb=cmplx(cbi,cbq)
     first=.false.
  endif

  nsuccess=0
  msgreceived=' '

! Estimate carrier phase.
  cca=sum(c(1:1+41)*conjg(cb))
  ccb=sum(c(1+56*6:1+56*6+41)*conjg(cb))
  cfac=ccb*conjg(cca)
  phase0=atan2(imag(cca+ccb),real(cca+ccb))

! Remove phase error - want constellation rotated so that sample points lie on I/Q axes
  cfac=cmplx(cos(phase0),sin(phase0))
  c=c*conjg(cfac)

! matched filter -
  softbits(1)=sum(imag(c(1:6))*pp(7:12))+sum(imag(c(864-5:864))*pp(1:6))
  softbits(2)=sum(real(c(1:12))*pp)
  do i=2,72
    softbits(2*i-1)=sum(imag(c(1+(i-1)*12-6:1+(i-1)*12+5))*pp)
    softbits(2*i)=sum(real(c(7+(i-1)*12-6:7+(i-1)*12+5))*pp)
  enddo

! sync word hard error weight is used as a discriminator for
! frames that have reasonable probability of decoding
  hardbits=0
  do i=1,144
    if( softbits(i) .ge. 0.0 ) then
      hardbits(i)=1
    endif
  enddo
  nbadsync1=(8-sum( (2*hardbits(1:8)-1)*s8 ) )/2
  nbadsync2=(8-sum( (2*hardbits(1+56:8+56)-1)*s8 ) )/2
  nbadsync=nbadsync1+nbadsync2
  if( nbadsync .gt. 4 ) then
    return
  endif

! normalize the softsymbols before submitting to decoder
  sav=sum(softbits)/144
  s2av=sum(softbits*softbits)/144
  ssig=sqrt(s2av-sav*sav)
  softbits=softbits/ssig

  sigma=0.75
  llr(1:48)=softbits(9:9+47)
  llr(49:128)=softbits(65:65+80-1)
  llr=2.0*llr/(sigma*sigma)

  max_iterations=10
!  call timer('bpdec144 ',0)
  call bpdecode144(llr,max_iterations,decoded,niterations)
!  call timer('bpdec144 ',1)
  if( niterations .ge. 0.0 ) then
    call extractmessage144(decoded,msgreceived,nhashflag,recent_calls,nrecent)
    if( nhashflag .gt. 0 ) then         !Hash codes match, so print it
      nsuccess=1
    endif
  endif

  return
end subroutine msk144decodeframe


code:
program msk144sd
!
! A simple decoder for slow msk144.
! Can be used as a (slow) brute-force multi-decoder by looping
! over a set of carrier frequencies.
!
  use options
  use timer_module, only: timer
  use timer_impl, only: init_timer
  use readwav

  parameter (NRECENT=10)
  parameter (NSPM=864)
  parameter (NPATTERNS=4)

  character ch
  character*80 line
  character*500 infile
  character*12 mycall,hiscall
  character*6 mygrid
  character(len=500) optarg
  character*22 msgreceived
  character*12 recent_calls(NRECENT)

  complex cdat(30*375)
  complex c(NSPM)
  complex ct(NSPM)

  real softbits(144)
  real xmc(NPATTERNS)

  logical :: display_help=.false.

  type(wav_header) :: wav

  integer iavmask(8)
  integer iavpatterns(8,NPATTERNS)
  integer npkloc(10)

  integer*2 id2(30*12000)
  integer*2 ichunk(7*1024)

  data iavpatterns/ &
       1,1,1,1,0,0,0,0, &
       0,1,1,1,1,0,0,0, &
       0,0,1,1,1,1,0,0, &
       1,1,1,1,1,1,0,0/
  data xmc/2.0,4.5,2.5,3.0/

  type (option) :: long_options(2) = [ &
       option ('frequency',.true.,'f','rxfreq',''), &
       option ('help',.false.,'h','Display this help message','') &
       ]
  t0=0.0
  ntol=100
  nrxfreq=1500

  do
     call getopt('f:h',long_options,ch,optarg,narglen,nstat,noffset,nremain,.true.)
     if( nstat .ne. 0 ) then
        exit
     end if
     select case (ch)
     case ('f')
        read (optarg(:narglen), *) nrxfreq
     case ('h')
        display_help = .true.
     end select
  end do

  if(display_help .or. nstat.lt.0 .or. nremain.lt.1) then
     print *, ''
     print *, 'Usage: msk144sd [OPTIONS] file1 [file2 ...]'
     print *, ''
     print *, '       decode pre-recorded .WAV file(s)'
     print *, ''
     print *, 'OPTIONS:'
     do i = 1, size (long_options)
        call long_options(i) % print (6)
     end do
     go to 999
  endif

  call init_timer ('timer.out')
  call timer('msk144  ',0)
  ndecoded=0
  do ifile=noffset+1,noffset+nremain
     call get_command_argument(ifile,optarg,narglen)
     infile=optarg(:narglen)
     call timer('read    ',0)
     call wav%read (infile)
     i1=index(infile,'.wav')
     if( i1 .eq. 0 ) i1=index(infile,'.WAV')
     read(infile(i1-6:i1-1),*,err=998) nutc
     inquire(FILE=infile,SIZE=isize)
     npts=min((isize-216)/2,360000)
     read(unit=wav%lun) id2(1:npts)
     close(unit=wav%lun)
     call timer('read    ',1)

!   do if=1,89  ! brute force multi-decoder
     fo=nrxfreq
!     fo=(if-1)*25.0+300.0
     call msksddc(id2,npts,fo,cdat)
     np=npts/32
     ntol=200  ! actual ntol is ntol/32=6.25 Hz. Detection window is 12.5 Hz wide
     fc=1500.0
     call msk144spd(cdat,np,ntol,ndecodesuccess,msgreceived,fc,fest,tdec,navg,ct, &
                    softbits,recent_calls,nrecent)
     nsnr=0  ! need an snr estimate
     if( ndecodesuccess .eq. 1 ) then
       fest=fo+fest-fc   ! fudging because spd thinks input signal is at 1500 Hz
       goto 900
     endif
! If short ping decoder doesn't find a decode
     npat=NPATTERNS
     do iavg=1,npat
       iavmask=iavpatterns(1:8,iavg)
       navg=sum(iavmask)
       deltaf=4.0/real(navg)  ! search increment for frequency sync
       npeaks=4
       ntol=200
       fc=1500.0
       call msk144sync(cdat(1:6*NSPM),6,ntol,deltaf,iavmask,npeaks,fc,           &
          fest,npkloc,nsyncsuccess,xmax,c)
       if( nsyncsuccess .eq. 0 ) cycle

       do ipk=1,npeaks
         do is=1,3
           ic0=npkloc(ipk)
           if(is.eq.2) ic0=max(1,ic0-1)
           if(is.eq.3) ic0=min(NSPM,ic0+1)
           ct=cshift(c,ic0-1)
           call msk144decodeframe(ct,softbits,msgreceived,ndecodesuccess,      &
                                  recent_calls,nrecent)
           if(ndecodesuccess .gt. 0) then
              tdec=tsec+xmc(iavg)*tframe
              fest=fo+(fest-fc)/32.0
              goto 900
           endif
         enddo                         !Slicer dither
       enddo                            !Peak loop
     enddo

!   enddo
900 continue
    if( ndecodesuccess .gt. 0 ) then
       write(*,1020) nutc,nsnr,tdec,nint(fest),' % ',msgreceived,navg
1020   format(i6.6,i4,f5.1,i5,a3,a22,i4)
    endif
  enddo

  call timer('msk144  ',1)
  call timer('msk144  ',101)
  go to 999

998 print*,'Cannot read from file:'
  print*,infile

999 continue
end program msk144sd

subroutine msksddc(id2,npts,fc,cdat)

! The msk144 detector/demodulator/decoder will decode signals
! with carrier frequency, fc, in the range fN/4 +/- 0.03333*fN.
!
! For slow MSK144 with nslow=32:
!  fs=12000/32=375 Hz, fN=187.5 Hz
!
! This routine accepts input samples with fs=12000 Hz. It
! downconverts and decimates by 32 to center a signal with input carrier
! frequency fc at new carrier frequency 1500/32=46.875 Hz.
! The analytic signal is returned.

  parameter (NFFT1=30*12000,NFFT2=30*375)
  integer*2 id2(npts)
  complex cx(0:NFFT1)
  complex cdat(30*375)

  dt=1.0/12000.0
  df=1.0/(NFFT1*dt)
  icenter=int(fc/df+0.5)
  i46p875=int(46.875/df+0.5)
  ishift=icenter-i46p875
  cx=cmplx(0.0,0.0)
  cx(1:npts)=id2
  call four2a(cx,NFFT1,1,-1,1)
  cx=cshift(cx,ishift)
  cx(1)=0.5*cx(1)
  cx(2*i46p875+1:)=cmplx(0.0,0.0)
  call four2a(cx,NFFT2,1,1,1)
  cdat(1:npts/32)=cx(0:npts/32-1)/NFFT1
  return

end subroutine msksddc


code:
subroutine msk144signalquality(cframe,snr,freq,t0,softbits,msg,dxcall,       &
   btrain,datadir,nbiterrors,eyeopening,pcoeffs)

  character*22 msg,msgsent
  character*12 dxcall
  character*12 training_dxcall
  character*12 trained_dxcall
  character*6 mygrid
  character*512 pcoeff_filename
  character*8 date
  character*10 time
  character*5 zone
  character*512 datadir

  complex cframe(864)
  complex cross(864)
  complex cross_avg(864)
  complex canalytic(1024)
  complex cmodel(1024)

  integer i4tone(144)
  integer hardbits(144)
  integer msgbits(144)
  integer values(8)

  logical*1 bcontest
  logical*1 btrain
  logical*1 first
  logical*1 currently_training
  logical*1 msg_has_dxcall
  logical*1 is_training_frame

  real softbits(144)
  real waveform(0:863)
  real d(1024)
  real phase(864)
  real twopi,freq,phi,dphi0,dphi1,dphi
  real*8 x(145),y(145),pp(145),sigmay(145),a(5),chisqr
  real*8 pcoeffs(5)

  parameter (NFREQLOW=500,NFREQHIGH=2500)
  data first/.true./
  save cross_avg,wt_avg,first,currently_training,   &
       navg,tlast,training_dxcall,trained_dxcall

  if (first) then
    navg=0
    cross=cmplx(0.0,0.0)
    cross_avg=cmplx(0.0,0.0)
    wt_avg=0.0
    tlast=0.0
    trained_dxcall(1:12)=' '
    training_dxcall(1:12)=' '
    currently_training=.false.
    first=.false.
  endif

  if( (currently_training .and. (dxcall .ne. training_dxcall)) .or. &
      (navg .gt. 10 )) then !reset and retrain
    navg=0
    cross=cmplx(0.0,0.0)
    cross_avg=cmplx(0.0,0.0)
    wt_avg=0.0
    tlast=0.0
    trained_dxcall(1:12)=' '
    currently_training=.false.
    training_dxcall(1:12)=' '
    trained_dxcall(1:12)=' '
!write(*,*) 'reset to untrained state '
  endif

  indx_dxcall=index(msg,trim(dxcall))
  msg_has_dxcall = indx_dxcall .ge. 4

  if( btrain .and. msg_has_dxcall .and. (.not. currently_training) ) then
    currently_training=.true.
    training_dxcall=trim(dxcall)
    trained_dxcall(1:12)=' '
!write(*,*) 'start training on call ',training_dxcall
  endif

  if( msg_has_dxcall .and. currently_training ) then
    trained_dxcall(1:12)=' '
    training_dxcall=dxcall
  endif

! use decoded message to figure out how many bit errors in the frame
  do i=1, 144
    hardbits(i)=0
    if(softbits(i) .gt. 0 ) hardbits(i)=1
  enddo

! generate tones from decoded message
  mygrid="EN50"
  ichk=0
  bcontest=.false.
  call genmsk144(msg,mygrid,ichk,bcontest,msgsent,i4tone,itype)

! reconstruct message bits from tones
  msgbits(1)=0
  do i=1,143
    if( i4tone(i) .eq. 0 ) then
      if( mod(i,2) .eq. 1 ) then
        msgbits(i+1)=msgbits(i)
      else
        msgbits(i+1)=mod(msgbits(i)+1,2)
      endif
    else
      if( mod(i,2) .eq. 1 ) then
        msgbits(i+1)=mod(msgbits(i)+1,2)
      else
        msgbits(i+1)=msgbits(i)
      endif
    endif
  enddo

  nbiterrors=0
  do i=1,144
    if( hardbits(i) .ne. msgbits(i) ) nbiterrors=nbiterrors+1
  enddo

  nplus=0
  nminus=0
  eyetop=1
  eyebot=-1
  do i=1,144
    if( msgbits(i) .eq. 1 ) then
      if( softbits(i) .lt. eyetop ) eyetop=softbits(i)
    else
      if( softbits(i) .gt. eyebot ) eyebot=softbits(i)
    endif
  enddo
  eyeopening=eyetop-eyebot

  is_training_frame =                                                          &
      (snr.gt.5.0 .and.(nbiterrors.lt.7))    .and.                             &
      (abs(t0-tlast) .gt. 0.072)             .and.                             &
      msg_has_dxcall
  if( currently_training .and. is_training_frame ) then
      twopi=8.0*atan(1.0)
      nsym=144
      if( i4tone(41) .lt. 0 ) nsym=40
      dphi0=twopi*(freq-500)/12000.0
      dphi1=twopi*(freq+500)/12000.0
      phi=-twopi/8
      indx=0
      do i=1,nsym
        if( i4tone(i) .eq. 0 ) then
          dphi=dphi0
        else
          dphi=dphi1
        endif
        do j=1,6
          waveform(indx)=cos(phi);
          indx=indx+1
          phi=mod(phi+dphi,twopi)
        enddo
      enddo
! convert the passband waveform to complex baseband
      npts=864
      nfft=1024
      d=0
      d(1:864)=waveform(0:863)
      call analytic(d,npts,nfft,canalytic,pcoeffs,.false.) ! don't equalize the model
      call tweak1(canalytic,nfft,-freq,cmodel)
      call four2a(cframe(1:864),864,1,-1,1)
      call four2a(cmodel(1:864),864,1,-1,1)

! Cross spectra from different messages can be averaged
! as long as all messages originate from dxcall.
      cross=cmodel(1:864)*conjg(cframe)/1000.0
      cross=cshift(cross,864/2)
      cross_avg=cross_avg+10**(snr/20.0)*cross
      wt_avg=wt_avg+10**(snr/20.0)
      navg=navg+1
      tlast=t0
      phase=atan2(imag(cross_avg),real(cross_avg))
      df=12000.0/864.0
      nm=145
      do i=1,145
        x(i)=(i-73)*df/1000.0
      enddo
      y=phase((864/2-nm/2):(864/2+nm/2))
      sigmay=wt_avg/abs(cross_avg((864/2-nm/2):(864/2+nm/2)))
      mode=1
      npts=145
      nterms=5
      call polyfit(x,y,sigmay,npts,nterms,mode,a,chisqr)
      pp=a(1)+x*(a(2)+x*(a(3)+x*(a(4)+x*a(5))))
      rmsdiff=sum( (pp-phase((864/2-nm/2):(864/2+nm/2)))**2 )/145.0
!write(*,*) 'training ',navg,sqrt(chisqr),rmsdiff
      if( (sqrt(chisqr).lt.1.8) .and. (rmsdiff.lt.0.5) .and. (navg.ge.5) ) then
        trained_dxcall=dxcall
        call date_and_time(date,time,zone,values)
        write(pcoeff_filename,'(i2.2,i2.2,i2.2,"_",i2.2,i2.2,i2.2)')    &
          values(1)-2000,values(2),values(3),values(5),values(6),values(7)
        pcoeff_filename=trim(trained_dxcall)//"_"//trim(pcoeff_filename)//".pcoeff"
        l1=index(datadir,char(0))-1
        datadir(l1+1:l1+1)="/"
        pcoeff_filename=datadir(1:l1+1)//trim(pcoeff_filename)
!write(*,*) 'trained - writing coefficients to: ',pcoeff_filename
        open(17,file=pcoeff_filename,status='new')
        write(17,'(i4,2f10.2,3i5,5e25.16)') navg,sqrt(chisqr),rmsdiff,NFREQLOW,NFREQHIGH,nterms,a
        do i=1, 145
          write(17,*) x(i),pp(i),y(i),sigmay(i)
        enddo
        do i=1,864
          write(17,*) i,real(cframe(i)),imag(cframe(i)),real(cross_avg(i)),imag(cross_avg(i))
        enddo
        close(17)
        training_dxcall(1:12)=' '
        currently_training=.false.
        btrain=.false.
        navg=0
      endif
  endif

  return
  end subroutine msk144signalquality


code:
program msk144sim

  use wavhdr
  parameter (NMAX=15*12000)
  real pings(0:NMAX-1)
  real waveform(0:NMAX-1)
  character*6 mygrid
  character arg*8,msg*22,msgsent*22,fname*40
  real wave(0:NMAX-1)              !Simulated received waveform
  real*8 twopi,freq,phi,dphi0,dphi1,dphi
  type(hdr) h                          !Header for .wav file
  integer*2 iwave(0:NMAX-1)
  integer itone(144)                   !Message bits
  logical*1 bcontest
  data mygrid/"EN50wc"/

  nargs=iargc()
  if(nargs.ne.6) then
     print*,'Usage:   msk144sim       message      freq width nslow snr nfiles'
     print*,'Example: msk144sim "K1ABC W9XYZ EN37" 1500  0.12   1    2    1'
     print*,'         msk144sim "K1ABC W9XYZ EN37" 1500  2.5   32   15    1'
     go to 999
  endif
  call getarg(1,msg)
  call getarg(2,arg)
  read(arg,*) freq
  call getarg(3,arg)
  read(arg,*) width
  call getarg(4,arg)
  read(arg,*) nslow
  call getarg(5,arg)
  read(arg,*) snrdb
  call getarg(6,arg)
  read(arg,*) nfiles

!sig is the peak amplitude of the ping.
  sig=sqrt(2.0)*10.0**(0.05*snrdb)
  h=default_header(12000,NMAX)
  i1=len(trim(msg))-5
  bcontest=.false.
  if(msg(i1:i1+1).eq.'R ') bcontest=.true.
  ichk=0
  call genmsk144(msg,mygrid,ichk,bcontest,msgsent,itone,itype)
  twopi=8.d0*atan(1.d0)

  nsym=144
  nsps=6*nslow
  if( itone(41) .lt. 0 ) nsym=40
  baud=2000.d0/nslow
  dphi0=twopi*(freq-0.25d0*baud)/12000.d0
  dphi1=twopi*(freq+0.25d0*baud)/12000.d0
  phi=0.0
  k=0
  nreps=NMAX/(nsym*nsps)
  print*,nsym,nslow,nsps,baud,freq
  do jrep=1,nreps
    do i=1,nsym
      if( itone(i) .eq. 0 ) then
        dphi=dphi0
      else
        dphi=dphi1
      endif
      do j=1,nsps
        waveform(k)=cos(phi)
        k=k+1
        phi=mod(phi+dphi,twopi)
      enddo
    enddo
  enddo

  if(itype.lt.1 .or. itype.gt.7) then
     print*,'Illegal message'
     go to 999
  endif

  if(nslow.eq.1) call makepings(pings,NMAX,width,sig)

!  call sgran()
  do ifile=1,nfiles                  !Loop over requested number of files
     write(fname,1002) ifile         !Output filename
1002 format('000000_',i6.6)
     open(10,file=fname(1:13)//'.wav',access='stream',status='unknown')

     wave=0.0
     iwave=0
     fac=sqrt(6000.0/2500.0)
     do i=0,NMAX-1
        xx=gran()
        if(nslow.eq.1) wave(i)=pings(i)*waveform(i) + fac*xx
        if(nslow.gt.1) wave(i)=sig*waveform(i) + fac*xx
        iwave(i)=30.0*wave(i)
     enddo

     write(10) h,iwave               !Save the .wav file
     close(10)

  enddo

999 end program msk144sim


code:
subroutine msk144spd(cbig,n,ntol,nsuccess,msgreceived,fc,fret,tret,navg,ct,   &
  softbits,recent_calls,nrecent)

! MSK144 short-ping-decoder

  use timer_module, only: timer

  parameter (NSPM=864, MAXSTEPS=100, NFFT=NSPM, MAXCAND=5, NPATTERNS=6)
  character*22 msgreceived
  character*12 recent_calls(nrecent)
  complex cbig(n)
  complex cdat(3*NSPM)                    !Analytic signal
  complex c(NSPM)
  complex ct(NSPM)
  complex ctmp(NFFT)
  integer, dimension(1) :: iloc
  integer indices(MAXSTEPS)
  integer npkloc(10)
  integer navpatterns(3,NPATTERNS)
  integer navmask(3)
  integer nstart(MAXCAND)
  logical ismask(NFFT)
  real detmet(-2:MAXSTEPS+3)
  real detmet2(-2:MAXSTEPS+3)
  real detfer(MAXSTEPS)
  real rcw(12)
  real ferrs(MAXCAND)
  real snrs(MAXCAND)
  real softbits(144)
  real tonespec(NFFT)
  real tpat(NPATTERNS)
  real*8 dt, df, fs, pi, twopi
  logical first
  data first/.true./
  data navpatterns/ &
       0,1,0, &
       1,0,0, &
       0,0,1, &
       1,1,0, &
       0,1,1, &
       1,1,1/
  data tpat/1.5,0.5,2.5,1.0,2.0,1.5/

  save df,first,fs,pi,twopi,dt,tframe,rcw

  if(first) then
     nmatchedfilter=1
! define half-sine pulse and raised-cosine edge window
     pi=4d0*datan(1d0)
     twopi=8d0*datan(1d0)
     fs=12000.0
     dt=1.0/fs
     df=fs/NFFT
     tframe=NSPM/fs

     do i=1,12
       angle=(i-1)*pi/12.0
       rcw(i)=(1-cos(angle))/2
     enddo

     first=.false.
  endif

  ! fill the detmet, detferr arrays
  nstep=(n-NSPM)/216  ! 72ms/4=18ms steps
  detmet=0
  detmet2=0
  detfer=-999.99
  nfhi=2*(fc+500)
  nflo=2*(fc-500)
  ihlo=nint((nfhi-2*ntol)/df) + 1
  ihhi=nint((nfhi+2*ntol)/df) + 1
  illo=nint((nflo-2*ntol)/df) + 1
  ilhi=nint((nflo+2*ntol)/df) + 1
  i2000=nint(nflo/df) + 1
  i4000=nint(nfhi/df) + 1
  do istp=1,nstep
    ns=1+216*(istp-1)
    ne=ns+NSPM-1
    if( ne .gt. n ) exit
    ctmp=cmplx(0.0,0.0)
    ctmp(1:NSPM)=cbig(ns:ne)

! Coarse carrier frequency sync - seek tones at 2000 Hz and 4000 Hz in
! squared signal spectrum.

    ctmp=ctmp**2
    ctmp(1:12)=ctmp(1:12)*rcw
    ctmp(NSPM-11:NSPM)=ctmp(NSPM-11:NSPM)*rcw(12:1:-1)
    call four2a(ctmp,NFFT,1,-1,1)
    tonespec=abs(ctmp)**2

    ismask=.false.
    ismask(ihlo:ihhi)=.true.  ! high tone search window
    iloc=maxloc(tonespec,ismask)
    ihpk=iloc(1)
    deltah=-real( (ctmp(ihpk-1)-ctmp(ihpk+1)) / (2*ctmp(ihpk)-ctmp(ihpk-1)-ctmp(ihpk+1)) )
    ah=tonespec(ihpk)
    ahavp=(sum(tonespec,ismask)-ah)/count(ismask)
    trath=ah/(ahavp+0.01)
    ismask=.false.
    ismask(illo:ilhi)=.true.   ! window for low tone
    iloc=maxloc(tonespec,ismask)
    ilpk=iloc(1)
    deltal=-real( (ctmp(ilpk-1)-ctmp(ilpk+1)) / (2*ctmp(ilpk)-ctmp(ilpk-1)-ctmp(ilpk+1)) )
    al=tonespec(ilpk)
    alavp=(sum(tonespec,ismask)-al)/count(ismask)
    tratl=al/(alavp+0.01)
    fdiff=(ihpk+deltah-ilpk-deltal)*df
    ferrh=(ihpk+deltah-i4000)*df/2.0
    ferrl=(ilpk+deltal-i2000)*df/2.0
    if( ah .ge. al ) then
      ferr=ferrh
    else
      ferr=ferrl
    endif
    detmet(istp)=max(ah,al)
    detmet2(istp)=max(trath,tratl)
    detfer(istp)=ferr
  enddo  ! end of detection-metric and frequency error estimation loop

  call indexx(detmet(1:nstep),nstep,indices) !find median of detection metric vector
  xmed=detmet(indices(nstep/4))
  detmet=detmet/xmed ! noise floor of detection metric is 1.0
  ndet=0

  do ip=1,MAXCAND ! Find candidates
    iloc=maxloc(detmet(1:nstep))
    il=iloc(1)
    if( (detmet(il) .lt. 3.0) ) exit
    if( abs(detfer(il)) .le. ntol ) then
      ndet=ndet+1
      nstart(ndet)=1+(il-1)*216+1
      ferrs(ndet)=detfer(il)
      snrs(ndet)=12.0*log10(detmet(il))/2-9.0
    endif
    detmet(il)=0.0
  enddo

  if( ndet .lt. 3 ) then
    do ip=1,MAXCAND-ndet ! Find candidates
      iloc=maxloc(detmet2(1:nstep))
      il=iloc(1)
      if( (detmet2(il) .lt. 12.0) ) exit
      if( abs(detfer(il)) .le. ntol ) then
        ndet=ndet+1
        nstart(ndet)=1+(il-1)*216+1
        ferrs(ndet)=detfer(il)
        snrs(ndet)=12.0*log10(detmet2(il))/2-9.0
      endif
      detmet2(il)=0.0
    enddo
  endif

  nsuccess=0
  msgreceived=' '
  npeaks=2
  ntol0=8
  deltaf=2.0
  do icand=1,ndet  ! Try to sync/demod/decode each candidate.
    ib=max(1,nstart(icand)-NSPM)
    ie=ib-1+3*NSPM
    if( ie .gt. n ) then
      ie=n
      ib=ie-3*NSPM+1
    endif
    cdat=cbig(ib:ie)
    fo=fc+ferrs(icand)
    do iav=1,NPATTERNS
      navmask=navpatterns(1:3,iav)
      call msk144sync(cdat,3,ntol0,deltaf,navmask,npeaks,fo,fest,npkloc,   &
           nsyncsuccess,xmax,c)

      if( nsyncsuccess .eq. 0 ) cycle

      do ipk=1,npeaks
        do is=1,3
          ic0=npkloc(ipk)
          if( is.eq.2) ic0=max(1,ic0-1)
          if( is.eq.3) ic0=min(NSPM,ic0+1)
          ct=cshift(c,ic0-1)
          call msk144decodeframe(ct,softbits,msgreceived,ndecodesuccess,    &
               recent_calls,nrecent)
          if( ndecodesuccess .gt. 0 ) then
            tret=(nstart(icand)+NSPM/2)/fs
            fret=fest
            navg=sum(navmask)
            nsuccess=1
            return
          endif
        enddo
      enddo
    enddo
  enddo       ! candidate loop

  return
end subroutine msk144spd


code:
subroutine msk144sync(cdat,nframes,ntol,delf,navmask,npeaks,fc,fest,   &
     npklocs,nsuccess,xmax,c)

!$ use omp_lib

  parameter (NSPM=864)
  complex cdat(NSPM*nframes)
  complex cdat2(NSPM*nframes,8)
  complex c(NSPM)                    !Coherently averaged complex data
  complex cs(NSPM,8)
  complex cb(42)                     !Complex waveform for sync word

  integer s8(8)
  integer iloc(1)
  integer npklocs(npeaks)
  integer navmask(nframes)                 ! defines which frames to average

  real cbi(42),cbq(42)
  real pkamps(npeaks)
  real xcc(0:NSPM-1)
  real xccs(0:NSPM-1,8)
  real xm(8)
  real bf(8)
  real pp(12)                        !Half-sine pulse shape
  logical first
  data first/.true./
  data s8/0,1,1,1,0,0,1,0/
  save first,cb,fs,pi,twopi,dt,s8,pp

  if(first) then
     pi=4.0*atan(1.0)
     twopi=8.0*atan(1.0)
     fs=12000.0
     dt=1.0/fs

     do i=1,12                       !Define half-sine pulse
       angle=(i-1)*pi/12.0
       pp(i)=sin(angle)
     enddo

! Define the sync word waveforms
     s8=2*s8-1
     cbq(1:6)=pp(7:12)*s8(1)
     cbq(7:18)=pp*s8(3)
     cbq(19:30)=pp*s8(5)
     cbq(31:42)=pp*s8(7)
     cbi(1:12)=pp*s8(2)
     cbi(13:24)=pp*s8(4)
     cbi(25:36)=pp*s8(6)
     cbi(37:42)=pp(1:6)*s8(8)
     cb=cmplx(cbi,cbq)

     first=.false.
  endif

  nfreqs=2*nint(ntol/delf) + 1
  xm=0.0
  bf=0.0
  nthreads=1
!$ nthreads=min(8,int(OMP_GET_MAX_THREADS(),4))
  nstep=nfreqs/nthreads

!$OMP PARALLEL NUM_THREADS(nthreads) PRIVATE(id,if1,if2)
  id=1
!$ id=OMP_GET_THREAD_NUM() + 1            !Thread id = 1,2,...
  if1=-nint(ntol/delf) + (id-1)*nstep
  if2=if1+nstep-1
  if(id.eq.nthreads) if2=nint(ntol/delf)
  call msk144_freq_search(cdat,fc,if1,if2,delf,nframes,navmask,cb,    &
       cdat2(1,id),xm(id),bf(id),cs(1,id),xccs(1,id))
!$OMP END PARALLEL

  xmax=xm(1)
  fest=fc+bf(1)
  c=cs(1:NSPM,1)
  xcc=xccs(0:NSPM-1,1)
  if(nthreads.gt.1) then
     do i=2,nthreads
        if(xm(i).gt.xmax) then
           xmax=xm(i)
           fest=fc+bf(i)
           c=cs(1:NSPM,i)
           xcc=xccs(0:NSPM-1,i)
        endif
     enddo
  endif

! Find npeaks largest peaks
  do ipk=1,npeaks
     iloc=maxloc(xcc)
     ic2=iloc(1)
     npklocs(ipk)=ic2
     pkamps(ipk)=xcc(ic2-1)
     xcc(max(0,ic2-7):min(NSPM-1,ic2+7))=0.0
  enddo

  nsuccess=0
  if(xmax.ge.1.3) nsuccess=1

  return
end subroutine msk144sync


code:
subroutine msk40_freq_search(cdat,fc,if1,if2,delf,nframes,navmask,cb,    &
     cdat2,xmax,bestf,cs,xccs)

  parameter (NSPM=240,NZ=7*NSPM)
  complex cdat(NZ)
  complex cdat2(NZ)
  complex c(NSPM)                    !Coherently averaged complex data
  complex ct2(2*NSPM)
  complex cs(NSPM)
  complex cb(42)                     !Complex waveform for sync word
  complex cc(0:NSPM-1)
  real xcc(0:NSPM-1)
  real xccs(0:NSPM-1)
  integer navmask(nframes)           !Tells which frames to average

  navg=sum(navmask)
  n=nframes*NSPM
!  fac=1.0/(48.0*sqrt(float(navg)))
  fac=1.0/(24.0*sqrt(float(navg)))

  do ifr=if1,if2                     !Find freq that maximizes sync
     ferr=ifr*delf
     call tweak1(cdat,n,-(fc+ferr),cdat2)
     c=0
     do i=1,nframes
        ib=(i-1)*NSPM+1
        ie=ib+NSPM-1
        if( navmask(i) .eq. 1 ) c=c+cdat2(ib:ie)
     enddo

     cc=0
     ct2(1:NSPM)=c
     ct2(NSPM+1:2*NSPM)=c

     do ish=0,NSPM-1
        cc(ish)=dot_product(ct2(1+ish:42+ish),cb(1:42))
     enddo

     xcc=abs(cc)
     xb=maxval(xcc)*fac
     if(xb.gt.xmax) then
        xmax=xb
        bestf=ferr
        cs=c
        xccs=xcc
     endif
  enddo

!  write(71,3001) fc,delf,if1,if2,nframes,bestf,xmax
!3001 format(2f8.3,3i5,2f8.3)

  return
end subroutine msk40_freq_search


code:
subroutine msk40decodeframe(c,mycall,hiscall,xsnr,bswl,nhasharray,             &
                            recent_calls,nrecent,msgreceived,nsuccess)
!  use timer_module, only: timer

  parameter (NSPM=240)
  character*4 rpt(0:15)
  character*6 mycall,hiscall,mycall0,hiscall0
  character*22 hashmsg,msgreceived
  character*12 recent_calls(nrecent)
  complex cb(42)
  complex cfac,cca
  complex c(NSPM)
  integer*1 cw(32)
  integer*1 decoded(16)
  integer s8r(8),hardbits(40)
  integer nhasharray(nrecent,nrecent)
  real*8 dt, fs, pi, twopi
  real cbi(42),cbq(42)
  real pp(12)
  real softbits(40)
  real llr(32)
  logical first
  logical*1 bswl
  data first/.true./
  data s8r/1,0,1,1,0,0,0,1/
  data mycall0/'dummy'/,hiscall0/'dummy'/
  data rpt/"-03 ","+00 ","+03 ","+06 ","+10 ","+13 ","+16 ", &
           "R-03","R+00","R+03","R+06","R+10","R+13","R+16", &
           "RRR ","73  "/
  save first,cb,fs,pi,twopi,dt,s8r,pp,rpt,mycall0,hiscall0,ihash

  if(first) then
! define half-sine pulse and raised-cosine edge window
     pi=4d0*datan(1d0)
     twopi=8d0*datan(1d0)
     fs=12000.0
     dt=1.0/fs

     do i=1,12
       angle=(i-1)*pi/12.0
       pp(i)=sin(angle)
     enddo

! define the sync word waveforms
     s8r=2*s8r-1
     cbq(1:6)=pp(7:12)*s8r(1)
     cbq(7:18)=pp*s8r(3)
     cbq(19:30)=pp*s8r(5)
     cbq(31:42)=pp*s8r(7)
     cbi(1:12)=pp*s8r(2)
     cbi(13:24)=pp*s8r(4)
     cbi(25:36)=pp*s8r(6)
     cbi(37:42)=pp(1:6)*s8r(8)
     cb=cmplx(cbi,cbq)
     first=.false.
  endif

  if(mycall.ne.mycall0 .or. hiscall.ne.hiscall0) then
    hashmsg=trim(mycall)//' '//trim(hiscall)
    if( hashmsg .ne. ' ' .and. hiscall .ne. ''  ) then ! protect against blank mycall/hiscall
      call fmtmsg(hashmsg,iz)
      call hash(hashmsg,22,ihash)
      ihash=iand(ihash,4095)
    else
      ihash=9999  ! so that it can never match a received hash
    endif
    mycall0=mycall
    hiscall0=hiscall
  endif

  nsuccess=0
  msgreceived=' '

! Estimate carrier phase.
  cca=sum(c(1:1+41)*conjg(cb))
  phase0=atan2(imag(cca),real(cca))

! Remove phase error - want constellation rotated so that sample points lie on I/Q axes
  cfac=cmplx(cos(phase0),sin(phase0))
  c=c*conjg(cfac)

! Matched filter.
  softbits(1)=sum(imag(c(1:6))*pp(7:12))+sum(imag(c(NSPM-5:NSPM))*pp(1:6))
  softbits(2)=sum(real(c(1:12))*pp)
  do i=2,20
    softbits(2*i-1)=sum(imag(c(1+(i-1)*12-6:1+(i-1)*12+5))*pp)
    softbits(2*i)=sum(real(c(7+(i-1)*12-6:7+(i-1)*12+5))*pp)
  enddo

! Sync word hard error weight is used to reject frames that
! are unlikely to decode.
  hardbits=0
  do i=1,40
    if( softbits(i) .ge. 0.0 ) then
      hardbits(i)=1
    endif
  enddo
  nbadsync1=(8-sum( (2*hardbits(1:8)-1)*s8r ) )/2
  nbadsync=nbadsync1
  if( nbadsync .gt. 3 ) then
    return
  endif

! Normalize the softsymbols before submitting to decoder.
  sav=sum(softbits)/40
  s2av=sum(softbits*softbits)/40
  ssig=sqrt(s2av-sav*sav)
  softbits=softbits/ssig

  sigma=0.75
!  if(xsnr.lt.0.0) sigma=0.75-0.0875*xsnr
  if(xsnr.lt.0.0) sigma=0.75-0.11*xsnr
  llr(1:32)=softbits(9:40)
  llr=2.0*llr/(sigma*sigma)

  max_iterations=5
  call bpdecode40(llr,max_iterations,decoded,niterations)

  if( niterations .ge. 0.0 ) then
    call encode_msk40(decoded,cw)
    nhammd=0
    cord=0.0
    do i=1,32
      if( cw(i) .ne. hardbits(i+8) ) then
        nhammd=nhammd+1
        cord=cord+abs(softbits(i+8))
      endif
    enddo

    imsg=0
    do i=1,16
      imsg=ishft(imsg,1)+iand(1,decoded(17-i))
    enddo
    nrxrpt=iand(imsg,15)
    nrxhash=(imsg-nrxrpt)/16

    if(nhammd.le.4 .and. cord .lt. 0.65 .and.                                  &
       nrxhash.eq.ihash .and. nrxrpt.ge.7) then
!write(*,*) 'decodeframe 1',nbadsync,nhammd,cord,nrxhash,nrxrpt,ihash,xsnr,sigma
      nsuccess=1
      write(msgreceived,'(a1,a,1x,a,a1,1x,a4)') "<",trim(mycall),   &
                                    trim(hiscall),">",rpt(nrxrpt)
      return
    elseif(bswl .and. nhammd.le.4 .and. cord.lt.0.65 .and. nrxrpt.ge.7 ) then
      do i=1,nrecent
        do j=i+1,nrecent
          if( nrxhash .eq. nhasharray(i,j) ) then
            nsuccess=2
            write(msgreceived,'(a1,a,1x,a,a1,1x,a4)') "<",trim(recent_calls(i)),   &
                                  trim(recent_calls(j)),">",rpt(nrxrpt)
!write(*,*) 'decodeframe 2',nbadsync,nhammd,cord,nrxhash,nrxrpt,ihash,xsnr,sigma
          elseif( nrxhash .eq. nhasharray(j,i) ) then
            nsuccess=2
            write(msgreceived,'(a1,a,1x,a,a1,1x,a4)') "<",trim(recent_calls(j)),   &
                                  trim(recent_calls(i)),">",rpt(nrxrpt)
!write(*,*) 'decodeframe 3',nbadsync,nhammd,cord,nrxhash,nrxrpt,ihash,xsnr,sigma
          endif
        enddo
      enddo
      if(nsuccess.eq.0) then
        nsuccess=3
!write(*,*) 'decodeframe 4',bswl,nbadsync,nhammd,cord,nrxhash,nrxrpt,ihash,xsnr,sigma,nsuccess
        write(msgreceived,'(a1,i4.4,a1,1x,a4)') "<",nrxhash,">",rpt(nrxrpt)
      endif
    endif
  endif

  return
end subroutine msk40decodeframe


code:
subroutine msk40spd(cbig,n,ntol,mycall,hiscall,bswl,nhasharray,recent_calls,   &
                    nrecent,nsuccess,msgreceived,fc,fret,tret,navg)
! msk40 short-ping-decoder

  use timer_module, only: timer

  parameter (NSPM=240, MAXSTEPS=150, NFFT=NSPM, MAXCAND=5, NPATTERNS=6)
  character*6 mycall,hiscall
  character*22 msgreceived
  character*12 recent_calls(nrecent)
  complex cbig(n)
  complex cdat(3*NSPM)                    !Analytic signal
  complex c(NSPM)
  complex ct(NSPM)
  complex ctmp(NFFT)
  integer, dimension(1) :: iloc
  integer indices(MAXSTEPS)
  integer npkloc(10)
  integer navpatterns(3,NPATTERNS)
  integer navmask(3)
  integer nstart(MAXCAND)
  integer nhasharray(nrecent,nrecent)
  logical ismask(NFFT)
  logical*1 bswl
  real detmet(-2:MAXSTEPS+3)
  real detmet2(-2:MAXSTEPS+3)
  real detfer(MAXSTEPS)
  real rcw(12)
  real ferrs(MAXCAND)
  real snrs(MAXCAND)
  real tonespec(NFFT)
  real tpat(NPATTERNS)
  real*8 dt, df, fs, pi, twopi
  logical first
  data first/.true./
  data navpatterns/ &
       0,1,0, &
       1,0,0, &
       0,0,1, &
       1,1,0, &
       0,1,1, &
       1,1,1/
  data tpat/1.5,0.5,2.5,1.0,2.0,1.5/
  save df,first,fs,pi,twopi,dt,tframe,rcw

  if(first) then
     nmatchedfilter=1
! define half-sine pulse and raised-cosine edge window
     pi=4d0*datan(1d0)
     twopi=8d0*datan(1d0)
     fs=12000.0
     dt=1.0/fs
     df=fs/NFFT
     tframe=NSPM/fs

     do i=1,12
       angle=(i-1)*pi/12.0
       rcw(i)=(1-cos(angle))/2
     enddo

     first=.false.
  endif

  ! fill the detmet, detferr arrays
  nstep=(n-NSPM)/60  ! 20ms/4=5ms steps
  detmet=0
  detmet2=0
  detfer=-999.99
  nfhi=2*(fc+500)
  nflo=2*(fc-500)
  ihlo=nint((nfhi-2*ntol)/df)+1
  ihhi=nint((nfhi+2*ntol)/df)+1
  illo=nint((nflo-2*ntol)/df)+1
  ilhi=nint((nflo+2*ntol)/df)+1
  i2000=nint(nflo/df)+1
  i4000=nint(nfhi/df)+1
  do istp=1,nstep
    ns=1+60*(istp-1)
    ne=ns+NSPM-1
    if( ne .gt. n ) exit
    ctmp=cmplx(0.0,0.0)
    ctmp(1:NSPM)=cbig(ns:ne)

! Coarse carrier frequency sync - seek tones at 2000 Hz and 4000 Hz in
! squared signal spectrum.

    ctmp=ctmp**2
    ctmp(1:12)=ctmp(1:12)*rcw
    ctmp(NSPM-11:NSPM)=ctmp(NSPM-11:NSPM)*rcw(12:1:-1)
    call four2a(ctmp,NFFT,1,-1,1)
    tonespec=abs(ctmp)**2

    ismask=.false.
    ismask(ihlo:ihhi)=.true.  ! high tone search window
    iloc=maxloc(tonespec,ismask)
    ihpk=iloc(1)
    deltah=-real( (ctmp(ihpk-1)-ctmp(ihpk+1)) / (2*ctmp(ihpk)-ctmp(ihpk-1)-ctmp(ihpk+1)) )
    ah=tonespec(ihpk)
    ahavp=(sum(tonespec,ismask)-ah)/count(ismask)
    trath=ah/(ahavp+0.01)
    ismask=.false.
    ismask(illo:ilhi)=.true.   ! window for low tone
    iloc=maxloc(tonespec,ismask)
    ilpk=iloc(1)
    deltal=-real( (ctmp(ilpk-1)-ctmp(ilpk+1)) / (2*ctmp(ilpk)-ctmp(ilpk-1)-ctmp(ilpk+1)) )
    al=tonespec(ilpk)
    alavp=(sum(tonespec,ismask)-al)/count(ismask)
    tratl=al/(alavp+0.01)
    fdiff=(ihpk+deltah-ilpk-deltal)*df
    ferrh=(ihpk+deltah-i4000)*df/2.0
    ferrl=(ilpk+deltal-i2000)*df/2.0
    if( ah .ge. al ) then
      ferr=ferrh
    else
      ferr=ferrl
    endif
    detmet(istp)=max(ah,al)
    detmet2(istp)=max(trath,tratl)
    detfer(istp)=ferr
  enddo  ! end of detection-metric and frequency error estimation loop

  call indexx(detmet(1:nstep),nstep,indices) !find median of detection metric vector
  xmed=detmet(indices(nstep/4))
  detmet=detmet/xmed ! noise floor of detection metric is 1.0
  ndet=0

  do ip=1,MAXCAND ! Find candidates
    iloc=maxloc(detmet(1:nstep))
    il=iloc(1)
    if( (detmet(il) .lt. 3.5) ) exit
    if( abs(detfer(il)) .le. ntol ) then
      ndet=ndet+1
      nstart(ndet)=1+(il-1)*60+1
      ferrs(ndet)=detfer(il)
      snrs(ndet)=12.0*log10(detmet(il))/2-9.0
    endif
    detmet(il)=0.0
  enddo

  if( ndet .lt. 3 ) then
    do ip=1,MAXCAND-ndet ! Find candidates
      iloc=maxloc(detmet2(1:nstep))
      il=iloc(1)
      if( (detmet2(il) .lt. 12.0) ) exit
      if( abs(detfer(il)) .le. ntol ) then
        ndet=ndet+1
        nstart(ndet)=1+(il-1)*60+1
        ferrs(ndet)=detfer(il)
        snrs(ndet)=12.0*log10(detmet2(il))/2-9.0
      endif
      detmet2(il)=0.0
    enddo
  endif

  nsuccess=0
  msgreceived=' '
  npeaks=2
  ntol0=29
  deltaf=7.2
  do icand=1,ndet  ! Try to sync/demod/decode each candidate.
    ib=max(1,nstart(icand)-NSPM)
    ie=ib-1+3*NSPM
    if( ie .gt. n ) then
      ie=n
      ib=ie-3*NSPM+1
    endif
    cdat=cbig(ib:ie)
    fo=fc+ferrs(icand)
    xsnr=snrs(icand)
    do iav=1,NPATTERNS
      navmask=navpatterns(1:3,iav)
      call msk40sync(cdat,3,ntol0,deltaf,navmask,npeaks,fo,fest,npkloc,       &
                     nsyncsuccess,c)
      if( nsyncsuccess .eq. 0 ) cycle

      do ipk=1,npeaks
        do is=1,3
          ic0=npkloc(ipk)
          if( is.eq.2) ic0=max(1,ic0-1)
          if( is.eq.3) ic0=min(NSPM,ic0+1)
          ct=cshift(c,ic0-1)
          call msk40decodeframe(ct,mycall,hiscall,xsnr,bswl,nhasharray,        &
                                recent_calls,nrecent,msgreceived,ndecodesuccess)
          if( ndecodesuccess .gt. 0 ) then
!write(*,*) icand, iav, ipk, is, tret, fret, msgreceived
            tret=(nstart(icand)+NSPM/2)/fs
            fret=fest
            navg=sum(navmask)
            nsuccess=ndecodesuccess
            return
          endif
        enddo
      enddo
    enddo
  enddo       ! candidate loop

  return
end subroutine msk40spd


code:
subroutine msk40sync(cdat,nframes,ntol,delf,navmask,npeaks,fc,fest,   &
     npklocs,nsuccess,c)

  !$ use omp_lib

  parameter (NSPM=240)
  complex cdat(NSPM*nframes)
  complex cdat2(NSPM*nframes,8)
  complex c(NSPM)                    !Coherently averaged complex data
  complex cs(NSPM,8)
  complex cb(42)                     !Complex waveform for sync word

  integer s8r(8)
  integer iloc(1)
  integer npklocs(npeaks)
  integer navmask(nframes)                 ! defines which frames to average

  real cbi(42),cbq(42)
  real pkamps(npeaks)
  real xcc(0:NSPM-1)
  real xccs(0:NSPM-1,8)
  real xm(8)
  real bf(8)
  real pp(12)                        !Half-sine pulse shape
  logical first
  data first/.true./
  data s8r/1,0,1,1,0,0,0,1/
  save first,cb,fs,pi,twopi,dt,s8r,pp

  if(first) then
     pi=4.0*atan(1.0)
     twopi=8.0*atan(1.0)
     fs=12000.0
     dt=1.0/fs

     do i=1,12                       !Define half-sine pulse
       angle=(i-1)*pi/12.0
       pp(i)=sin(angle)
     enddo

! Define the sync word waveforms
     s8r=2*s8r-1
     cbq(1:6)=pp(7:12)*s8r(1)
     cbq(7:18)=pp*s8r(3)
     cbq(19:30)=pp*s8r(5)
     cbq(31:42)=pp*s8r(7)
     cbi(1:12)=pp*s8r(2)
     cbi(13:24)=pp*s8r(4)
     cbi(25:36)=pp*s8r(6)
     cbi(37:42)=pp(1:6)*s8r(8)
     cb=cmplx(cbi,cbq)

     first=.false.
  endif

  nfreqs=2*nint(ntol/delf) + 1
  xm=0.0
  bf=0.0
  nthreads=1
  !$ nthreads=min(8,int(OMP_GET_MAX_THREADS(),4))
  nstep=nfreqs/nthreads

  !$OMP PARALLEL NUM_THREADS(nthreads) PRIVATE(id,if1,if2)
  id=1
  !$ id=OMP_GET_THREAD_NUM() + 1            !Thread id = 1,2,...
  if1=-nint(ntol/delf) + (id-1)*nstep
  if2=if1+nstep-1
  if(id.eq.nthreads) if2=nint(ntol/delf)
  call msk40_freq_search(cdat,fc,if1,if2,delf,nframes,navmask,cb,    &
       cdat2(1,id),xm(id),bf(id),cs(1,id),xccs(1,id))
!  write(73,3002) id,if1,if2,nfreqs,nthreads,bf(id),xm(id)
!3002 format(5i5,2f10.3)
  !$OMP END PARALLEL

  xmax=xm(1)
  fest=fc+bf(1)
  c=cs(1:NSPM,1)
  xcc=xccs(0:NSPM-1,1)
  if(nthreads.gt.1) then
     do i=2,nthreads
        if(xm(i).gt.xmax) then
           xmax=xm(i)
           fest=fc+bf(i)
           c=cs(1:NSPM,i)
           xcc=xccs(0:NSPM-1,i)
        endif
     enddo
  endif

! Find npeaks largest peaks
  do ipk=1,npeaks
     iloc=maxloc(xcc)
     ic2=iloc(1)
     npklocs(ipk)=ic2
     pkamps(ipk)=xcc(ic2-1)
     xcc(max(0,ic2-7):min(NSPM-1,ic2+7))=0.0
  enddo

  nsuccess=0
  if( xmax .ge. 1.3 ) nsuccess=1

  return
end subroutine msk40sync


code:
program mskber

! Generate an MSK waveform, pass it through an AWGN channel, apply coherent
! MSK receiver, and count number of errors at each Eb/No.

  parameter (MAXSYM=1000*1000)
  parameter (NSPS=5)                        !Samples per symbol
  real ct(-NSPS:NSPS*MAXSYM-1)              !cos(pi*t/2T)
  real st(-NSPS:NSPS*MAXSYM-1)              !sin(pi*t/2T)
  real r(0:MAXSYM-1)                        !Random numbers to set test bits
  real xsym(0:MAXSYM-1)                     !Soft Rx symbols

  complex xt(-NSPS:NSPS*MAXSYM-1)           !Complex baseband Tx waveform
  complex nt(-NSPS:NSPS*MAXSYM-1)           !Generated AWGN channel noise
  complex yt(-NSPS:NSPS*MAXSYM-1)           !Received signal, yt = xt + fac*nt
  complex cwave(-NSPS:NSPS*MAXSYM-1)        !Audio waveform, Tx real part
  complex z

  integer sym0(0:MAXSYM-1)                  !Generated test bits
  integer sym(0:MAXSYM-1)                   !Hard-copy received bits
  integer sym1(0:7)

  character*12 arg
  data sym1/1,1,0,0,0,1,1,1/

  nargs=iargc()
  if(nargs.ne.2) then
     print*,'Usage:  mskber nsym EbNo'
     go to 999
  endif
  call getarg(1,arg)
  read(arg,*) nsym
  call getarg(2,arg)
  read(arg,*) EbNo

  pi=4.0*atan(1.0)

  do i=-NSPS,NSPS*nsym-1                  !Define ct, st arrays
     t=i*pi/(2.0*NSPS)
     ct(i)=cos(t)
     st(i)=sin(t)
  enddo
  fac=1.0/sqrt(float(NSPS))

  do iEbNo=0,10                                  !Loop over a range of Eb/No
     sym0=0
     call random_number(r)
     where(r(0:nsym-1).gt.0.5) sym0(0:nsym-1)=1  !Generate random data bits
     if(nsym.eq.8) sym0(0:nsym-1)=sym1
     call mskmod(sym0,nsym,NSPS,ct,st,xt,cwave)  !Generate Tx waveform

     do i=-NSPS,NSPS*nsym-1                      !Generate Gaussian noise
        xx=0.707*gran()
        yy=0.707*gran()
        nt(i)=cmplx(xx,yy)
     enddo
     fac_noise=10.0**(-iEbNo/20.0)
     if(EbNo.ne.0.0) fac_noise=10.0**(-EbNo/20.0)
     yt=xt + fac_noise*nt                        !Rx signal, with noise

     call mskdemod(yt,nsym,NSPS,ct,st,xsym)      !MSK demodulator

     sym=0
     where(xsym.gt.0.0) sym=1

     if(nsym.le.160 .and. EbNo.ne.0.0) then
        write(*,1012) sym0(0:nsym-1)
        if(nsym.gt.50) write(*,1012)
        write(*,1012) sym(0:nsym-1)
1012    format(50i1)
        do i=-nsps,nsps*nsym-1
           phi=i*2.0*pi*1500/12000.0
           z=cwave(i)*cmplx(cos(phi),sin(phi))    !Mix back to baseband
           write(51,1014) float(i)/nsps,xt(i),abs(xt(i)),cwave(i),z
1014       format(8f8.4)
        enddo
     endif

! Count the hard errors
     nerr=count(sym(0:nsym-1).ne.sym0(0:nsym-1))
     thber=0.5*erfc(10.0**(iEbNo/20.0))
     xEbNo=iEbNo
     if(EbNo.ne.0.0) xEbNo=EbNo
     write(*,1000) xEbNo,thber,float(nerr)/nsym
1000 format(f6.1,2f10.6)
     if(EbNo.ne.0.0) exit
  enddo

999 end program mskber

subroutine mskmod(sym,nsym,nsps,ct,st,xt,cwave)

! Generate MSK Tx waveform at baseband.

  integer sym(0:nsym-1)                   !Hard-copy received bits
  complex xt(-nsps:nsps*nsym-1)           !Complex baseband Tx waveform
  complex cwave(-nsps:nsps*nsym-1)        !Audio waveform, fc=1500 Hz.
  real ct(-nsps:nsps*nsym-1)              !cos(pi*t/2T)
  real st(-nsps:nsps*nsym-1)              !sin(pi*t/2T)
  real ai(-nsps:nsps*nsym-1)              !Rectangular pulses for even symbols
  real aq(-nsps:nsps*nsym-1)              !Rectangular pulses for odd symbols

  ai=0.
  aq=0.
  fac=1.0/sqrt(float(nsps))
  do j=0,nsym-1,2
     ia=(j-1)*nsps
     ib=ia+2*nsps-1
     ai(ia:ib)=2*sym(j)-1                !Even bits as rectangular pulses
     aq(ia+nsps:ib+nsps)=2*sym(j+1)-1    !Odd bits as rectangular pulses
  enddo
  xt=fac*cmplx(ai*ct,aq*st)                  !Baseband Tx waveform

  twopi=8.0*atan(1.0)
  do i=-nsps,nsps*nsym-1
     phi=i*twopi*1500/12000.0
     cwave(i)=xt(i)*cmplx(cos(phi),-sin(phi))
  enddo

  return
end subroutine mskmod

subroutine mskdemod(yt,nsym,nsps,ct,st,xsym)

! MSK demodulator
! Rx phase must be known and stable; symbol sync must be established.

  complex yt(-nsps:nsps*nsym-1)           !Received signal
  real ct(-nsps:nsps*nsym-1)              !cos(pi*t/2T)
  real st(-nsps:nsps*nsym-1)              !sin(pi*t/2T)
  real xe(-nsps:nsps*nsym-1)              !Temp array for received even symbols
  real xo(-nsps:nsps*nsym-1)              !Temp array for received odd symbols
  real xsym(0:nsym-1)                     !Soft Rx symbols

  iz=nsps*(nsym+1)
  xe(-nsps:nsps*nsym-1)=real(yt)*ct       !Apply matched filters
  xo(-nsps:nsps*nsym-1)=aimag(yt)*st
  do j=0,nsym-1,2
     ia=(j-1)*nsps
     ib=ia+2*nsps-1
     xsym(j)=sum(xe(ia:ib))               !Integrate over 2 symbol lengths
     xsym(j+1)=sum(xo(ia+nsps:ib+nsps))
  enddo

  return
end subroutine mskdemod


code:
subroutine mskdt(d,npts,ty,yellow,nyel)

  parameter (NFFT=1024,NH=NFFT/2)
  real d(npts)
  real x(0:NFFT-1)
  real green(703)
  real yellow(703)                    !703 = 30*12000/512
  real ty(703)
  real y2(175)
  real ty2(175)
  integer indx(703)
  logical ok
  complex c(0:NH)
  equivalence (x,c)

  df=12000.0/NFFT
  i1=nint(300.0/df)
  i2=nint(800.0/df)
  i3=nint(2200.0/df)
  i4=nint(2700.0/df)
  nblks=npts/NH - 1

  do j=1,nblks
     ib=(j+1)*NH
     ia=ib-NFFT+1
     x=d(ia:ib)
     call four2a(x,NFFT,1,-1,0)             !r2c FFT
     sqlow=0.
     do i=i1,i2
        sqlow=sqlow + real(c(i))**2 + aimag(c(i))**2
     enddo
     sqmid=0.
     do i=i2,i3
        sqmid=sqmid + real(c(i))**2 + aimag(c(i))**2
     enddo
     sqhigh=0.
     do i=i3,i4
        sqhigh=sqhigh + real(c(i))**2 + aimag(c(i))**2
     enddo
     green(j)=db(sqlow+sqmid+sqhigh)
     yellow(j)=db(sqmid/(sqlow+sqhigh))
     ty(j)=j*512.0/12000.0
  enddo

  npct=20
  call pctile(green,nblks,npct,base)
  green(1:nblks)=green(1:nblks) - base - 0.3
  call pctile(yellow,nblks,npct,base)
  yellow(1:nblks)=yellow(1:nblks) - base - 0.6
  call indexx(yellow,nblks,indx)

  do j=1,nblks/4
     k=indx(nblks+1-j)
     ty(j)=ty(k)
     yellow(j)=yellow(k)
     if(yellow(j).lt.1.5) exit
  enddo
  nyel=j-1
  k=1
  y2(1)=yellow(1)
  ty2(1)=ty(1)
  do j=2,nyel
     ok=.true.
     do i=1,j-1
        if(abs(ty(i)-ty(j)).lt.0.117) ok=.false.
     enddo
     if(ok) then
        k=k+1
        y2(k)=yellow(j)
        ty2(k)=ty(j)
     endif
  enddo
  nyel=k
  yellow(1:nyel)=y2(1:nyel)
  ty(1:nyel)=ty2(1:nyel)

  return
end subroutine mskdt


code:
subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall,   &
     bshmsg,bcontest,btrain,pcoeffs,bswl,datadir,line)

! Real-time decoder for MSK144.
! Analysis block size = NZ = 7168 samples, t_block = 0.597333 s
! Called from hspec() at half-block increments, about 0.3 s

  parameter (NZ=7168)                !Block size
  parameter (NSPM=864)               !Number of samples per message frame
  parameter (NFFT1=8192)             !FFT size for making analytic signal
  parameter (NPATTERNS=4)            !Number of frame averaging patterns to try
  parameter (NRECENT=10)             !Number of recent calls to remember
  parameter (NSHMEM=50)              !Number of recent SWL messages to remember

  character*4 decsym                 !"&" for mskspd or "^" for long averages
  character*22 msgreceived           !Decoded message
  character*22 msglast,msglastswl   !Used for dupechecking
  character*80 line                  !Formatted line with UTC dB T Freq Msg
  character*12 mycall,hiscall
  character*6 mygrid
  character*12 recent_calls(NRECENT)
  character*22 recent_shmsgs(NSHMEM)
  character*512 datadir

  complex cdat(NFFT1)                !Analytic signal
  complex c(NSPM)                    !Coherently averaged complex data
  complex ct(NSPM)

  integer*2 id2(NZ)                  !Raw 16-bit data
  integer iavmask(8)
  integer iavpatterns(8,NPATTERNS)
  integer npkloc(10)
  integer nhasharray(NRECENT,NRECENT)
  integer nsnrlast,nsnrlastswl

  real d(NFFT1)
  real pow(8)
  real softbits(144)
  real xmc(NPATTERNS)
  real*8 pcoeffs(5)

  logical*1 bshmsg,bcontest,btrain,bswl
  logical*1 first
  logical*1 bshdecode
  logical*1 seenb4
  logical*1 bflag
  logical*1 bvar

  data first/.true./
  data iavpatterns/ &
       1,1,1,1,0,0,0,0, &
       0,0,1,1,1,1,0,0, &
       1,1,1,1,1,0,0,0, &
       1,1,1,1,1,1,1,0/
  data xmc/2.0,4.5,2.5,3.5/     !Used to set time at center of averaging mask
  save first,tsec0,nutc00,pnoise,cdat,msglast,msglastswl,     &
       nsnrlast,nsnrlastswl,recent_calls,nhasharray,recent_shmsgs

  if(first) then
     tsec0=tsec
     nutc00=nutc0
     pnoise=-1.0
     do i=1,nrecent
       recent_calls(i)(1:12)=' '
     enddo
     do i=1,nshmem
       recent_shmsgs(i)(1:22)=' '
     enddo
     msglast='                      '
     msglastswl='                      '
     nsnrlast=-99
     nsnrlastswl=-99
     first=.false.
  endif

  fc=nrxfreq

! Dupe checking setup
  if(nutc00.ne.nutc0 .or. tsec.lt.tsec0) then ! reset dupe checker
    msglast='                      '
    msglastswl='                      '
    nsnrlast=-99
    nsnrlastswl=-99
    nutc00=nutc0
  endif

  tframe=float(NSPM)/12000.0
  line=char(0)
  msgreceived='                      '
  max_iterations=10
  niterations=0
  d(1:NZ)=id2
  rms=sqrt(sum(d(1:NZ)*d(1:NZ))/NZ)
  if(rms.lt.1.0) go to 999
  fac=1.0/rms
  d(1:NZ)=fac*d(1:NZ)
  d(NZ+1:NFFT1)=0.
  bvar=.true.
  if( btrain ) bvar=.false.   ! if training, turn off rx eq
  call analytic(d,NZ,NFFT1,cdat,pcoeffs,bvar)

! Calculate average power for each frame and for the entire block.
! If decode is successful, largest power will be taken as signal+noise.
! If no decode, entire-block average will be used to update noise estimate.
  pmax=-99
  do i=1,8
     ib=(i-1)*NSPM+1
     ie=ib+NSPM-1
     pow(i)=real(dot_product(cdat(ib:ie),cdat(ib:ie)))*rms**2
     pmax=max(pmax,pow(i))
  enddo
  pavg=sum(pow)/8.0

! Short ping decoder uses squared-signal spectrum to determine where to
! center a 3-frame analysis window and attempts to decode each of the
! 3 frames along with 2- and 3-frame averages.
  np=8*NSPM
  call msk144spd(cdat,np,ntol,ndecodesuccess,msgreceived,fc,fest,tdec,navg,ct, &
                 softbits,recent_calls,nrecent)
  if(ndecodesuccess.eq.0 .and. (bshmsg.or.bswl)) then
     call msk40spd(cdat,np,ntol,mycall(1:6),hiscall(1:6),bswl,nhasharray,      &
              recent_calls,nrecent,ndecodesuccess,msgreceived,fc,fest,tdec,navg)
  endif
  if( ndecodesuccess .ge. 1 ) then
    tdec=tsec+tdec
    ipk=0
    is=0
    goto 900
  endif

! If short ping decoder doesn't find a decode,
! Fast - short ping decoder only.
! Normal - try 4-frame averages
! Deep - try 4-, 5- and 7-frame averages.
  npat=NPATTERNS
  if( ndepth .eq. 1 ) npat=0
  if( ndepth .eq. 2 ) npat=2
  do iavg=1,npat
     iavmask=iavpatterns(1:8,iavg)
     navg=sum(iavmask)
     deltaf=10.0/real(navg)  ! search increment for frequency sync
     npeaks=2
     call msk144sync(cdat(1:8*NSPM),8,ntol,deltaf,iavmask,npeaks,fc,           &
          fest,npkloc,nsyncsuccess,xmax,c)
     if( nsyncsuccess .eq. 0 ) cycle

     do ipk=1,npeaks
        do is=1,3
           ic0=npkloc(ipk)
           if(is.eq.2) ic0=max(1,ic0-1)
           if(is.eq.3) ic0=min(NSPM,ic0+1)
           ct=cshift(c,ic0-1)
           call msk144decodeframe(ct,softbits,msgreceived,ndecodesuccess,      &
                                  recent_calls,nrecent)
           if(ndecodesuccess .gt. 0) then
              tdec=tsec+xmc(iavg)*tframe
              goto 900
           endif
        enddo                         !Slicer dither
     enddo                            !Peak loop
  enddo

  msgreceived=' '

! no decode - update noise level used for calculating displayed snr.
  if( pnoise .lt. 0 ) then         ! initialize noise level
     pnoise=pavg
  elseif( pavg .gt. pnoise ) then  ! noise level is slow to rise
     pnoise=0.9*pnoise+0.1*pavg
  elseif( pavg .lt. pnoise ) then  ! and quick to fall
     pnoise=pavg
  endif
  go to 999

900 continue
! Successful decode - estimate snr
  if( pnoise .gt. 0.0 ) then
    snr0=10.0*log10(pmax/pnoise-1.0)
  else
    snr0=0.0
  endif
  nsnr=nint(snr0)

  bshdecode=.false.
  if( msgreceived(1:1) .eq. '<' ) bshdecode=.true.

  if(.not. bshdecode) then
    call msk144signalquality(ct,snr0,fest,tdec,softbits,msgreceived,hiscall,   &
                          btrain,datadir,ncorrected,eyeopening,pcoeffs)
  endif

  decsym=' &  '
  if( btrain ) decsym=' ^  '
  if( msgreceived(1:1).eq.'<') then
    ncorrected=0
    eyeopening=0.0
  endif

  if( nsnr .lt. -8 ) nsnr=-8
  if( nsnr .gt. 24 ) nsnr=24

! Dupe check.
  bflag=ndecodesuccess.eq.1 .and.                                              &
        (msgreceived.ne.msglast .or. nsnr.gt.nsnrlast .or. tsec.lt.tsec0)
  if(bflag) then
     msglast=msgreceived
     nsnrlast=nsnr
     if(.not. bshdecode) then
        call update_hasharray(recent_calls,nrecent,nhasharray)
        if(bcontest) call fix_contest_msg(mygrid,msgreceived)
     endif
     write(line,1020) nutc0,nsnr,tdec,nint(fest),decsym,msgreceived,           &
          navg,ncorrected,eyeopening,char(0)
1020 format(i6.6,i4,f5.1,i5,a4,a22,i2,i3,f5.1,a1)
  elseif(bswl .and. ndecodesuccess.ge.2) then
    seenb4=.false.
    do i=1,nshmem
      if( msgreceived .eq. recent_shmsgs(i) ) then
        seenb4=.true.
      endif
    enddo
    call update_recent_shmsgs(msgreceived,recent_shmsgs,nshmem)
    bflag=seenb4 .and.                                                        &
      (msgreceived.ne.msglastswl .or. nsnr.gt.nsnrlastswl .or. tsec.lt.tsec0) &
      .and. nsnr.gt.-6
    if(bflag) then
      msglastswl=msgreceived
      nsnrlastswl=nsnr
      write(line,1020) nutc0,nsnr,tdec,nint(fest),decsym,msgreceived,    &
          navg,ncorrected,eyeopening,char(0)
    endif
  endif
999 tsec0=tsec

  return
end subroutine mskrtd

subroutine update_recent_shmsgs(message,msgs,nsize)
  character*22 msgs(nsize)
  character*22 message
  logical*1 seen

  seen=.false.
  do i=1,nsize
    if( msgs(i) .eq. message ) seen=.true.
  enddo

  if( .not. seen ) then
    do i=nsize,2,-1
      msgs(i)=msgs(i-1)
    enddo
    msgs(1)=message
  endif

  return
end subroutine update_recent_shmsgs


code:
subroutine noisegen(d4,nmax)

  real*4 d4(4,nmax)

  call init_random_seed()       ! seed Fortran RANDOM_NUMBER generator
  call sgran()                  ! see C rand generator (used in gran)

  do i=1,nmax
     d4(1,i)=gran()
     d4(2,i)=gran()
     d4(3,i)=gran()
     d4(4,i)=gran()
  enddo

  return
end subroutine noisegen