WSJT-X/lib/fast9.f90

subroutine fast9(id2,narg,line)

! Decoder for "fast9" modes, JT9E to JT9H.

  parameter (NMAX=30*12000,NSAVE=500)
  integer*2 id2(0:NMAX)
  integer narg(0:14)
  integer*1 i1SoftSymbols(207)
  integer*1 i1save(207,NSAVE)
  integer indx(NSAVE)
  integer*8 count0,count1,clkfreq
  real s1(720000)                      !To reserve space.  Logically s1(nq,jz)
  real s2(240,340)                     !Symbol spectra at quarter-symbol steps
  real ss2(0:8,85)                     !Folded symbol spectra
  real ss3(0:7,69)                     !Folded spectra without sync symbols
  real s(1500)
  real ccfsave(NSAVE)
  real t0save(NSAVE)
  real t1save(NSAVE)
  real freqSave(NSAVE)
  real t(6)
  character*22 msg                     !Decoded message
  character*80 line(100)
  data nsubmode0/-1/,ntot/0/
  save s1,nsubmode0,ntot

! Parameters from GUI are in narg():
  nutc=narg(0)                         !UTC
  npts=min(narg(1),NMAX)               !Number of samples in id2 (12000 Hz)
  nsubmode=narg(2)                     !0=A 1=B 2=C 3=D 4=E 5=F 6=G 7=H
  if(nsubmode.lt.4) go to 900
  newdat=narg(3)                       !1==> new data, compute symbol spectra
  minsync=narg(4)                      !Lower sync limit
  npick=narg(5)
  t0=0.001*narg(6)
  t1=0.001*narg(7)
  maxlines=narg(8)                     !Max # of decodes to return to caller
  nmode=narg(9)
  nrxfreq=narg(10)                     !Targer Rx audio frequency (Hz)
  ntol=narg(11)                        !Search range, +/- ntol (Hz)

  tmid=npts*0.5/12000.0
  line(1:100)(1:1)=char(0)
  s=0
  s2=0
  nsps=60 * 2**(7-nsubmode)            !Samples per sysbol
  nfft=2*nsps                          !FFT size
  nh=nfft/2
  nq=nfft/4
  istep=nsps/4                         !Symbol spectra at quarter-symbol steps
  jz=npts/istep
  df=12000.0/nfft                      !FFT bin width
  db1=db(2500.0/df)
  nfa=max(200,nrxfreq-ntol)            !Lower frequency limit
  nfb=min(nrxfreq+ntol,2500)           !Upper frequency limit
  nline=0
  t=0.

  if(newdat.eq.1 .or. nsubmode.ne.nsubmode0) then
     call system_clock(count0,clkfreq)
     call spec9f(id2,npts,nsps,s1,jz,nq)          !Compute symbol spectra, s1 
     call system_clock(count1,clkfreq)
     t(1)=t(1)+float(count1-count0)/float(clkfreq)
  endif

  nsubmode0=nsubmode
  tmsg=nsps*85.0/12000.0
  limit=2000
  nlen0=0
  i1=0
  i2=0
  ccfsave=0.
  do ilength=1,14
     nlen=1.4142136**(ilength-1)
     if(nlen.gt.jz/340) nlen=jz/340
     if(nlen.eq.nlen0) cycle
     nlen0=nlen

     db0=db(float(nlen))
     jlen=nlen*340
     jstep=jlen/4                      !### Is this about right? ###
     if(nsubmode.ge.6) jstep=jlen/2

     do ja=1,jz-jlen,jstep
        jb=ja+jlen-1
        call system_clock(count0,clkfreq)
        call foldspec9f(s1,nq,jz,ja,jb,s2)        !Fold symbol spectra into s2
        call system_clock(count1,clkfreq)
        t(2)=t(2)+float(count1-count0)/float(clkfreq)

! Find sync; put sync'ed symbol spectra into ss2 and ss3
! Might want to do a peakup in DT and DF, then re-compute symbol spectra.

        call system_clock(count0,clkfreq)
        call sync9f(s2,nq,nfa,nfb,ss2,ss3,lagpk,ipk,ccfbest)
        call system_clock(count1,clkfreq)
        t(3)=t(3)+float(count1-count0)/float(clkfreq)

        i1=i1+1
        if(ccfbest.lt.30.0) cycle
        call system_clock(count0,clkfreq)
        call softsym9f(ss2,ss3,i1SoftSymbols)     !Compute soft symbols
        call system_clock(count1,clkfreq)
        t(4)=t(4)+float(count1-count0)/float(clkfreq)

        i2=i2+1
        ccfsave(i2)=ccfbest
        i1save(1:207,i2)=i1SoftSymbols
        t0=(ja-1)*istep/12000.0
        t1=(jb-1)*istep/12000.0
        t0save(i2)=t0
        t1save(i2)=t1
        freq=ipk*df
        freqSave(i2)=freq
     enddo
  enddo
  nsaved=i2

  ccfsave(1:nsaved)=-ccfsave(1:nsaved)
  call system_clock(count0,clkfreq)
  indx=0
  call indexx(ccfsave,nsaved,indx)
  call system_clock(count1,clkfreq)
  t(5)=t(5)+float(count1-count0)/float(clkfreq)

  ccfsave(1:nsaved)=-ccfsave(1:nsaved)

  do iter=1,2
!     do isave=1,nsaved
     do isave=1,50
        i2=indx(isave)
        if(i2.lt.1 .or. i2.gt.nsaved) cycle      !### Why needed? ###
        t0=t0save(i2)
        t1=t1save(i2)
        if(iter.eq.1 .and. t1.lt.tmid) cycle
        if(iter.eq.2 .and. t1.ge.tmid) cycle
        ccfbest=ccfsave(i2)
        i1SoftSymbols=i1save(1:207,i2)
        freq=freqSave(i2)
        call system_clock(count0,clkfreq)
        call jt9fano(i1SoftSymbols,limit,nlim,msg)      !Invoke Fano decoder
        call system_clock(count1,clkfreq)
        t(6)=t(6)+float(count1-count0)/float(clkfreq)

        i=t0*12000.0
        kz=(t1-t0)/0.02
        smax=0.
        do k=1,kz
           sq=0.
           do n=1,240
              i=i+1
              x=id2(i)
              sq=sq+x*x
           enddo
           s(k)=sq/240.
           smax=max(s(k),smax)
        enddo
        call pctile(s,kz,35,base)
        snr=smax/(1.1*base) - 1.0
        nsnr=-20
        if(snr.gt.0.0) nsnr=nint(db(snr))

!        write(72,3002) nutc,iter,isave,nlen,tmid,t0,t1,ccfbest,   &
!             nint(freq),nlim,msg
!3002    format(i6.6,i1,i4,i3,4f6.1,i5,i7,1x,a22)

        if(msg.ne.'                      ') then

! Display multiple decodes only if they differ:
           do n=1,nline
              if(index(line(n),msg).gt.1) go to 100
           enddo
!### Might want to use decoded message to get a complete estimate of S/N.
           nline=nline+1
           write(line(nline),1000) nutc,nsnr,t0,nint(freq),msg,char(0)
1000       format(i6.6,i4,f5.1,i5,1x,'@',1x,a22,a1)
           ntot=ntot+1
!           write(70,5001) nsaved,isave,nline,maxlines,ntot,nutc,msg
!5001       format(5i5,i7.6,1x,a22)
           if(nline.ge.maxlines) go to 900
        endif
100     continue
     enddo
  enddo

900 continue
!  write(*,6001) t,t(6)/sum(t)
!6001 format(7f10.3)

  return
end subroutine fast9