code:
program dbpsksim

  parameter (ND=121)                    !Data symbols: LDPC (120,60), r=1/2
  parameter (NN=ND)                     !Total symbols (121)
  parameter (NSPS=28800)                !Samples per symbol at 12000 sps
  parameter (NZ=NSPS*NN)                !Samples in waveform (3484800)
  parameter (NFFT1=65536,NH1=NFFT1/2)
  parameter (NFFT2=128,NH2=NFFT2/2)

  character*8 arg
  complex c(0:NZ-1)                     !Complex waveform
  complex c2(0:NFFT1-1)                 !Short spectra
  complex cr(0:NZ-1)
  complex ct(0:NZ-1)
  complex cz(0:NFFT2-1)
  complex z0,z,zp
  real s(-NH1+1:NH1)
  real s2(-NH2+1:NH2)
  real xnoise(0:NZ-1)                   !Generated random noise
  real ynoise(0:NZ-1)                   !Generated random noise
  real rxdata(120),llr(120)
  integer id(NN)                        !Encoded NRZ data (values +/-1)
  integer id1(NN)                       !Recovered data (1st pass)
  integer id2(NN)                       !Recovered data (2nd pass)
!  integer icw(NN)
  integer*1 msgbits(60),decoded(60),codeword(120),apmask(120),cw(120)
  data msgbits/0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,&
       0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,1,0,0,1,0,1,1,0,1,0/

  nnn=0
  nargs=iargc()
  if(nargs.ne.6) then
     print*,'Usage:   dbpsksim f0(Hz) delay(ms) fspread(Hz) ndiff iters snr(dB)'
     print*,'Example: dbpsksim 1500 0 0 10 -35'
     print*,'Set snr=0 to cycle through a range'
     go to 999
  endif
  call getarg(1,arg)
  read(arg,*) f0                        !Low tone frequency
  call getarg(2,arg)
  read(arg,*) delay
  call getarg(3,arg)
  read(arg,*) fspread
  call getarg(4,arg)
  read(arg,*) ndiff
  call getarg(5,arg)
  read(arg,*) iters
  call getarg(6,arg)
  read(arg,*) snrdb

  twopi=8.d0*atan(1.d0)
  fs=12000.d0
  dt=1.0/fs
  ts=NSPS*dt
  baud=1.d0/ts
  txt=NZ*dt
  bandwidth_ratio=2500.0/6000.0
  write(*,1000) baud,5*baud,txt,delay,fspread,ndiff
1000 format('Baud:',f6.3,'  BW:',f4.1,'  TxT:',f6.1,'  Delay:',f5.2,   &
          '  fSpread:',f5.2,'  ndiff:',i2/)

  write(*,1004)
1004 format('  SNR    err   ber    fer   fsigma'/35('-'))

  call encode120(msgbits,codeword)                !Encode the test message
  isna=-28
  isnb=-40
  if(snrdb.ne.0.0) then
     isna=nint(snrdb)
     isnb=isna
  endif
  do isnr=isna,isnb,-1
     snrdb=isnr
     sig=sqrt(bandwidth_ratio) * 10.0**(0.05*snrdb)
     if(snrdb.gt.90.0) sig=1.0
     nhard=0
     nhardc=0
     nfe1=0
     nfe2=0
     sqf=0.
     do iter=1,iters
        nnn=nnn+1
        id(1)=1                                  !First bit is always 1
        id(2:NN)=2*codeword-1
        call genbpsk(id,f0,ndiff,0,c)            !Generate the 4-FSK waveform
        if(delay.ne.0.0 .or. fspread.ne.0.0) call watterson(c,delay,fspread)
        c=sig*c                                  !Scale to requested SNR
        if(snrdb.lt.90) then
           do i=0,NZ-1                           !Generate gaussian noise
              xnoise(i)=gran()
              ynoise(i)=gran()
           enddo
           c=c + cmplx(xnoise,ynoise)            !Add noise to signal
        endif

! First attempt at finding carrier frequency fc: 64k FFTs ==> avg power spectra
        nspec=NZ/NFFT1
        df1=12000.0/NFFT1
        s=0.
        do k=1,nspec
           ia=(k-1)*NSPS
           ib=ia+NSPS-1
           c2(0:NSPS-1)=c(ia:ib)
           c2(NSPS:)=0.
           call four2a(c2,NFFT1,1,-1,1)
           do i=0,NFFT1-1
              j=i
              if(j.gt.NH1) j=j-NFFT1
              s(j)=s(j) + real(c2(i))**2 + aimag(c2(i))**2
           enddo
        enddo
        s=1.e-6*s
        smax=0.
        ipk=0
        ia=(1400.0)/df1
        ib=(1600.0)/df1
        do i=ia,ib
           f=i*df1
           if(s(i).gt.smax) then
              smax=s(i)
              ipk=i
              fc=f
           endif
        enddo
        a=(s(ipk+1)-s(ipk-1))/2.0
        b=(s(ipk+1)+s(ipk-1)-2.0*s(ipk))/2.0
        dx=-a/(2.0*b)
        fc=fc + df1*dx                           !Estimated carrier frequency
        sqf=sqf + (fc-f0)**2

! The following is for testing SNR calibration:
!        sp5n=(s(ipk-2)+s(ipk-1)+s(ipk)+s(ipk+1)+s(ipk+2))  !Sig + 5*noise
!        base=(sum(s)-sp5n)/(NFFT1-5.0)                     !Noise per bin
!        psig=sp5n-5*base                                   !Sig only
!        pnoise=(2500.0/df1)*base                           !Noise in 2500 Hz
!        xsnrdb=db(psig/pnoise)

        call genbpsk(id,fc,ndiff,1,cr)           !Generate reference carrier
        c=c*conjg(cr)                            !Mix signal to baseband

        z0=1.0
        do j=1,NN                                !Demodulate
           ia=(j-1)*NSPS
           ib=ia+NSPS-1
           z=sum(c(ia:ib))
           cz(j-1)=z
           zp=z*conjg(z0)
           p=1.e-4*real(zp)
           id1(j)=-1
           if(p.ge.0.0) id1(j)=1
           if(j.ge.2) rxdata(j-1)=p
           z0=z
        enddo

        rxav=sum(rxdata)/120
        rx2av=sum(rxdata*rxdata)/120
        rxsig=sqrt(rx2av-rxav*rxav)
        rxdata=rxdata/rxsig
        ss=0.84
        llr=2.0*rxdata/(ss*ss)
        apmask=0
        max_iterations=10
        call bpdecode120(llr,apmask,max_iterations,decoded,niterations,cw)

! Count frame errors
        if(niterations.lt.0 .or. count(msgbits.ne.decoded).gt.0) nfe1=nfe1+1

! Find carrier frequency from squared cz array.
        cz(121:)=0.
        cz=cz*cz
        call four2a(cz,NFFT2,1,-1,1)
        s2max=0.
        do i=0,NFFT2-1
           j=i
           if(i.gt.NH2) j=j-NFFT2
           s2(j)=real(cz(i))**2 + aimag(cz(i))**2
           if(s2(j).gt.s2max) then
              s2max=s2(j)
              jpk=j
           endif
!           write(16,1200) j*baud/NFFT2,1.e-12*s2(j)
!1200       format(2f12.3)
        enddo
        a=(s2(jpk+1)-s2(jpk-1))/2.0
        b=(s2(jpk+1)+s2(jpk-1)-2.0*s2(jpk))/2.0
        dx=-a/(2.0*b)
        fc2=0.5*(jpk+dx)*baud/NFFT2

        call genbpsk(id,fc2,ndiff,1,cr)         !Generate new ref carrier at fc2
        c=c*conjg(cr)
        z0=1.0
        do j=1,NN                                !Demodulate
           ia=(j-1)*NSPS
           ib=ia+NSPS-1
           z=sum(c(ia:ib))
           if(j.eq.1) z0=z
           zp=z*conjg(z0)
           p=1.e-4*real(zp)
           id2(j)=-1
           if(p.ge.0.0) id2(j)=1
           if(j.ge.2) rxdata(j-1)=p
           ierr=0
           if(id2(j).ne.id(j)) ierr=1
           id3=-1
           if(real(z).ge.0.0) id3=1
           if(j.ge.2 .and. id3.ne.id(j)) nhardc=nhardc+1
           if(j.ge.2 .and. ndiff.eq.0) rxdata(j-1)=real(z)
           z0=z
        enddo
        nhard=nhard + count(id2.ne.id)           !Count hard errors

        rxav=sum(rxdata)/120
        rx2av=sum(rxdata*rxdata)/120
        rxsig=sqrt(rx2av-rxav*rxav)
        rxdata=rxdata/rxsig
        ss=0.84
        llr=2.0*rxdata/(ss*ss)               !Soft symbols
        apmask=0
        max_iterations=10
        decoded=0
        call bpdecode120(llr,apmask,max_iterations,decoded,niterations,cw)
!        if(niterations.lt.0) then
!           llr=-llr
!           call bpdecode120(llr,apmask,max_iterations,decoded,niterations,cw)
!           if(niterations.ge.0) nhard=NN*iters-nhard
!        endif
        if(niterations.ge.0) call chkcrc10(decoded,nbadcrc)
        if(niterations.lt.0 .or. count(msgbits.ne.decoded).gt.0 .or.        &
             nbadcrc.ne.0) nfe2=nfe2+1
     enddo

     if(ndiff.eq.0) nhard=nhardc
     fsigma=sqrt(sqf/iters)
     ber=float(nhard)/(NN*iters)
     fer=float(nfe2)/iters
     write(*,1050)  snrdb,nhard,ber,fer,fsigma
     write(14,1050)  snrdb,nhard,ber,fer,fsigma
1050 format(f6.1,i5,f8.4,f7.3,f8.2)
  enddo

999 end program dbpsksim


code:
subroutine dopspread(c,fspread)

  parameter (NFFT=268800,NH=NFFT/2)
  complex c(0:NFFT-1)
  complex cspread(0:NFFT-1)

  df=12000.0/nfft
  twopi=8*atan(1.0)
  cspread(0)=1.0
  cspread(NH)=0.
  b=6.0                                     !Lorenzian 3/28 onward
  do i=1,NH
     f=i*df
     x=b*f/fspread
     z=0.
     a=0.
     if(x.lt.3.0) then                          !Cutoff beyond x=3
        a=sqrt(1.111/(1.0+x*x)-0.1)             !Lorentzian
        call random_number(r1)
        phi1=twopi*r1
        z=a*cmplx(cos(phi1),sin(phi1))
     endif
     cspread(i)=z
     z=0.
     if(x.lt.50.0) then
        call random_number(r2)
        phi2=twopi*r2
        z=a*cmplx(cos(phi2),sin(phi2))
     endif
     cspread(NFFT-i)=z
  enddo

  izh=fspread/df
  do i=-izh,izh
     f=i*df
     j=i
     if(j.lt.0) j=j+nfft
     s=real(cspread(j))**2 + aimag(cspread(j))**2
!          write(23,3000) f,s,cspread(j)
!3000      format(f10.3,3f12.6)
  enddo

  call four2a(cspread,NFFT,1,1,1)             !Transform to time domain

  sum=0.
  do i=0,NFFT-1
     p=real(cspread(i))**2 + aimag(cspread(i))**2
     sum=sum+p
  enddo
  avep=sum/NFFT
  fac=sqrt(1.0/avep)
  cspread=fac*cspread                   !Normalize to constant avg power
  c=cspread*c                           !Apply Rayleigh fading to c()

  do i=0,NFFT-1
     p=real(cspread(i))**2 + aimag(cspread(i))**2
