code:
subroutine sfrsd(mrsym,mrprob,mr2sym,mr2prob,ntrials,correct,indexes,   &
     param,ntry)

  integer mrsym(0:62),mrprob(0:62),mr2sym(0:62),mr2prob(0:62)
  integer correct(0:62),indexes(0:62),probs(0:62),thresh0(0:62)
  integer rxdat(0:62),rxdat2(0:62),rxprob(0:62),rxprob2(0:62)
  integer workdat(0:62),era_pos(0:50)
  integer perr(0:7,0:7)
  integer param(0:7)
  real ratio0(0:62)

  call init_rs_int()
  do i=0,62
     rxdat(i)=mrsym(62-i)
     rxdat2(i)=mr2sym(62-i)
     rxprob(i)=mrprob(62-i)
     rxprob2(i)=mr2prob(62-i)
  enddo

  do i=0,62
     indexes(i)=i
     probs(i)=rxprob(i)
  enddo

  do ip=1,62
     do k=0,63-ip
        if(probs(k).lt.probs(k+1)) then
           ntmp=probs(k)
           probs(k)=probs(k+1)
           probs(k+1)=ntmp
           ntmp=indexes(k)
           indexes(k)=indexes(k+1)
           indexes(k+1)=ntmp
        endif
     enddo
  enddo

  era_pos=0
  numera=0
  workdat=rxdat
  call decode_rs_int()
  if(nerr.ge.0) then
     correct=workdat
     param=0
     return
  endif

  call random_seed()

  ncandidates=0
  nsum=0
  do i=0,62
     nsum=nsum+rxprob(i)
     j=indexes(62-i)
     ratio0(i)=float(rxprob2(j))/(float(rxprob(j))+0.01)
     ii=int(7.999*ratio0(i))
     jj=(62-i)/8
     thresh0(i)=nint(1.3*perr(jj,ii))
  enddo
  if(nsum.eq.0) return

  do k=0,ntrials
     era_pos=0
     workdat=rxdat
     numera=0
     do i=0,62
        j=indexes(62-i)
        thresh=thresh0(i)
        ir=rand()
        if(...) then
           era_pos(numera)=j
           numera=numera+1
        endif
     enddo

     call decode_rs_int()
     if(nerr.ge.0) then
        ncandidates=ncandidates+1
        nhard=0
        nsoft=0
        nsofter=0
        do i=0,62
           if(workdat(i).ne.rxdat(i)) then
              nhard=nhard+1
              nsofter=nsofter+rxprob(i)
              if(workdat(i).ne.rxdat2(i)) nsoft=nsoft+rxprob(i)
           else
              nsofter=nsofter-rxprob(i)
           endif
        enddo
        nsoft=63*nsoft/nsum
        nsofter=63*nsofter/nsum
        ntotal=nsoft+nhard
        if(ntotal.lt.ntotal_min) then
           nsoft_min=nsoft
           nhard_min=nhard
           nsofter_min=nsofter
           ntotal_min=ntotal
           correct=workdat
           nera_best=numera
           ntry=k
        endif
        if(ntotal_min.lt.72 .and. nhard_min.lt.42) exit
     endif
     if(k.eq.ntrials-1) ntry=k+1
  enddo

  if(ntotal_min.ge.76 .or. nhard.ge.44) nhard_min=-1

  param(0)=ncandidates
  param(1)=nhard_min
  param(2)=nsoft_min
  param(3)=nera_best
  param(4)=nsofter_min
  if(param(0).eq.0) param(2)=-1

  return
end subroutine sfrsd


code:
#include "init_random_seed.h"

/* Fortran wrapper to seed the C library rand */
void sgran_(void)
{
  init_random_seed();
}

license: mit