!     write(24,3010) i,p,cspread(i)
!3010 format(i8,3f12.6)
  enddo

  return
end subroutine dopspread

license: mit


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

include "ldpc_204_68_params.f90"

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

save first,gen

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

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

return
end subroutine encode204


code:
subroutine extractmessage168(decoded,msgreceived,ncrcflag,recent_calls,nrecent)
  use iso_c_binding, only: c_loc,c_size_t
  use crc
  use packjt

  character*22 msgreceived
  character*12 call1,call2
  character*12 recent_calls(nrecent)
  integer*1 decoded(84)
  integer*1, target::  i1Dec8BitBytes(11)
  integer*4 i4Dec6BitWords(12)

! Collapse 84 decoded bits to 11 bytes. Bytes 1-9 are the message, byte 10 and first half of byte 11 is the crc
  do ibyte=1,9
      itmp=0
    do ibit=1,8
      itmp=ishft(itmp,1)+iand(1,decoded((ibyte-1)*8+ibit))
    enddo
    i1Dec8BitBytes(ibyte)=itmp
  enddo
! Need to pack the crc into bytes 10 and 11 for crc12_check
  i1Dec8BitBytes(10)=decoded(73)*8+decoded(74)*4+decoded(75)*2+decoded(76)
  i1Dec8BitBytes(11)=decoded(77)*128+decoded(78)*64+decoded(79)*2*32+decoded(80)*16
  i1Dec8BitBytes(11)=i1Dec8BitBytes(11)+decoded(81)*8+decoded(82)*4+decoded(83)*2+decoded(84)

  if( crc12_check(c_loc (i1Dec8BitBytes), 11) ) then
! CRC12 checks out --- unpack 72-bit message
    do ibyte=1,12
      itmp=0
      do ibit=1,6
        itmp=ishft(itmp,1)+iand(1,decoded((ibyte-1)*6+ibit))
      enddo
      i4Dec6BitWords(ibyte)=itmp
    enddo
    call unpackmsg144(i4Dec6BitWords,msgreceived,call1,call2)
    ncrcflag=1
    if( call1(1:2) .ne. 'CQ' .and. call1(1:2) .ne. '  ' ) then
      call update_recent_calls(call1,recent_calls,nrecent)
    endif
    if( call2(1:2) .ne. '  ' ) then
      call update_recent_calls(call2,recent_calls,nrecent)
    endif
  else
    msgreceived=' '
    ncrcflag=-1
  endif
  return
  end subroutine extractmessage168


code:
  INTEGER FFTW_R2HC
  PARAMETER (FFTW_R2HC=0)
  INTEGER FFTW_HC2R
  PARAMETER (FFTW_HC2R=1)
  INTEGER FFTW_DHT
  PARAMETER (FFTW_DHT=2)
  INTEGER FFTW_REDFT00
  PARAMETER (FFTW_REDFT00=3)
  INTEGER FFTW_REDFT01
  PARAMETER (FFTW_REDFT01=4)
  INTEGER FFTW_REDFT10
  PARAMETER (FFTW_REDFT10=5)
  INTEGER FFTW_REDFT11
  PARAMETER (FFTW_REDFT11=6)
  INTEGER FFTW_RODFT00
  PARAMETER (FFTW_RODFT00=7)
  INTEGER FFTW_RODFT01
  PARAMETER (FFTW_RODFT01=8)
  INTEGER FFTW_RODFT10
  PARAMETER (FFTW_RODFT10=9)
  INTEGER FFTW_RODFT11
  PARAMETER (FFTW_RODFT11=10)
  INTEGER FFTW_FORWARD
  PARAMETER (FFTW_FORWARD=-1)
  INTEGER FFTW_BACKWARD
  PARAMETER (FFTW_BACKWARD=+1)
  INTEGER FFTW_MEASURE
  PARAMETER (FFTW_MEASURE=0)
  INTEGER FFTW_DESTROY_INPUT
  PARAMETER (FFTW_DESTROY_INPUT=1)
  INTEGER FFTW_UNALIGNED
  PARAMETER (FFTW_UNALIGNED=2)
  INTEGER FFTW_CONSERVE_MEMORY
  PARAMETER (FFTW_CONSERVE_MEMORY=4)
  INTEGER FFTW_EXHAUSTIVE
  PARAMETER (FFTW_EXHAUSTIVE=8)
  INTEGER FFTW_PRESERVE_INPUT
  PARAMETER (FFTW_PRESERVE_INPUT=16)
  INTEGER FFTW_PATIENT
  PARAMETER (FFTW_PATIENT=32)
  INTEGER FFTW_ESTIMATE
  PARAMETER (FFTW_ESTIMATE=64)
  INTEGER FFTW_ESTIMATE_PATIENT
  PARAMETER (FFTW_ESTIMATE_PATIENT=128)
  INTEGER FFTW_BELIEVE_PCOST
  PARAMETER (FFTW_BELIEVE_PCOST=256)
  INTEGER FFTW_DFT_R2HC_ICKY
  PARAMETER (FFTW_DFT_R2HC_ICKY=512)
  INTEGER FFTW_NONTHREADED_ICKY
  PARAMETER (FFTW_NONTHREADED_ICKY=1024)
  INTEGER FFTW_NO_BUFFERING
  PARAMETER (FFTW_NO_BUFFERING=2048)
  INTEGER FFTW_NO_INDIRECT_OP
  PARAMETER (FFTW_NO_INDIRECT_OP=4096)
  INTEGER FFTW_ALLOW_LARGE_GENERIC
  PARAMETER (FFTW_ALLOW_LARGE_GENERIC=8192)
  INTEGER FFTW_NO_RANK_SPLITS
  PARAMETER (FFTW_NO_RANK_SPLITS=16384)
  INTEGER FFTW_NO_VRANK_SPLITS
  PARAMETER (FFTW_NO_VRANK_SPLITS=32768)
  INTEGER FFTW_NO_VRECURSE
  PARAMETER (FFTW_NO_VRECURSE=65536)
  INTEGER FFTW_NO_SIMD
  PARAMETER (FFTW_NO_SIMD=131072)


code:
subroutine four2a(a,nfft,ndim,isign,iform)

! IFORM = 1, 0 or -1, as data is
! complex, real, or the first half of a complex array.  Transform
! values are returned in array DATA.  They are complex, real, or
! the first half of a complex array, as IFORM = 1, -1 or 0.

! The transform of a real array (IFORM = 0) dimensioned N(1) by N(2)
! by ... will be returned in the same array, now considered to
! be complex of dimensions N(1)/2+1 by N(2) by ....  Note that if
! IFORM = 0 or -1, N(1) must be even, and enough room must be
! reserved.  The missing values may be obtained by complex conjugation.

! The reverse transformation of a half complex array dimensioned
! N(1)/2+1 by N(2) by ..., is accomplished by setting IFORM
! to -1.  In the N array, N(1) must be the true N(1), not N(1)/2+1.
! The transform will be real and returned to the input array.

! This version of four2a makes calls to the FFTW library to do the
! actual computations.

  parameter (NPMAX=2100)                 !Max numberf of stored plans
  parameter (NSMALL=16384)               !Max size of "small" FFTs
  complex a(nfft)                        !Array to be transformed
  complex aa(NSMALL)                     !Local copy of "small" a()
  integer nn(NPMAX),ns(NPMAX),nf(NPMAX)  !Params of stored plans
  integer*8 nl(NPMAX),nloc               !More params of plans
  integer*8 plan(NPMAX)                  !Pointers to stored plans
  logical found_plan
  data nplan/0/                          !Number of stored plans
  common/patience/npatience,nthreads     !Patience and threads for FFTW plans
  include 'fftw3.f90'                    !FFTW definitions
  save plan,nplan,nn,ns,nf,nl

  if(nfft.lt.0) go to 999

  nloc=loc(a)

  found_plan = .false.
  !$omp critical(four2a_setup)
  do i=1,nplan
     if(nfft.eq.nn(i) .and. isign.eq.ns(i) .and.                     &
          iform.eq.nf(i) .and. nloc.eq.nl(i)) then
        found_plan = .true.
        exit
     end if
  enddo

  if(i.ge.NPMAX) stop 'Too many FFTW plans requested.'

  if (.not. found_plan) then
     nplan=nplan+1
     i=nplan

     nn(i)=nfft
     ns(i)=isign
     nf(i)=iform
     nl(i)=nloc

! Planning: FFTW_ESTIMATE, FFTW_ESTIMATE_PATIENT, FFTW_MEASURE,
!            FFTW_PATIENT,  FFTW_EXHAUSTIVE
     nflags=FFTW_ESTIMATE
     if(npatience.eq.1) nflags=FFTW_ESTIMATE_PATIENT
     if(npatience.eq.2) nflags=FFTW_MEASURE
     if(npatience.eq.3) nflags=FFTW_PATIENT
     if(npatience.eq.4) nflags=FFTW_EXHAUSTIVE

     if(nfft.le.NSMALL) then
        jz=nfft
        if(iform.eq.0) jz=nfft/2
        aa(1:jz)=a(1:jz)
     endif

     !$omp critical(fftw) ! serialize non thread-safe FFTW3 calls
     if(isign.eq.-1 .and. iform.eq.1) then
        call sfftw_plan_dft_1d(plan(i),nfft,a,a,FFTW_FORWARD,nflags)
     else if(isign.eq.1 .and. iform.eq.1) then
        call sfftw_plan_dft_1d(plan(i),nfft,a,a,FFTW_BACKWARD,nflags)
     else if(isign.eq.-1 .and. iform.eq.0) then
        call sfftw_plan_dft_r2c_1d(plan(i),nfft,a,a,nflags)
     else if(isign.eq.1 .and. iform.eq.-1) then
        call sfftw_plan_dft_c2r_1d(plan(i),nfft,a,a,nflags)
     else
        stop 'Unsupported request in four2a'
     endif
     !$omp end critical(fftw)

     if(nfft.le.NSMALL) then
        jz=nfft
        if(iform.eq.0) jz=nfft/2
        a(1:jz)=aa(1:jz)
     endif
  end if
  !$omp end critical(four2a_setup)

  call sfftw_execute(plan(i))
  return

999 continue

  !$omp critical(four2a)
  do i=1,nplan
! The test is only to silence a compiler warning:
     if(ndim.ne.-999) then
        !$omp critical(fftw) ! serialize non thread-safe FFTW3 calls
        call sfftw_destroy_plan(plan(i))
        !$omp end critical(fftw)
     end if
  enddo

  nplan=0
  !$omp end critical(four2a)

  return
end subroutine four2a


code:
program fsk4hf

! Simulate characteristics of a potential mode using LDPC (168,84) code,
! 4-FSK modulation, and 30 s T/R sequences.

  parameter (KK=84)                     !Information bits (72 + CRC12)
  parameter (ND=84)                     !Data symbols: LDPC (168,84), r=1/2
  parameter (NS=12)                     !Sync symbols (3 @ 4x4 Costas arrays)
  parameter (NR=2)                      !Ramp up/down
  parameter (NN=NR+NS+ND)               !Total symbols (98)
  parameter (NSPS=2688/84)              !Samples per symbol (32)
  parameter (NZ=NSPS*NN)                !Samples in baseband waveform (3760)

  character*8 arg
  complex c0(0:NZ-1)                    !Complex waveform
  complex c(0:NZ-1)                     !Complex waveform
  real xnoise(0:NZ-1)                   !Generated random noise
  real ynoise(0:NZ-1)                   !Generated random noise
  real rxdata(2*ND),llr(2*ND)           !Soft symbols
  real s(0:NSPS,NN)
  real savg(0:NSPS)
  real ps(0:3)
  integer id(ND)                        !Symbol values (0-3), data only
  integer id1(ND)                       !Recovered data values
  integer*1 msgbits(KK),decoded(KK),apmask(ND),cw(ND)
  data msgbits/0,0,1,0,0,1,1,1,1,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,1,1,0,0,0,1, &
       1,1,1,0,1,1,1,1,1,1,1,0,0,1,0,0,1,1,0,1,0,1,1,1,0,1,1,0,1,1,         &
       1,1,0,1,0,1,1,0,0,0,0,0,1,0,0,0,0,0,1,0,1,0,1,0/

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

  twopi=8.0*atan(1.0)
  fs=12000.0/84.0                        !Sample rate = 142.857... Hz
  dt=1.0/fs                              !Sample interval (s)
  tt=NSPS*dt                             !Duration of "itone" symbols (s)
  baud=1.0/tt                            !Keying rate for "itone" symbols (baud)
  txt=NZ*dt                              !Transmission length (s)
  bandwidth_ratio=2500.0/(fs/2.0)
  write(*,1000) f0,delay,fspread,iters,baud,4*baud,txt
1000 format('f0:',f5.1,'  Delay:',f4.1,'  fSpread:',f5.2,                    &
          '  Iters:',i6/'Baud:',f7.3,'  BW:',f5.1,'  TxT:',f5.1,f5.2/)
  write(*,1004)
1004 format(/'  SNR    sym   bit    ser     ber    fer   fsigma'/50('-'))

  call genfsk4hf(msgbits,f0,id,c0)       !Generate baseband waveform
  isna=-10
  isnb=-30
  if(snrdb.ne.0.0) then
     isna=nint(snrdb)
     isnb=isna
  endif
  do isnr=isna,isnb,-1                   !Loop over SNR range
     snrdb=isnr
     sig=sqrt(bandwidth_ratio) * 10.0**(0.05*snrdb)
     if(snrdb.gt.90.0) sig=1.0
     nhard=0
     nbit=0
     nfe=0
     sqf=0.
     do iter=1,iters                     !Loop over requested iterations
        c=c0
        if(delay.ne.0.0 .or. fspread.ne.0.0) then
           call watterson(c,NZ,fs,delay,fspread)
        endif
        c=sig*c                          !Scale to requested SNR
        if(snrdb.lt.90) then
           do i=0,NZ-1                   !Generate gaussian noise
              xnoise(i)=gran()
              ynoise(i)=gran()
           enddo
           c=c + cmplx(xnoise,ynoise)    !Add AWGN noise
        endif
        df=fs/(2*NSPS)
        i0=nint(f0/df)
        call spec4(c,s,savg)
        do i=0,NSPS
           write(12,3001) i*df,savg(i),db(savg(i))
3001       format(3f15.3)
        enddo

        do j=1,ND
           nlo=0
           nhi=0
           k=j+5
           if(j.ge.43) k=j+9
           ps=s(i0:i0+6:2,k)
           ps=sqrt(ps)  !###
           rlo=max(ps(1),ps(3))-max(ps(0),ps(2))
           rhi=max(ps(2),ps(3))-max(ps(0),ps(1))
           if(rlo.ge.0.0) nlo=1
           if(rhi.ge.0.0) nhi=1
           rxdata(2*j-1)=rhi
           rxdata(2*j)=rlo
           id1(j)=2*nhi+nlo
        enddo
!        write(*,1001) id(1:70)
!        write(*,1001) id1(1:70)
!1001    format(70i1)
        nhard=nhard+count(id.ne.id1)
        nbit=nbit + count(iand(id,1).ne.iand(id1,1)) +               &
             count(iand(id,2).ne.iand(id1,2))

        rxav=sum(rxdata)/ND
        rx2av=sum(rxdata*rxdata)/ND
        rxsig=sqrt(rx2av-rxav*rxav)
        rxdata=rxdata/rxsig
        ss=0.84
        llr=2.0*rxdata/(ss*ss)
        apmask=0
        max_iterations=40
        ifer=0
        call bpdecode168(llr,apmask,max_iterations,decoded,niterations,cw)
        nbadcrc=0
        if(niterations.ge.0) call chkcrc12(decoded,nbadcrc)
        if(niterations.lt.0 .or. count(msgbits.ne.decoded).gt.0 .or.        &
             nbadcrc.ne.0) ifer=1
        nfe=nfe+ifer
     enddo

     fsigma=sqrt(sqf/iters)
     ser=float(nhard)/(ND*iters)
     ber=float(nbit)/(2*ND*iters)
     fer=float(nfe)/iters
     write(*,1050)  snrdb,nhard,nbit,ser,ber,fer,fsigma
!     write(60,1050)  snrdb,nhard,ber,fer,fsigma
1050 format(f6.1,2i6,2f8.4,f7.3,f8.2)
  enddo

999 end program fsk4hf


code:
program fsk4sim

  parameter (ND=60)                      !Data symbols: LDPC (120,60), r=1/2
  parameter (NN=ND)                      !Total symbols (60)
  parameter (NSPS=57600)                 !Samples per symbol at 12000 sps
  parameter (NZ=NSPS*NN)                 !Samples in waveform (3456000)

  character*8 arg
  complex c(0:NZ-1)                      !Complex waveform
  complex cr(0:NZ-1)
  complex cs(NSPS,NN)
  complex cps(0:3)
  complex ct(0:2*NN-1)
  complex z,w,zsum
  real r(0:NZ-1)
  real s(NSPS,NN)
  real savg(NSPS)
  real tmp(NN)                          !For generating random data
  real xnoise(0:NZ-1)                   !Generated random noise
  real ps(0:3)
  integer id(NN)                         !Encoded 2-bit data (values 0-3)
  integer id2(NN)                        !Recovered data
  equivalence (r,cr)

  nnn=0
  nargs=iargc()
  if(nargs.ne.6) then
     print*,'Usage: fsk8sim f0 delay(ms) fspread(Hz) nts iters snr(dB)'
     go to 999
  endif
  call getarg(1,arg)
  read(arg,*) f0                        !Low tone frequency
  call getarg(2,arg)
  read(arg,*) delay
  call getarg(3,arg)
  read(arg,*) fspread
  call getarg(4,arg)
  read(arg,*) nts
  call getarg(5,arg)
  read(arg,*) iters
  call getarg(6,arg)
  read(arg,*) snrdb

  twopi=8.d0*atan(1.d0)
  fs=12000.d0
  dt=1.0/fs
  ts=NSPS*dt
  baud=1.d0/ts
  txt=NZ*dt
  bandwidth_ratio=2500.0/6000.0
  write(*,1000) baud,5*baud,txt,delay,fspread,nts
1000 format('Baud:',f6.3,'  BW:',f5.1,'  TxT:',f5.1,'  Delay:',f5.2,   &
          '  fSpread:',f5.2,'  nts:',i3/)

  write(*,1004)
1004 format('  SNR    Sym  Bit   SER     BER    Sym  Bit   SER     BER'/59('-'))

  isna=-25
  isnb=-40
  if(snrdb.ne.0.0) then
     isna=nint(snrdb)
     isnb=isna
  endif
  do isnr=isna,isnb,-1
     snrdb=isnr
     sig=sqrt(2*bandwidth_ratio) * 10.0**(0.05*snrdb)
     if(snrdb.gt.90.0) sig=1.0
     nhard1=0
     nhard2=0
     nbit1=0
     nbit2=0
     nh2=0
     nb2=0
     do iter=1,iters
        nnn=nnn+1
        id=0
        call random_number(tmp)
        where(tmp.ge.0.25 .and. tmp.lt.0.50) id=1
        where(tmp.ge.0.50 .and. tmp.lt.0.75) id=2
        where(tmp.ge.0.75) id=3

        call genfsk4(id,f0,nts,c)                !Generate the 4-FSK waveform
        call watterson(c,delay,fspread)
        if(sig.ne.1.0) c=sig*c                   !Scale to requested SNR
        if(snrdb.lt.90) then
           do i=0,NZ-1                           !Generate gaussian noise
              xnoise(i)=gran()
           enddo
        endif
        r(0:NZ-1)=real(c(0:NZ-1)) + xnoise       !Add noise to signal

        call snr2_wsprlf(r,freq,snr2500,width,1)
        write(*,3001) freq,snr2500,width
3001    format(40x,3f10.3)

        df=12000.0/(2*NSPS)
!        i0=nint(f0/df)
!        i0=nint((1500.0+freq)/df)
        i0=nint((f0+freq)/df)
        call spec4(r,cs,s,savg)

        do j=1,NN
           nlo=0
           nhi=0
           ps=s(i0:i0+6*nts:2*nts,j)
           cps=cs(i0:i0+6*nts:2*nts,j)
           if(max(ps(1),ps(3)).ge.max(ps(0),ps(2))) nlo=1
           if(max(ps(2),ps(3)).ge.max(ps(0),ps(1))) nhi=1
           id2(j)=2*nhi+nlo
           z=cps(id2(j))
           ct(j-1)=z
        enddo
        nh1=count(id.ne.id2)
        nb1=count(iand(id,1).ne.iand(id2,1)) + count(iand(id,2).ne.iand(id2,2))

        ct(NN:)=0.
        call four2a(ct,2*NN,1,-1,1)
        df2=baud/(2*NN)
        ct=cshift(ct,NN)
        ppmax=0.
        dfpk=0.
        do i=0,2*NN-1
           f=(i-NN)*df2
           pp=real(ct(i))**2 + aimag(ct(i))**2
           if(pp.gt.ppmax) then
              ppmax=pp
              dfpk=f
           endif
        enddo

        zsum=0.
        do j=1,NN
           phi=(j-1)*twopi*dfpk*ts
           w=cmplx(cos(phi),sin(phi))
           cps=cs(i0:i0+6*nts:2*nts,j)*conjg(w)
           z=cps(id2(j))
           ct(j)=z
           zsum=zsum+z
           write(12,1042) j,id(j),id2(j),20*ps,atan2(aimag(z),real(z)),  &
                atan2(aimag(zsum),real(zsum)),zsum
1042       format(3i2,6f8.3,2f8.1)
        enddo

        phi0=atan2(aimag(zsum),real(zsum))
        zsum=0.
        do j=1,NN
           phi=(j-1)*twopi*dfpk*ts + phi0
           w=cmplx(cos(phi),sin(phi))
           nlo=0
           nhi=0
           cps=cs(i0:i0+6*nts:2*nts,j)*conjg(w)
           ps=real(cps)
           if(max(ps(1),ps(3)).ge.max(ps(0),ps(2))) nlo=1
           if(max(ps(2),ps(3)).ge.max(ps(0),ps(1))) nhi=1
           id2(j)=2*nhi+nlo
           z=cps(id2(j))
           ct(j)=z
           zsum=zsum+z
        enddo

        nh2=count(id.ne.id2)
        nb2=count(iand(id,1).ne.iand(id2,1)) + count(iand(id,2).ne.iand(id2,2))
        nhard1=nhard1+nh1
        nhard2=nhard2+nh2
        nbit1=nbit1+nb1
        nbit2=nbit2+nb2

        fdiff=1500.0+freq - f0
        write(13,1040)  snrdb,snr2500,f0,fdiff,width,nh1,nb1,nh2,nb2
1040    format(2f7.1,f9.2,f7.2,f6.1,2(i8,i6))
40      continue
     enddo

     ser1=float(nhard1)/(NN*iters)
     ser2=float(nhard2)/(NN*iters)
     ber1=float(nbit1)/(2*NN*iters)
     ber2=float(nbit2)/(2*NN*iters)
     write(*,1050)  snrdb,nhard1,nbit1,ser1,ber1,nhard2,nbit2,ser2,ber2
     write(14,1050) snrdb,nhard1,nbit1,ser1,ber1,nhard2,nbit2,ser2,ber2