code:
subroutine sh65(cx,n5,mode65,ntol,xdf,nspecial,snrdb)
  parameter(NFFT=2048,NH=NFFT/2,MAXSTEPS=150)
  complex cx(90000)
  complex c(0:NFFT-1)
  real s(-NH+1:NH)
  real s2(-NH+1:NH,MAXSTEPS)
  real ss(-NH+1:NH,8)
  real sigmax(8)
  integer ipk(8)

  s=0.
  ss=0.

  jstep=NFFT/4
  nblks=n5/jstep - 3
  ia=-jstep+1
  do iblk=1,nblks
     ia=ia+jstep
     ib=ia+NFFT-1
     c=cx(ia:ib)
     call four2a(c,nfft,1,1,1)            !c2c FFT
     do i=0,NFFT-1
        j=i
        if(j.gt.NH) j=j-NFFT
        p=real(c(i))**2 + aimag(c(i))**2
        s(j)=s(j) + p
        s2(j,iblk)=p
     enddo
     n=mod(iblk-1,8) +1
     ss(-NH+1:NH,n)=ss(-NH+1:NH,n) + s2(-NH+1:NH,iblk)
  enddo

  s=1.e-6*s
  ss=1.e-6*ss
  df=1378.1285/NFFT
  nfac=40*mode65
  dtstep=0.25/df

! Define range of frequencies to be searched
  fa=-ntol
  fb=ntol
  ia2=max(-NH+1,nint(fa/df))
! Upper tone is above sync tone by 4*nfac*df Hz
  ib2=min(NH,nint(fb/df + 4.1*nfac))

! Find strongest line in each of the 4 phases, repeating for each drift rate.
  sbest=0.
  snrbest=0.
  nbest=1
  ipk=0

  do n=1,8
     sigmax(n)=0.
     do i=ia2,ib2
        sig=ss(i,n)
        if(sig.ge.sigmax(n)) then
           ipk(n)=i
           sigmax(n)=sig
           if(sig.ge.sbest) then
              sbest=sig
              nbest=n
           endif
        endif
     enddo
  enddo
  n2best=nbest+4
  if(n2best.gt.8) n2best=nbest-4
  xdf=min(ipk(nbest),ipk(n2best))*df
  nspecial=0
  if(abs(xdf).le.ntol) then
     idiff=abs(ipk(nbest)-ipk(n2best))
     xk=float(idiff)/nfac
     k=nint(xk)
     iderr=nint((xk-k)*nfac)
!     maxerr=nint(0.008*abs(idiff) + 0.51)
     maxerr=nint(0.02*abs(idiff) + 0.51)     !### Better test ??? ###
     if(abs(iderr).le.maxerr .and. k.ge.2 .and. k.le.4) nspecial=k
     snrdb=-30.0
     if(nspecial.gt.0) then
        call sh65snr(ss(ia2,nbest),ib2-ia2+1,snr1)
        call sh65snr(ss(ia2,n2best),ib2-ia2+1,snr2)
        snr=0.5*(snr1+snr2)
        snrdb=db(snr) - db(2500.0/df) - db(sqrt(nblks/4.0)) + 8.0
     endif
     if(snr1.lt.4.0 .or. snr2.lt.4.0 .or. snr.lt.5.0) nspecial=0
  endif

  return
end subroutine sh65


code:
subroutine sh65snr(x,nz,snr)

  real x(nz)

  ipk=0 !Shut up compiler warnings. -db
  smax=-1.e30
  do i=1,nz
     if(x(i).gt.smax) then
        ipk=i
        smax=x(i)
     endif
     s=s+x(i)
  enddo

  s=0.
  ns=0
  do i=1,nz
     if(abs(i-ipk).ge.3) then
        s=s+x(i)
        ns=ns+1
     endif
  enddo
  ave=s/ns

  sq=0.
  do i=1,nz
     if(abs(i-ipk).ge.3) then
        sq=sq+(x(i)-ave)**2
        ns=ns+1
     endif
  enddo
  rms=sqrt(sq/(nz-2))
  snr=(smax-ave)/rms

  return
end subroutine sh65snr


code:
subroutine shell(n,a)
  integer n
  real a(n)
  integer i,j,inc
  real v

  inc=1
1 inc=3*inc+1
  if(inc.le.n) go to 1
2  inc=inc/3

  do i=inc+1,n
     v=a(i)
     j=i
3    if(a(j-inc).gt.v) then
        a(j)=a(j-inc)
        j=j-inc
        if(j.le.inc) go to 4
        go to 3
     endif
4    a(j)=v
  enddo

  if(inc.gt.1) go to 2

  return
end subroutine shell