1050 format(f6.1,2(2i5,2f8.4))
  enddo
  write(*,1060) NN*iters,2*NN*iters
1060 format(59('-')/'Max:  ',2i5)

999 end program fsk4sim


code:
subroutine genbpsk(id,f00,ndiff,nref,c)

  parameter (ND=121)                     !Data symbols: LDPC (120,60), r=1/2
  parameter (NN=ND)                      !Total symbols (121)
  parameter (NSPS=28800)                 !Samples per symbol at 12000 sps
  parameter (NZ=NSPS*NN)                 !Samples in waveform (3456000)

  complex c(0:NZ-1)                      !Complex waveform
  real*8 twopi,dt,fs,baud,f0,dphi,phi
  integer id(NN)                         !Encoded NRZ data (values +/-1)
  integer ie(NN)                         !Differentially encoded data

  f0=f00
  twopi=8.d0*atan(1.d0)
  fs=12000.d0
  dt=1.0/fs
  baud=1.d0/(NSPS*dt)

  if(ndiff.ne.0) then
     ie(1)=1                             !First bit is always 1
     do i=2,NN                           !Differentially encode
        ie(i)=id(i)*ie(i-1)
     enddo
  endif

! Generate the BPSK waveform
  phi=0.d0
  k=-1
  do j=1,NN
     dphi=twopi*f0*dt
     x=id(j)
     if(ndiff.ne.0) x=ie(j)                !Differential
     if(nref.ne.0) x=1.0                   !Generate reference carrier
     do i=1,NSPS
        k=k+1
        phi=phi+dphi
        if(phi.gt.twopi) phi=phi-twopi
        xphi=phi
        c(k)=x*cmplx(cos(xphi),sin(xphi))
     enddo
  enddo

  return
end subroutine genbpsk


code:
subroutine genfsk4(id,f00,nts,c)

  parameter (ND=60)                      !Data symbols: LDPC (120,60), r=1/2
  parameter (NN=ND)                      !Total symbols (60)
  parameter (NSPS=57600)                 !Samples per symbol at 12000 sps
  parameter (NZ=NSPS*NN)                 !Samples in waveform (3456000)
  parameter (NFFT=NZ)                    !Full length FFT

  complex c(0:NFFT-1)                    !Complex waveform
  real*8 twopi,dt,fs,baud,f0,dphi,phi
  integer id(NN)                         !Encoded 2-bit data (values 0-3)

  f0=f00
  twopi=8.d0*atan(1.d0)
  fs=12000.d0
  dt=1.0/fs
  baud=1.d0/(NSPS*dt)

! Generate the 4-FSK waveform
  x=0.
  c=0.
  phi=0.d0
  k=-1
  do j=1,NN
     dphi=twopi*(f0 + nts*id(j)*baud)*dt
     do i=1,NSPS
        k=k+1
        phi=phi+dphi
        if(phi.gt.twopi) phi=phi-twopi
        xphi=phi
        c(k)=cmplx(cos(xphi),sin(xphi))
     enddo
  enddo

  return
end subroutine genfsk4


code:
subroutine genfsk4hf(msgbits,f0,id,c)

  parameter (KK=84)                     !Information bits (72 + CRC12)
  parameter (ND=84)                     !Data symbols: LDPC (168,84), r=1/2
  parameter (NS=12)                     !Sync symbols (3 @ 4x4 Costas arrays)
  parameter (NR=2)                      !Ramp up/down
  parameter (NN=NR+NS+ND)               !Total symbols (98)
  parameter (NSPS=2688/84)              !Samples per symbol (32)
  parameter (NZ=NSPS*NN)                !Samples in baseband waveform (3136)

  complex c(0:NZ-1)                     !Complex waveform
  integer id0(NN)                       !2-bit data (values 0-3), all symbols
  integer id(ND)                        !2-bit data (values 0-3), data only
  integer*1 msgbits(KK),codeword(2*ND)
  integer icos4(4)                      !4x4 Costas array
  data icos4/0,1,3,2/

  twopi=8.0*atan(1.0)
  fs=12000.0/84.0
  dt=1.0/fs
  baud=1.0/(NSPS*dt)
  call encode168(msgbits,codeword)      !Encode the test message
  id0(1)=0                               !Ramp-up
  id0(2:5)=icos4                         !First Costas array
  id0(48:51)=icos4                       !Second
  id0(94:97)=icos4                       !Third
  id0(98)=0                              !Ramp down
  j=5
  do i=1,84                             !Data symbols
     id(i)=2*codeword(2*i-1) + codeword(2*i)
     j=j+1
     if(i.eq.43) j=j+4
     id0(j)=id(i)
  enddo

! Generate the 4-FSK waveform, low tone at f=0
  c=0.
  phi=0.d0
  k=-1
  do j=1,NN
     dphi=twopi*(f0+id0(j)*baud)*dt
     do i=1,NSPS
        k=k+1
        phi=phi+dphi
        if(phi.gt.twopi) phi=phi-twopi
        c(k)=cmplx(cos(phi),sin(phi))
     enddo
  enddo

  return
end subroutine genfsk4hf


code:
subroutine genmskhf(msgbits,id,icw,cbb,csync)

!Encode an MSK-HF message, produce baseband waveform and sync vector.

  parameter (KK=84)                     !Information bits (72 + CRC12)
  parameter (ND=168)                    !Data symbols: LDPC (168,84), r=1/2
  parameter (NS=65)                     !Sync symbols (2 x 26 + Barker 13)
  parameter (NR=3)                      !Ramp up/down
  parameter (NN=NR+NS+ND)               !Total symbols (236)
  parameter (NSPS=16)                   !Samples per MSK symbol (16)
  parameter (N2=2*NSPS)                 !Samples per OQPSK symbol (32)
  parameter (NZ=NSPS*NN)                !Samples in baseband waveform (3760)

  complex cbb(0:NZ-1)
  complex csync(0:NZ-1)
  real x(0:NZ-1)
  real y(0:NZ-1)
  real pp(N2)
  logical first
  integer*1 msgbits(KK),codeword(ND)
  integer icw(ND)
  integer id(NS+ND)
  integer isync(26)                          !Long sync vector
  integer ib13(13)                           !Barker 13 code
  data ib13/1,1,1,1,1,-1,-1,1,1,-1,1,-1,1/
  data first/.true./
  save first,isync,twopi,pp

  if(first) then
     n=z'2c1aeb1'
     do i=1,26
        isync(i)=-1
        if(iand(n,1).eq.1) isync(i)=1
        n=n/2
     enddo

     twopi=8.0*atan(1.0)
     do i=1,N2                             !Half-sine shaped pulse
        pp(i)=sin(0.5*(i-1)*twopi/N2)
     enddo
     first=.false.
  endif

  call encode168(msgbits,codeword)      !Encode the test message
  icw=2*codeword - 1

! Message structure: R1 26*(S1+D1) S13 26*(D1+S1) R1
! Generate QPSK without any offset; then shift the y array to get OQPSK.

! Do the I channel first: results in array x
  n=0
  k=0
  ia=0
  ib=NSPS-1
  x(ia:ib)=0.                           !Ramp up (half-symbol; shape TBD)
  do j=1,26                             !Insert group of 26*(S1+D1)
     ia=ib+1
     ib=ia+N2-1
     n=n+1
     id(n)=2*isync(j)
     x(ia:ib)=isync(j)*pp               !Insert Sync bit
     ia=ib+1
     ib=ia+N2-1
     k=k+1
     n=n+1
     id(n)=icw(k)
     x(ia:ib)=id(n)*pp                  !Insert data bit
  enddo

  do j=1,13                             !Insert Barker 13 code
     ia=ib+1
     ib=ia+N2-1
     n=n+1
     id(n)=2*ib13(j)
     x(ia:ib)=ib13(j)*pp
  enddo

  do j=1,26                             !Insert group of 26*(S1+D1)
     ia=ib+1
     ib=ia+N2-1
     k=k+1
     n=n+1
     id(n)=icw(k)
     x(ia:ib)=id(n)*pp                  !Insert data bit
     ia=ib+1
     ib=ia+N2-1
     n=n+1
     id(n)=2*isync(j)
     x(ia:ib)=isync(j)*pp               !Insert Sync bit
  enddo
  ia=ib+1
  ib=ia+NSPS-1
  x(ia:ib)=0.                           !Ramp down (half-symbol; shape TBD)

! Now do the Q channel: results in array y
  ia=0
  ib=NSPS-1
  y(ia:ib)=0.                           !Ramp up  (half-symbol; shape TBD)
  do j=1,116
     ia=ib+1
     ib=ia+N2-1
     k=k+1
     n=n+1
     id(n)=icw(k)
     y(ia:ib)=id(n)*pp
  enddo
  ia=ib+1
  ib=ia+NSPS-1
  y(ia:ib)=0.                          !Ramp down (half-symbol; shape TBD)
  y=cshift(y,-NSPS)                    !Shift Q array to get OQPSK
  cbb=cmplx(x,y)                       !Complex baseband waveform

  ib=NSPS-1
  ib2=NSPS-1+64*N2
  do j=1,26                            !Zero all data symbols in x
     ia=ib+1+N2
     ib=ia+N2-1
     x(ia:ib)=0.
     ia2=ib2+1+N2
     ib2=ia2+N2-1
     x(ia2:ib2)=0.
  enddo
  csync=x

  return
end subroutine genmskhf


code:
subroutine genwspr5(msg,msgsent,itone)

! Encode a WSPR-LF message, producing array itone().

  use crc
  include 'wsprlf_params.f90'

  character*22 msg,msgsent
  character*60 cbits
  integer*1,target :: idat(9)
  integer*1 msgbits(KK),codeword(ND)
  logical first
  integer icw(ND)
  integer id(NS+ND)
  integer jd(NS+ND)
  integer isync(48)                          !Long sync vector
  integer ib13(13)                           !Barker 13 code
  integer itone(NN)
  integer*8 n8
  data ib13/1,1,1,1,1,-1,-1,1,1,-1,1,-1,1/
  data first/.true./
  save first,isync

  if(first) then
     n8=z'cbf089223a51'
     do i=1,48
        isync(i)=-1
        if(iand(n8,1).eq.1) isync(i)=1
        n8=n8/2
     enddo
     first=.false.
  endif

  idat=0
  call wqencode(msg,ntype0,idat)             !Source encoding
  id7=idat(7)
  if(id7.lt.0) id7=id7+256
  id7=id7/64
  icrc=crc10(c_loc(idat),9)                  !Compute the 10-bit CRC
  idat(8)=icrc/256                           !Insert CRC into idat(8:9)
  idat(9)=iand(icrc,255)
  call wqdecode(idat,msgsent,itype)

  write(cbits,1004) idat(1:6),id7,icrc
1004 format(6b8.8,b2.2,b10.10)
  read(cbits,1006) msgbits
1006 format(60i1)

!  call chkcrc10(msgbits,nbadcrc)
!  print*,msgsent,itype,crc10_check(c_loc(idat),9),nbadcrc

  call encode300(msgbits,codeword)      !Encode the test message
  icw=2*codeword - 1                    !NRZ codeword

! Message structure:
! I channel:  R1 48*(S1+D1) S13 48*(D1+S1) R1
! Q channel:  R1 D204 R1
! Generate QPSK with no offset, then shift the y array to get OQPSK.

! I channel:
  n=0
  k=0
  do j=1,48                             !Insert group of 48*(S1+D1)
     n=n+1
     id(n)=2*isync(j)
     k=k+1
     n=n+1
     id(n)=icw(k)
  enddo

  do j=1,13                             !Insert Barker 13 code
     n=n+1
     id(n)=2*ib13(j)
  enddo

  do j=1,48                             !Insert group of 48*(S1+D1)
     k=k+1
     n=n+1
     id(n)=icw(k)
     n=n+1
     id(n)=2*isync(j)
  enddo

! Q channel
  do j=1,204
     k=k+1
     n=n+1
     id(n)=icw(k)
  enddo

! Map I and Q to tones.
  n=0
  jz=(NS+ND+1)/2
  do j=1,jz-1
     jd(2*j-1)=id(j)/abs(id(j))
     jd(2*j)=id(j+jz)/abs(id(j+jz))
  enddo
  jd(NS+ND)=id(jz)/abs(id(jz))
  itone=0
  do j=1,jz-1
     itone(2*j+1)=(jd(2*j)*jd(2*j-1)+1)/2;
     itone(2*j+2)=-(jd(2*j)*jd(2*j+1)-1)/2;
  enddo
  itone(NS+ND+2)=jd(NS+ND)                       !### Is this correct ??? ###

  return
end subroutine genwspr5


code:
subroutine genwspr_fsk8(msg,msgsent,itone)

! Encode a WSPR-LF 8-FSK message, producing array itone().

  use crc
  include 'wspr_fsk8_params.f90'

  character*22 msg,msgsent
  character*60 cbits
  integer*1,target :: idat(9)
  integer*1 msgbits(KK),codeword(3*ND)
  integer itone(NN)
  integer icos7(0:6)
  data icos7/2,5,6,0,4,1,3/                  !Costas 7x7 tone pattern

  idat=0
  call wqencode(msg,ntype0,idat)             !Source encoding
  id7=idat(7)
  if(id7.lt.0) id7=id7+256
  id7=id7/64
  icrc=crc10(c_loc(idat),9)                  !Compute the 10-bit CRC
  idat(8)=icrc/256                           !Insert CRC into idat(8:9)
  idat(9)=iand(icrc,255)
  call wqdecode(idat,msgsent,itype)

  write(cbits,1004) idat(1:6),id7,icrc
1004 format(6b8.8,b2.2,b10.10)
  read(cbits,1006) msgbits
1006 format(60i1)

!  call chkcrc10(msgbits,nbadcrc)
!  print*,msgsent,itype,crc10_check(c_loc(idat),9),nbadcrc

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

! Message structure: S7 D100 S7
  itone(1:7)=icos7
  itone(NN-6:NN)=icos7
  do j=1,ND
     i=3*j -2
     itone(j+7)=codeword(i)*4 + codeword(i+1)*2 + codeword(i+2)
  enddo

  return
end subroutine genwspr_fsk8


code:
subroutine genwsprcpm(msg,msgsent,itone)

! Encode a WSPRCPM message, producing array itone().
!
  use crc
  include 'wsprcpm_params.f90'

  character*22 msg,msgsent
  character*64 cbits
  character*32 sbits
  integer      iuniqueword0
  integer*1,target :: idat(9)
  integer*1 msgbits(68),codeword(ND)
  logical first
  integer icw(ND)
  integer id(NS+ND)
  integer jd(NS+ND)
  integer ipreamble(16)                      !Freq estimation preamble
  integer isync(200)                          !Long sync vector
  integer itone(NN)
  data ipreamble/1,1,1,1,0,0,0,0,0,0,0,0,1,1,1,1/
  data first/.true./
  data iuniqueword0/z'30C9E8AD'/
  save first,isync,ipreamble

  if(first) then
    write(sbits,'(b32.32)') iuniqueword0
    read(sbits,'(32i1)') isync(1:32)
    read(sbits,'(32i1)') isync(33:64)
    read(sbits,'(32i1)') isync(65:96)
    read(sbits,'(32i1)') isync(97:128)
    read(sbits,'(32i1)') isync(129:160)
    read(sbits,'(32i1)') isync(161:192)
    read(sbits,'(8i1)') isync(193:200)
    first=.false.
  endif

  idat=0
  call wqencode(msg,ntype0,idat)             !Source encoding
  id7=idat(7)
  if(id7.lt.0) id7=id7+256
  id7=id7/64
write(*,*) 'idat ',idat
  icrc=crc14(c_loc(idat),9)
write(*,*) 'icrc: ',icrc
write(*,'(a6,b16.16)') 'icrc: ',icrc
  call wqdecode(idat,msgsent,itype)
  print*,msgsent,itype
  write(cbits,1004) idat(1:6),id7,iand(icrc,z'3FFF')
1004 format(6b8.8,b2.2,b14.14)
  msgbits=0
  read(cbits,1006) msgbits(1:64)
1006 format(64i1)

write(*,'(50i1,1x,14i1,1x,4i1)') msgbits

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

! Message structure:
! d100 p16 d100
  itone(1:100)=isync(1:100)+2*codeword(1:100)
  itone(101:116)=ipreamble+1
  itone(117:216)=isync(101:200)+2*codeword(101:200)
  itone=2*itone-3
  do i=1,216
    write(*,*) i,itone(i)
  enddo

  return
end subroutine genwsprcpm


code:
subroutine genwsprdpsk(msg,msgsent,imsgde)

! Encode a WSPRDPSK message, producing array txwave().
!
  use crc
  include 'wsprdpsk_params.f90'

  character*22 msg,msgsent
  character*64 cbits
  character*32 sbits
  integer      iuniqueword0
  integer*1,target :: idat(9)
  integer*1 msgbits(68),codeword(ND)
  logical first
  integer ipreamble(16)                      !Freq estimation preamble
  integer isync(32)                          !Long sync vector
  integer imsg(NN),imsgde(NN)
  data ipreamble/1,1,1,1,0,0,0,0,0,0,0,0,1,1,1,1/
  data first/.true./
  data iuniqueword0/z'30C9E8AD'/
  save first,isync,ipreamble

  if(first) then
    write(sbits,'(b32.32)') iuniqueword0
    read(sbits,'(32i1)') isync(1:32)
    first=.false.
  endif

  idat=0
  call wqencode(msg,ntype0,idat)             !Source encoding
  id7=idat(7)
  if(id7.lt.0) id7=id7+256
  id7=id7/64
write(*,*) 'idat ',idat
  icrc=crc14(c_loc(idat),9)
write(*,*) 'icrc: ',icrc
write(*,'(a6,b16.16)') 'icrc: ',icrc
  call wqdecode(idat,msgsent,itype)
  print*,msgsent,itype
  write(cbits,1004) idat(1:6),id7,iand(icrc,z'3FFF')
1004 format(6b8.8,b2.2,b14.14)
  msgbits=0
  read(cbits,1006) msgbits(1:64)
1006 format(64i1)

write(*,'(50i1,1x,14i1,1x,4i1)') msgbits

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

  imsg(1)=1                         !reference bit
  imsg(2:101)=codeword(1:100)
  imsg(102:132)=isync(1:31)         !only use 31 of the sync bits
  imsg(133:232)=codeword(101:200)
write(*,'(232i1)') imsg(1:232)
  imsgde(1)=1
  do i=2,232
    imsgde(i)=mod(imsgde(i-1)+imsg(i),2)
  enddo

write(*,*) '-------------'
write(*,'(232i1)') imsgde(1:232)
  return
end subroutine genwsprdpsk


code:
subroutine genwsprlf(msgbits,id,icw,cbb,csync,itone)

!Encode a WSPR-LF message, produce baseband waveform and sync vector.

  include 'wsprlf_params.f90'

  complex cbb(0:NZ-1)
  complex csync(0:NZ-1)
  real x(0:NZ-1)
  real y(0:NZ-1)
  real pp(N2)
  logical first
  integer*1 msgbits(KK),codeword(ND)
  integer icw(ND)
  integer id(NS+ND)
  integer jd(NS+ND)
  integer isync(48)                          !Long sync vector
  integer ib13(13)                           !Barker 13 code
  integer itone(NN)
  integer*8 n8
  data ib13/1,1,1,1,1,-1,-1,1,1,-1,1,-1,1/
  data first/.true./
  save first,isync,twopi,pp

  if(first) then
     n8=z'cbf089223a51'
     do i=1,48
        isync(i)=-1
        if(iand(n8,1).eq.1) isync(i)=1
        n8=n8/2
     enddo

     twopi=8.0*atan(1.0)
     do i=1,N2                             !Half-sine shaped pulse
        pp(i)=sin(0.5*(i-1)*twopi/N2)
     enddo
     first=.false.
  endif

  call encode300(msgbits,codeword)      !Encode the test message
  icw=2*codeword - 1

! Message structure: R1 48*(S1+D1) S13 48*(D1+S1) R1
! Generate QPSK without any offset; then shift the y array to get OQPSK.

! Do the I channel first: results in array x
  n=0
  k=0
  ia=0
  ib=NSPS-1
  x(ia:ib)=0.                           !Ramp up (half-symbol; shape TBD)
  do j=1,48                             !Insert group of 48*(S1+D1)
     ia=ib+1
     ib=ia+N2-1
     n=n+1
     id(n)=2*isync(j)
     x(ia:ib)=isync(j)*pp               !Insert Sync bit
     ia=ib+1
     ib=ia+N2-1
     k=k+1
     n=n+1
     id(n)=icw(k)
     x(ia:ib)=id(n)*pp                  !Insert data bit
  enddo

  do j=1,13                             !Insert Barker 13 code
     ia=ib+1
     ib=ia+N2-1
     n=n+1
     id(n)=2*ib13(j)
     x(ia:ib)=ib13(j)*pp
  enddo

  do j=1,48                             !Insert group of 48*(S1+D1)
     ia=ib+1
     ib=ia+N2-1
     k=k+1
     n=n+1
     id(n)=icw(k)
     x(ia:ib)=id(n)*pp                  !Insert data bit
     ia=ib+1
     ib=ia+N2-1
     n=n+1
     id(n)=2*isync(j)
     x(ia:ib)=isync(j)*pp               !Insert Sync bit
  enddo
  ia=ib+1
  ib=ia+NSPS-1
  x(ia:ib)=0.                           !Ramp down (half-symbol; shape TBD)

! Now do the Q channel: results in array y
  ia=0
  ib=NSPS-1
  y(ia:ib)=0.                           !Ramp up  (half-symbol; shape TBD)
  do j=1,204
     ia=ib+1
     ib=ia+N2-1
     k=k+1
     n=n+1
     id(n)=icw(k)
     y(ia:ib)=id(n)*pp
  enddo
  ia=ib+1
  ib=ia+NSPS-1
  y(ia:ib)=0.                          !Ramp down (half-symbol; shape TBD)
  y=cshift(y,-NSPS)                    !Shift Q array to get OQPSK
  cbb=cmplx(x,y)                       !Complex baseband waveform

  ib=NSPS-1
  ib2=NSPS-1+64*N2
  do j=1,48                            !Zero all data symbols in x
     ia=ib+1+N2
     ib=ia+N2-1
     x(ia:ib)=0.
     ia2=ib2+1+N2
     ib2=ia2+N2-1
     x(ia2:ib2)=0.
  enddo
  csync=x

! Map I and Q to tones.
  n=0
  jz=(NS+ND+1)/2
  do j=1,jz-1
     jd(2*j-1)=id(j)/abs(id(j))
     jd(2*j)=id(j+jz)/abs(id(j+jz))
  enddo
  jd(NS+ND)=id(jz)/abs(id(jz))
  itone=0
  do j=1,jz-1
     itone(2*j-1)=(jd(2*j)*jd(2*j-1)+1)/2;
     itone(2*j)=-(jd(2*j)*jd(2*j+1)-1)/2;
  enddo
  itone(NS+ND)=jd(NS+ND)                       !### Is this correct ??? ###

  return
end subroutine genwsprlf


code:
subroutine getfc1(c,fc1)

  parameter (KK=84)                     !Information bits (72 + CRC12)
  parameter (ND=168)                    !Data symbols: LDPC (168,84), r=1/2
  parameter (NS=65)                     !Sync symbols (2 x 26 + Barker 13)
  parameter (NR=3)                      !Ramp up/down
  parameter (NN=NR+NS+ND)               !Total symbols (236)
  parameter (NSPS=16)                   !Samples per MSK symbol (16)
  parameter (N2=2*NSPS)                 !Samples per OQPSK symbol (32)
  parameter (N13=13*N2)                 !Samples in central sync vector (416)
  parameter (NZ=NSPS*NN)                !Samples in baseband waveform (3760)
  parameter (NFFT1=4*NSPS,NH1=NFFT1/2)

  complex c(0:NZ-1)                     !Complex waveform
  complex c2(0:NFFT1-1)                 !Short spectra
  real s(-NH1+1:NH1)                    !Coarse spectrum

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

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

  return
end subroutine getfc1


code:
subroutine getfc1w(c,fs,fa,fb,fc1,xsnr)

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

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

  return
end subroutine getfc1w


code:
subroutine getfc2(c,csync,fc1,fc2,fc3)

  parameter (KK=84)                     !Information bits (72 + CRC12)
  parameter (ND=168)                    !Data symbols: LDPC (168,84), r=1/2
  parameter (NS=65)                     !Sync symbols (2 x 26 + Barker 13)
  parameter (NR=3)                      !Ramp up/down
  parameter (NN=NR+NS+ND)               !Total symbols (236)
  parameter (NSPS=16)                   !Samples per MSK symbol (16)
  parameter (N2=2*NSPS)                 !Samples per OQPSK symbol (32)
  parameter (N13=13*N2)                 !Samples in central sync vector (416)
  parameter (NZ=NSPS*NN)                !Samples in baseband waveform (3760)
  parameter (NFFT1=4*NSPS,NH1=NFFT1/2)

  complex c(0:NZ-1)                     !Complex waveform
  complex cs(0:NZ-1)                    !For computing spectrum
  complex csync(0:NZ-1)                 !Sync symbols only, from cbb
  real a(5)

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

! Filter, square, then FFT to get refined carrier frequency fc2.
  call four2a(cs,NZ,1,-1,1)          !To freq domain
  ia=nint(0.75*baud/df)
  cs(ia:NZ-1-ia)=0.                  !Save only freqs around fc1
  call four2a(cs,NZ,1,1,1)           !Back to time domain
  cs=cs/NZ
  cs=cs*cs                           !Square the data
  call four2a(cs,NZ,1,-1,1)          !Compute squared spectrum

! Find two peaks separated by baud
  pmax=0.
  fc2=0.
  ic=nint(baud/df)
  ja=nint(0.5*baud/df)
  do j=-ja,ja
     f2=j*df
     ia=nint((f2-0.5*baud)/df)
     if(ia.lt.0) ia=ia+NZ
     ib=nint((f2+0.5*baud)/df)
     p=real(cs(ia))**2 + aimag(cs(ia))**2 +                        &
          real(cs(ib))**2 + aimag(cs(ib))**2
     if(p.gt.pmax) then
        pmax=p
        fc2=0.5*f2
     endif
!           write(52,1200) f2,p,db(p)
!1200       format(f10.3,2f15.3)
  enddo

  a(1)=-fc1
  a(2:5)=0.
  call twkfreq1(c,NZ,fs,a,cs)         !Mix down by fc1
  cs=cs*conjg(csync)
  call four2a(cs,NZ,1,-1,1)          !To freq domain
  pmax=0.
  do i=0,NZ-1
     f=i*df
     if(i.gt.NZ/2) f=(i-NZ)*df
     p=real(cs(i))**2 + aimag(cs(i))**2
!     write(51,3001) f,p,db(p)
!3001 format(f10.3,e12.3,f10.3)
     if(p.gt.pmax) then
        pmax=p
        fc3=f
     endif
  enddo

  return
end subroutine getfc2


code:
subroutine getfc2w(c,csync,npeaks,fs,fc1,fpks)

  include 'wsprlf_params.f90'

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

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

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

  ia=nint(0.75*baud/df)
  cs(ia:NZ-1-ia)=0.                  !Save only freqs around fc1
  call four2a(cs,NZ,1,1,1)           !Back to time domain
  cs=cs/NZ
  cs=cs*cs                           !Square the data
  call four2a(cs,NZ,1,-1,1)          !Compute squared spectrum

! Find two peaks separated by baud
  pmax=0.
  fc2=0.
  ja=nint(0.3*baud/df)
  k=1
  do j=-ja,ja
     f2=j*df
     ia=nint((f2-0.5*baud)/df)
     if(ia.lt.0) ia=ia+NZ
     ib=nint((f2+0.5*baud)/df)
     p=real(cs(ia))**2 + aimag(cs(ia))**2 +                        &
          real(cs(ib))**2 + aimag(cs(ib))**2
     if(p.gt.pmax) then
        pmax=p
        fc2=0.5*f2
     endif
     freqs(k)=0.5*f2
     sp2(k)=p
     k=k+1
!           write(52,1200) f2,p,db(p)
!1200       format(f10.3,2f15.3)
  enddo

  do i=1,npeaks
    pkloc=maxloc(sp2)
    ipk=pkloc(1)
    fpks(i)=freqs(ipk)
    ipk0=max(1,ipk-1)
    ipk1=min(413,ipk+1)
!    ipk0=ipk
!    ipk1=ipk
    sp2(ipk0:ipk1)=0.0
!write(*,*) i,fpks(i),fc2
  enddo

  a(1)=-fc1
  a(2:5)=0.
  call twkfreq1(c,NZ,fs,a,cs)         !Mix down by fc1
  cs=cs*conjg(csync)
  call four2a(cs,NZ,1,-1,1)          !To freq domain
  pmax=0.
  do i=0,NZ-1
     f=i*df
     if(i.gt.NZ/2) f=(i-NZ)*df
     p=real(cs(i))**2 + aimag(cs(i))**2
!     write(51,3001) f,p,db(p)
!3001 format(f10.3,e12.3,f10.3)
     if(p.gt.pmax) then
        pmax=p
        fc3=f
     endif
  enddo

  return
end subroutine getfc2w


code:
#include <stdlib.h>
#include <math.h>

/* Generate gaussian random float with mean=0 and std_dev=1 */
float gran_()
{
  float fac,rsq,v1,v2;
  static float gset;
  static int iset;

  if(iset){
    /* Already got one */
    iset = 0;
    return gset;
  }
  /* Generate two evenly distributed numbers between -1 and +1
   * that are inside the unit circle
   */
  do {
    v1 = 2.0 * (float)rand() / RAND_MAX - 1;
    v2 = 2.0 * (float)rand() / RAND_MAX - 1;
    rsq = v1*v1 + v2*v2;
  } while(rsq >= 1.0 || rsq == 0.0);
  fac = sqrt(-2.0*log(rsq)/rsq);
  gset = v1*fac;
  iset++;
  return v2*fac;
}


code:
program ldpcsim120
! End to end test of the (120,60)/crc10 encoder and decoder.
use crc
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(9)
integer*1, target:: i1Dec8BitBytes(9)
integer*1 msgbits(60)
integer*1 apmask(120)
integer*1 cw(120)
integer*2 checksum
integer colorder(120)
integer nerrtot(120),nerrdec(120),nmpcbad(60)
logical checksumok,fsk,bpsk
real*8, allocatable ::  rxdata(:)
real, allocatable :: llr(:)
real dllr(120),llrd(120)

data colorder/   &
  0,1,2,21,3,4,5,6,7,8,20,10,9,11,12,23,13,28,14,31, &
  15,16,22,26,17,30,18,29,25,32,41,34,19,33,27,36,38,43,42,24, &
  37,39,45,40,35,44,47,46,50,51,53,48,52,56,54,57,55,49,58,61, &
  60,59,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79, &
  80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99, &
  100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119/

do i=1,NRECENT
  recent_calls(i)='            '
enddo
nerrtot=0
nerrdec=0
nmpcbad=0  ! Used to collect the number of errors in the message+crc part of the codeword

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

fsk=.false.
bpsk=.true.

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

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

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

! The message should be packed into the first 7 bytes
  i1Msg8BitBytes(1:6)=85
  i1Msg8BitBytes(7)=64
! The CRC will be put into the last 2 bytes
  i1Msg8BitBytes(8:9)=0
  checksum = crc10 (c_loc (i1Msg8BitBytes), 9)
! For reference, the next 3 lines show how to check the CRC
  i1Msg8BitBytes(8)=checksum/256
  i1Msg8BitBytes(9)=iand (checksum,255)
  checksumok = crc10_check(c_loc (i1Msg8BitBytes), 9)
  if( checksumok ) write(*,*) 'Good checksum'
write(*,*) i1Msg8BitBytes(1:9)

  mbit=0
  do i=1, 7
    i1=i1Msg8BitBytes(i)
    do ibit=1,8
      mbit=mbit+1
      msgbits(mbit)=iand(1,ishft(i1,ibit-8))
    enddo
  enddo
  i1=i1Msg8BitBytes(8) ! First 2 bits of crc10 are LSB of this byte
  do ibit=1,2
    msgbits(50+ibit)=iand(1,ishft(i1,ibit-2))
  enddo
  i1=i1Msg8BitBytes(9) ! Now shift in last 8 bits of the CRC
  do ibit=1,8
    msgbits(52+ibit)=iand(1,ishft(i1,ibit-8))
  enddo

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

  call encode120(msgbits,codeword)
  call init_random_seed()
  call sgran()

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

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

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

    llr=2.0*rxdata/(ss*ss)
    apmask=0

! max_iterations is max number of belief propagation iterations
    call bpdecode120(llr, apmask, max_iterations, decoded, niterations, cw)
    n2err=0
    do i=1,N
      if( cw(i)*(2*codeword(i)-1.0) .lt. 0 ) n2err=n2err+1
    enddo
!write(*,*) nerr,niterations,n2err
    damp=0.75
    ndither=0
    if( niterations .lt. 0 ) then
      do i=1, ndither
        do in=1,N
          dllr(in)=damp*gran()
        enddo
        llrd=llr+dllr
        call bpdecode120(llrd, apmask, max_iterations, decoded, niterations, cw)
        if( niterations .ge. 0 ) exit
      enddo
    endif

! If the decoder finds a valid codeword, niterations will be .ge. 0.
    if( niterations .ge. 0 ) then
! Check the CRC
      do ibyte=1,6
        itmp=0
        do ibit=1,8
          itmp=ishft(itmp,1)+iand(1,decoded((ibyte-1)*8+ibit))
        enddo
        i1Dec8BitBytes(ibyte)=itmp
      enddo
      i1Dec8BitBytes(7)=decoded(49)*128+decoded(50)*64
! Need to pack the received crc into bytes 8 and 9 for crc10_check
      i1Dec8BitBytes(8)=decoded(51)*2+decoded(52)
      i1Dec8BitBytes(9)=decoded(53)*128+decoded(54)*64+decoded(55)*32+decoded(56)*16
      i1Dec8BitBytes(9)=i1Dec8BitBytes(9)+decoded(57)*8+decoded(58)*4+decoded(59)*2+decoded(60)*1
      ncrcflag=0
      if( crc10_check( c_loc( i1Dec8BitBytes ), 9 ) ) ncrcflag=1

      if( ncrcflag .ne. 1 ) then
        nbadcrc=nbadcrc+1
      endif
      nueflag=0

      nerrmpc=0
      do i=1,K   ! find number of errors in message+crc part of codeword
        if( msgbits(i) .ne. decoded(i) ) then
          nueflag=1
          nerrmpc=nerrmpc+1
        endif
      enddo
      nmpcbad(nerrmpc)=nmpcbad(nerrmpc)+1  ! This histogram should inform our selection of CRC poly
      if( ncrcflag .eq. 1 .and. nueflag .eq. 0 ) then
        ngood=ngood+1
        nerrdec(nerr)=nerrdec(nerr)+1
      else if( ncrcflag .eq. 1 .and. nueflag .eq. 1 ) then
        nue=nue+1;
      endif
    endif
  enddo
  snr2500=db+10*log10(0.4166/2500.0)
  pberr=real(nberr)/(real(ntrials*N))
  write(*,"(f4.1,4x,f5.1,1x,i8,1x,i8,1x,i8,8x,f5.2,8x,e10.3)") db,snr2500,ngood,nue,nbadcrc,ss,pberr

enddo

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

end program ldpcsim120


code:
program ldpcsim168
! End to end test of the (168,84)/crc12 encoder and decoder.
use crc
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(11)
integer*1 msgbits(84)
integer*1 apmask(168), cw(168)
integer*2 checksum
integer*4 i4Msg6BitWords(13)
integer colorder(168)
integer nerrtot(168),nerrdec(168),nmpcbad(84)
logical checksumok,fsk,bpsk
real*8, allocatable ::  rxdata(:)
real, allocatable :: llr(:)

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

do i=1,NRECENT
  recent_calls(i)='            '
enddo
nerrtot=0
nerrdec=0
nmpcbad=0  ! Used to collect the number of errors in the message+crc part of the codeword

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

fsk=.false.
bpsk=.true.

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

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

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

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

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

  i1Msg8BitBytes(10:11)=0
  checksum = crc12 (c_loc (i1Msg8BitBytes), 11)
! For reference, the next 3 lines show how to check the CRC
  i1Msg8BitBytes(10)=checksum/256
  i1Msg8BitBytes(11)=iand (checksum,255)
  checksumok = crc12_check(c_loc (i1Msg8BitBytes), 11)
  if( checksumok ) write(*,*) 'Good checksum'

  mbit=0
  do i=1, 9
    i1=i1Msg8BitBytes(i)
    do ibit=1,8
      mbit=mbit+1
      msgbits(mbit)=iand(1,ishft(i1,ibit-8))
    enddo
  enddo
  i1=i1Msg8BitBytes(10) ! First 4 bits of crc12 are LSB of this byte
  do ibit=1,4
    msgbits(72+ibit)=iand(1,ishft(i1,ibit-4))
  enddo
  i1=i1Msg8BitBytes(11) ! Now shift in last 8 bits of the CRC
  do ibit=1,8
    msgbits(76+ibit)=iand(1,ishft(i1,ibit-8))
  enddo

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

  call encode168(msgbits,codeword)
  call init_random_seed()
  call sgran()

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

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

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

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

! max_iterations is max number of belief propagation iterations
    call bpdecode168(llr, apmask, max_iterations, decoded, niterations)
!    if( niterations .eq. -1 ) then
!      norder=3
!      call osd168(llr, norder, decoded, niterations, cw)
!    endif
! If the decoder finds a valid codeword, niterations will be .ge. 0.
    if( niterations .ge. 0 ) then
      call extractmessage168(decoded,msgreceived,ncrcflag,recent_calls,nrecent)
      if( ncrcflag .ne. 1 ) then
        nbadcrc=nbadcrc+1
      endif

      nueflag=0
      nerrmpc=0
      do i=1,K   ! find number of errors in message+crc part of codeword
        if( msgbits(i) .ne. decoded(i) ) then
          nueflag=1
          nerrmpc=nerrmpc+1
        endif
      enddo

write(37,*) niterations, ncrcflag, nueflag
      nmpcbad(nerrmpc)=nmpcbad(nerrmpc)+1
      if( ncrcflag .eq. 1 ) then
        if( nueflag .eq. 0 ) then
          ngood=ngood+1
          nerrdec(nerr)=nerrdec(nerr)+1
        else if( nueflag .eq. 1 ) then
          nue=nue+1;
        endif
      endif
    endif
  enddo
  snr2500=db+10*log10(10.417/2500.0)
  pberr=real(nberr)/(real(ntrials*N))
  write(*,"(f4.1,4x,f5.1,1x,i8,1x,i8,1x,i8,8x,f5.2,8x,e10.3)") db,snr2500,ngood,nue,nbadcrc,ss,pberr

enddo

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

end program ldpcsim168


code:
program msksim

! Simulate characteristics of a potential "MSK10" mode using LDPC (168,84)
! code, OQPDK modulation, and 30 s T/R sequences.

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

  parameter (KK=84)                     !Information bits (72 + CRC12)
  parameter (ND=168)                    !Data symbols: LDPC (168,84), r=1/2
  parameter (NS=65)                     !Sync symbols (2 x 26 + Barker 13)
  parameter (NR=3)                      !Ramp up/down
  parameter (NN=NR+NS+ND)               !Total symbols (236)
  parameter (NSPS=1152/72)              !Samples per MSK symbol (16)
  parameter (N2=2*NSPS)                 !Samples per OQPSK symbol (32)
  parameter (N13=13*N2)                 !Samples in central sync vector (416)
  parameter (NZ=NSPS*NN)                !Samples in baseband waveform (3776)
  parameter (NFFT1=4*NSPS,NH1=NFFT1/2)

  character*8 arg
  complex cbb(0:NZ-1)                   !Complex baseband waveform
  complex csync(0:NZ-1)                 !Sync symbols only, from cbb
  complex cb13(0:N13-1)                 !Barker 13 waveform
  complex c(0:NZ-1)                     !Complex waveform
  complex c0(0:NZ-1)                    !Complex waveform
  complex zz(NS+ND)                     !Complex symbol values (intermediate)
  complex z
  real xnoise(0:NZ-1)                   !Generated random noise
  real ynoise(0:NZ-1)                   !Generated random noise
  real rxdata(ND),llr(ND)               !Soft symbols
  real pp(2*NSPS)                       !Shaped pulse for OQPSK
  real a(5)                             !For twkfreq1
  real aa(20),bb(20)                    !Fitted polyco's
  integer id(NS+ND)                     !NRZ values (+/-1) for Sync and Data
  integer ierror(NS+ND)
  integer icw(NN)
  integer*1 msgbits(KK),decoded(KK),apmask(ND),cw(ND)
!  integer*1 codeword(ND)
  data msgbits/0,0,1,0,0,1,1,1,1,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,1,1,0,0,0,1, &
       1,1,1,0,1,1,1,1,1,1,1,0,0,1,0,0,1,1,0,1,0,1,1,1,0,1,1,0,1,1,         &
       1,1,0,1,0,1,1,0,0,0,0,0,1,0,0,0,0,0,1,0,1,0,1,0/

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

  twopi=8.0*atan(1.0)
  fs=12000.0/72.0                        !Sample rate = 166.6666667 Hz
  dt=1.0/fs                              !Sample interval (s)
  tt=NSPS*dt                             !Duration of "itone" symbols (s)
  ts=2*NSPS*dt                           !Duration of OQPSK symbols (s)
  baud=1.0/tt                            !Keying rate for "itone" symbols (baud)
  txt=NZ*dt                              !Transmission length (s)
  bandwidth_ratio=2500.0/(fs/2.0)
  write(*,1000) f0,delay,fspread,maxn,iters,baud,3*baud,txt
1000 format('f0:',f5.1,'  Delay:',f4.1,'  fSpread:',f5.2,'  maxn:',i3,   &
          '  Iters:',i6/'Baud:',f7.3,'  BW:',f5.1,'  TxT:',f5.1,f5.2/)
  write(*,1004)
1004 format(/'  SNR     err    ber    fer   fsigma'/37('-'))

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

  call genmskhf(msgbits,id,icw,cbb,csync)!Generate baseband waveform and csync
  cb13=csync(1680:2095)                  !Copy the Barker 13 waveform
  a=0.
  a(1)=f0
  call twkfreq1(cbb,NZ,fs,a,cbb)         !Mix to specified frequency

  isna=-10
  isnb=-30
  if(snrdb.ne.0.0) then
     isna=nint(snrdb)
     isnb=isna
  endif
  do isnr=isna,isnb,-1                   !Loop over SNR range
     snrdb=isnr
     sig=sqrt(bandwidth_ratio) * 10.0**(0.05*snrdb)
     if(snrdb.gt.90.0) sig=1.0
     nhard=0
     nhardsync=0
     nfe=0
     sqf=0.
     do iter=1,iters                     !Loop over requested iterations
        c=cbb
        if(delay.ne.0.0 .or. fspread.ne.0.0) then
           call watterson(c,NZ,fs,delay,fspread)
        endif
        c=sig*c                          !Scale to requested SNR
        if(snrdb.lt.90) then
           do i=0,NZ-1                   !Generate gaussian noise
              xnoise(i)=gran()
              ynoise(i)=gran()
           enddo
           c=c + cmplx(xnoise,ynoise)    !Add AWGN noise
        endif

        call getfc1(c,fc1)               !First approx for freq
        call getfc2(c,csync,fc1,fc2,fc3) !Refined freq
        sqf=sqf + (fc1+fc2-f0)**2

!NB: Measured performance is about equally good using fc2 or fc3 here:
        a(1)=-(fc1+fc2)
        a(2:5)=0.
        call twkfreq1(c,NZ,fs,a,c)       !Mix c down by fc1+fc2

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

!---------------------------------------------------------------- DT
! Not presently used:
        amax=0.
        jpk=0
        do j=-20*NSPS,20*NSPS            !Get jpk
           z=sum(c(1680+j:2095+j)*cb13)/208.0
           if(abs(z).gt.amax) then
              amax=abs(z)
              jpk=j
           endif
        enddo
        xdt=jpk/fs

        nterms=maxn
        c0=c
        do itry=1,10
           idf=itry/2
           if(mod(itry,2).eq.0) idf=-idf
           nhard0=0
           nhardsync0=0
           ifer=1
           a(1)=idf*0.01
           a(2:5)=0.
           call twkfreq1(c0,NZ,fs,a,c)       !Mix c0 into c
           call cpolyfit(c,pp,id,maxn,aa,bb,zz,nhs)
           call msksoftsym(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhardsync0)
           if(nhardsync0.gt.12) cycle
           rxav=sum(rxdata)/ND
           rx2av=sum(rxdata*rxdata)/ND
           rxsig=sqrt(rx2av-rxav*rxav)
           rxdata=rxdata/rxsig
           ss=0.84
           llr=2.0*rxdata/(ss*ss)
           apmask=0
           max_iterations=40
           ifer=0
           call bpdecode168(llr,apmask,max_iterations,decoded,niterations,cw)
           nbadcrc=0
           if(niterations.ge.0) call chkcrc12(decoded,nbadcrc)
           if(niterations.lt.0 .or. count(msgbits.ne.decoded).gt.0 .or.        &
                nbadcrc.ne.0) ifer=1
!           if(ifer.eq.0) write(67,1301) snrdb,itry,idf,niterations,    &
!                nhardsync0,nhard0
!1301       format(f6.1,5i6)
           if(ifer.eq.0) exit
        enddo                                !Freq dither loop
        nhard=nhard+nhard0
        nhardsync=nharsdync+nhardsync0
        nfe=nfe+ifer
     enddo

     fsigma=sqrt(sqf/iters)
     ber=float(nhard)/((NS+ND)*iters)
     fer=float(nfe)/iters
     write(*,1050)  snrdb,nhard,ber,fer,fsigma
!     write(60,1050)  snrdb,nhard,ber,fer,fsigma
1050 format(f6.1,i7,f8.4,f7.3,f8.2)
  enddo

999 end program msksim


code:
subroutine msksoftsym(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhardsync0)

  parameter (KK=84)                     !Information bits (72 + CRC12)
  parameter (ND=168)                    !Data symbols: LDPC (168,84), r=1/2
  parameter (NS=65)                     !Sync symbols (2 x 26 + Barker 13)
  parameter (NR=3)                      !Ramp up/down
  parameter (NN=NR+NS+ND)               !Total symbols (236)
  parameter (NSPS=16)                   !Samples per MSK symbol (16)
  parameter (N2=2*NSPS)                 !Samples per OQPSK symbol (32)
  parameter (N13=13*N2)                 !Samples in central sync vector (416)
  parameter (NZ=NSPS*NN)                !Samples in baseband waveform (3760)
  parameter (NFFT1=4*NSPS,NH1=NFFT1/2)

  complex zz(NS+ND)                     !Complex symbol values (intermediate)
  complex z,z0
  real rxdata(ND)                       !Soft symbols
  real aa(20),bb(20)                    !Fitted polyco's
  integer id(NS+ND)                     !NRZ values (+/-1) for Sync and Data
  integer ierror(NS+ND)

  n=0
  ierror=0
  do j=1,117
     xx=j*2.0/117.0 - 1.0
     yii=1.
     yqq=0.
     if(nterms.gt.0) then
        yii=aa(1)
        yqq=bb(1)
        do i=2,nterms
           yii=yii + aa(i)*xx**(i-1)
           yqq=yqq + bb(i)*xx**(i-1)
        enddo
     endif
     z0=cmplx(yii,yqq)
     z=zz(j)*conjg(z0)
     p=real(z)
     if(abs(id(j)).eq.2) then
        if(real(z)*id(j).lt.0) then              !Sync bit
           nhardsync0=nhardsync0+1
           ierror(j)=2
        endif
     else
        n=n+1                                    !Data bit
        rxdata(n)=p
        ierr=0
        if(id(j)*p.lt.0) then
           ierr=1
           ierror(j)=1
        endif
        nhard0=nhard0+ierr
     endif
  enddo

  do j=118,233
     xx=(j-116.5)*2.0/117.0 - 1.0
     yii=1.
     yqq=0.
     if(nterms.gt.0) then
        yii=aa(1)
        yqq=bb(1)
        do i=2,nterms
           yii=yii + aa(i)*xx**(i-1)
           yqq=yqq + bb(i)*xx**(i-1)
        enddo
     endif
     z0=cmplx(yii,yqq)
     z=zz(j)*conjg(z0)
     p=aimag(z)
     n=n+1
     rxdata(n)=p
     ierr=0
     if(id(j)*p.lt.0) then
        ierr=1
        ierror(j)=1
     endif
     nhard0=nhard0+ierr
  enddo

  return
end subroutine msksoftsym


code:
subroutine msksoftsymw(zz,aa,bb,id,nterms,ierror,rxdata,nhard0,nhardsync0)

  include 'wsprlf_params.f90'

  complex zz(NS+ND)                     !Complex symbol values (intermediate)
  complex z,z0
  real rxdata(ND)                       !Soft symbols
  real aa(20),bb(20)                    !Fitted polyco's
  integer id(NS+ND)                     !NRZ values (+/-1) for Sync and Data
  integer ierror(NS+ND)

  n=0
  ierror=0
  ierr=0
  jz=(NS+ND+1)/2
  do j=1,jz
     xx=j*2.0/jz - 1.0
     yii=1.
     yqq=0.
     if(nterms.gt.0) then
        yii=aa(1)
        yqq=bb(1)
        do i=2,nterms
           yii=yii + aa(i)*xx**(i-1)
           yqq=yqq + bb(i)*xx**(i-1)
        enddo
     endif
     z0=cmplx(yii,yqq)
     z=zz(j)*conjg(z0)
     p=real(z)
     if(abs(id(j)).eq.2) then
        if(real(z)*id(j).lt.0) then              !Sync bit
           nhardsync0=nhardsync0+1
           ierror(j)=2
        endif
     else
        n=n+1                                    !Data bit
        rxdata(n)=p
        ierr=0
        if(id(j)*p.lt.0) then
           ierr=1
           ierror(j)=1
        endif
        nhard0=nhard0+ierr
     endif
!     write(41,3301) j,id(j),ierror(j),ierr,n,p,p*id(j)
!3301 format(5i6,2f10.3)
  enddo

  do j=jz+1,NS+ND
     xx=(j-jz+0.5)*2.0/jz - 1.0
     yii=1.
     yqq=0.
     if(nterms.gt.0) then
        yii=aa(1)
        yqq=bb(1)
        do i=2,nterms
           yii=yii + aa(i)*xx**(i-1)
           yqq=yqq + bb(i)*xx**(i-1)
        enddo
     endif
     z0=cmplx(yii,yqq)
     z=zz(j)*conjg(z0)
     p=aimag(z)
     n=n+1
     if(n.gt.ND) exit
     rxdata(n)=p
     ierr=0
     if(id(j)*p.lt.0) then
        ierr=1
        ierror(j)=1
     endif
     nhard0=nhard0+ierr
!     write(41,3301) j,id(j),ierror(j),ierr,n,p,p*id(j)
  enddo

  return
end subroutine msksoftsymw


code:
subroutine osd204(llr,apmask,ndeep,decoded,cw,nhardmin,dmin)
!
! An ordered-statistics decoder for the (204,68) code.
!
include "ldpc_204_68_params.f90"

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

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

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

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

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

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

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

g2=transpose(genmrb)

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

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

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

cw=c0
ntotal=0
nrejected=0

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


code:
subroutine osd300(llr,apmask,ndeep,decoded,cw,nhardmin,dmin)
!
! An ordered-statistics decoder for the (300,60) code.
!
include "ldpc_300_60_params.f90"

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

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

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

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

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

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

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

g2=transpose(genmrb)

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

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

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

cw=c0
ntotal=0
nrejected=0

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


code:
subroutine polyfit4(x,y,sigmay,npts,nterms,mode,a,chisqr)

  parameter (MAXN=20)
  implicit real*8 (a-h,o-z)
  real x(npts), y(npts), sigmay(npts), a(nterms),chisqr
  real*8 sumx(2*MAXN-1), sumy(MAXN), array(MAXN,MAXN)

! Accumulate weighted sums
  nmax = 2*nterms-1
  sumx=0.
  sumy=0.
  chisq=0.
  do i=1,npts
     xi=x(i)
     yi=y(i)
     if(mode.lt.0) then
        weight=1./abs(yi)
     else if(mode.eq.0) then
        weight=1
     else
        weight=1./sigmay(i)**2
     end if
     xterm=weight
     do n=1,nmax
        sumx(n)=sumx(n)+xterm
        xterm=xterm*xi
     enddo
     yterm=weight*yi
     do n=1,nterms
        sumy(n)=sumy(n)+yterm
        yterm=yterm*xi
     enddo
     chisq=chisq+weight*yi**2
  enddo

! Construct matrices and calculate coefficients
  do j=1,nterms
     do k=1,nterms
        n=j+k-1
        array(j,k)=sumx(n)
     enddo
  enddo

  delta=determ4(array,nterms)
  if(delta.eq.0) then
     chisqr=0.
     a=0.
  else
     do l=1,nterms
        do j=1,nterms
           do k=1,nterms
              n=j+k-1
              array(j,k)=sumx(n)
           enddo
           array(j,l)=sumy(j)
        enddo
        a(l)=determ4(array,nterms)/delta
     enddo

! Calculate chi square

     do j=1,nterms
        chisq=chisq-2*a(j)*sumy(j)
        do k=1,nterms
           n=j+k-1
           chisq=chisq+a(j)*a(k)*sumx(n)
        enddo
     enddo
     free=npts-nterms
     chisqr=chisq/free
  end if

  return
end subroutine polyfit4

real*8 function determ4(array,norder)

  parameter (MAXN=20)
  implicit real*8 (a-h,o-z)
  real*8 array(MAXN,MAXN)

  determ4=1.
  do k=1,norder
     if (array(k,k).ne.0) go to 41
     do j=k,norder
        if(array(k,j).ne.0) go to 31
     enddo
     determ4=0.
     go to 60

31   do i=k,norder
        s8=array(i,j)
        array(i,j)=array(i,k)
        array(i,k)=s8
     enddo
     determ4=-1.*determ4
41   determ4=determ4*array(k,k)
     if(k.lt.norder) then
        k1=k+1
        do i=k1,norder
           do j=k1,norder
              array(i,j)=array(i,j)-array(i,k)*array(k,j)/array(k,k)
           enddo
        enddo
     end if
  enddo

60 return
end function determ4


code:
subroutine spec4(c,s,savg)

  parameter (KK=84)                     !Information bits (72 + CRC12)
  parameter (ND=84)                     !Data symbols: LDPC (168,84), r=1/2
  parameter (NS=12)                     !Sync symbols (3 @ 4x4 Costas arrays)
  parameter (NR=2)                      !Ramp up/down
  parameter (NN=NR+NS+ND)               !Total symbols (98)
  parameter (NSPS=2688/84)              !Samples per symbol (32)
  parameter (NZ=NSPS*NN)                !Samples in baseband waveform (3760)
  parameter (NFFT=2*NSPS,NH=NSPS)

  complex c(0:NZ-1)
  complex c1(0:NFFT-1)
  real s(0:NH,NN)
  real savg(0:NH)

  fs=12000.0/84.0
  df=fs/NFFT
  savg=0.
  do j=1,NN
     ia=(j-1)*NSPS
     ib=ia + NSPS-1
     c1(0:NSPS-1)=c(ia:ib)
     c1(NSPS:)=0.
     call four2a(c1,NFFT,1,-1,1)
     do k=1,NSPS
        s(k,j)=real(c1(k))**2 + aimag(c1(k))**2
     enddo
     savg=savg+s(0:NH,j)
  enddo
  s=s/NZ
  savg=savg/(NN*NZ)

  return
end subroutine spec4


code:
subroutine spec8(c,s,savg)

  include 'wspr_fsk8_params.f90'
  complex c(0:NMAXD-1)
  complex c1(0:NSPS-1)
  real s(0:NH2,NN)
  real savg(0:NH2)

  fs=12000.0/NDOWN
  df=fs/NSPS
  savg=0.
  do j=1,NN
     ia=(j-1)*NSPS
     ib=ia + NSPS-1
     c1(0:NSPS-1)=c(ia:ib)
     c1(NSPS:)=0.
     call four2a(c1,NSPS,1,-1,1)
     do i=0,NH2
        s(i,j)=real(c1(i))**2 + aimag(c1(i))**2
     enddo
     savg=savg+s(0:NH2,j)
  enddo
  s=s/NZ
  savg=savg/(NN*NZ)
!  do i=0,NH2
!     write(31,3101) i*df,savg(i)
!3101 format(f10.3,f12.3)
!  enddo

  return
end subroutine spec8