code:
subroutine to_contest_msg(msg0,msg)

! If the message has "R grid4" istead of "grid4", remove the "R "
! and substitute the diametrically opposite grid.

  character*6 g1,g2
  character*22 msg0,msg
  logical isgrid
  isgrid(g1)=g1(1:1).ge.'A' .and. g1(1:1).le.'R' .and. g1(2:2).ge.'A' .and. &
       g1(2:2).le.'R' .and. g1(3:3).ge.'0' .and. g1(3:3).le.'9' .and.       &
       g1(4:4).ge.'0' .and. g1(4:4).le.'9' .and. g1(1:4).ne.'RR73'

  i0=index(msg0,' R ') + 3                  !Check for ' R ' in message
  g1=msg0(i0:i0+3)//'  '
  if(isgrid(g1)) then                       !Check for ' R grid'
     call grid2deg(g1,dlong,dlat)
     dlong=dlong+180.0
     if(dlong.gt.180.0) dlong=dlong-360.0
     dlat=-dlat
     call deg2grid(dlong,dlat,g2)           !g2=antipodes grid
     msg=msg0(1:i0-3)//g2(1:4)              !Send message with g2
  else
     msg=msg0
  endif

  return
end subroutine to_contest_msg


code:
#include "config.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "tstrig.h"

int set_conf(RIG *my_rig, char *conf_parms);

int rig_control(rig_model_t my_model, int verbose);

int main (int argc, char *argv[])
{
  rig_model_t my_model = RIG_MODEL_DUMMY;
  int verbose=0;

  my_model=214;
  rig_control(my_model,verbose);
  return 0;
}

/*
gcc -c -Wall -I../include rig_control.c
gcc -c -Wall -I../include tstrig.c
gcc -o tstrig.exe -Wl,--enable-auto-import tstrig.o rig_control.o libhamlib.dll.a
strip tstrig.exe
*/

license: mit


code:
/*
 * rigctl_parse.h - (C) Stephane Fillod 2000-2010
 *
 * This program test/control a radio using Hamlib.
 * It takes commands in interactive mode as well as
 * from command line options.
 *
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License along
 *   with this program; if not, write to the Free Software Foundation, Inc.,
 *   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 *
 */

#ifndef RIGCTL_PARSE_H
#define RIGCTL_PARSE_H

#include <stdio.h>
#include <hamlib/rig.h>

/*
 * external prototype
 */

int dumpcaps (RIG *, FILE *);
int dumpconf (RIG *, FILE *);

/*
 * Prototypes
 */
void usage_rig(FILE *);
void version();
void list_models();
int dump_chan(FILE *, RIG*, channel_t*);
int print_conf_list(const struct confparams *cfp, rig_ptr_t data);
int set_conf(RIG *my_rig, char *conf_parms);

int rigctl_parse(RIG *my_rig, FILE *fin, FILE *fout, char *argv[], int argc);

#endif	/* RIGCTL_PARSE_H */


code:
subroutine tweak1(ca,jz,f0,cb)

! Shift frequency of analytic signal ca, with output to cb

  complex ca(jz),cb(jz)
  real*8 twopi
  complex*16 w,wstep
  complex w4
  data twopi/0.d0/
  save twopi

  if(twopi.eq.0.d0) twopi=8.d0*atan(1.d0)
  w=1.d0
  dphi=twopi*f0/12000.d0
  wstep=cmplx(cos(dphi),sin(dphi))
  do i=1,jz
     w=w*wstep
     w4=w
     cb(i)=w4*ca(i)
  enddo

  return
end subroutine tweak1


code:
subroutine twkfreq(c3,c4,npts,fsample,a)

  complex c3(npts)
  complex c4(npts)
  complex w,wstep
  real a(3)
  data twopi/6.283185307/

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

  return
end subroutine twkfreq


code:
subroutine twkfreq65(c4aa,n5,a)

  complex c4aa(n5)
  real a(5)
  complex w,wstep
  data twopi/6.283185307/

! Apply AFC corrections to the c4aa data
  w=1.0
  wstep=1.0
  x0=0.5*(n5+1)
  s=2.0/n5
  do i=1,n5
     x=s*(i-x0)
     if(mod(i,100).eq.1) then
        p2=1.5*x*x - 0.5
        dphi=(a(1) + x*a(2) + p2*a(3)) * (twopi/1378.125)
        wstep=cmplx(cos(dphi),sin(dphi))
     endif
     w=w*wstep
     c4aa(i)=w*c4aa(i)
  enddo

  return
end subroutine twkfreq65


code:
 subroutine unpackmsg144(dat,msg,c1,c2)
! special unpackmsg for MSK144 - returns call1 and call2 to enable
! maintenance of a recent-calls-heard list

   use packjt
   parameter (NBASE=37*36*10*27*27*27)
   parameter (NGBASE=180*180)
   integer dat(12)
   character c1*12,c2*12,grid*4,msg*22,grid6*6,psfx*4,junk2*4
   logical cqnnn

   cqnnn=.false.
   nc1=ishft(dat(1),22) + ishft(dat(2),16) + ishft(dat(3),10)+         &
        ishft(dat(4),4) + iand(ishft(dat(5),-2),15)

   nc2=ishft(iand(dat(5),3),26) + ishft(dat(6),20) +                   &
        ishft(dat(7),14) + ishft(dat(8),8) + ishft(dat(9),2) +         &
        iand(ishft(dat(10),-4),3)

   ng=ishft(iand(dat(10),15),12) + ishft(dat(11),6) + dat(12)

   if(ng.ge.32768) then
      call unpacktext(nc1,nc2,ng,msg)
      c1(1:12)=' '
      c2(1:12)=' '
      go to 100
   endif

   call unpackcall(nc1,c1,iv2,psfx)
   if(iv2.eq.0) then
 ! This is an "original JT65" message
      if(nc1.eq.NBASE+1) c1='CQ    '
      if(nc1.eq.NBASE+2) c1='QRZ   '
      nfreq=nc1-NBASE-3
      if(nfreq.ge.0 .and. nfreq.le.999) then
         write(c1,1002) nfreq
 1002    format('CQ ',i3.3)
         cqnnn=.true.
      endif
   endif

   call unpackcall(nc2,c2,junk1,junk2)
   call unpackgrid(ng,grid)

   if(iv2.gt.0) then
 ! This is a JT65v2 message
      do i=1,4
         if(ichar(psfx(i:i)).eq.0) psfx(i:i)=' '
      enddo

      n1=len_trim(psfx)
      n2=len_trim(c2)
      if(iv2.eq.1) msg='CQ '//psfx(:n1)//'/'//c2(:n2)//' '//grid
      if(iv2.eq.2) msg='QRZ '//psfx(:n1)//'/'//c2(:n2)//' '//grid
      if(iv2.eq.3) msg='DE '//psfx(:n1)//'/'//c2(:n2)//' '//grid
      if(iv2.eq.4) msg='CQ '//c2(:n2)//'/'//psfx(:n1)//' '//grid
      if(iv2.eq.5) msg='QRZ '//c2(:n2)//'/'//psfx(:n1)//' '//grid
      if(iv2.eq.6) msg='DE '//c2(:n2)//'/'//psfx(:n1)//' '//grid
      if(iv2.eq.7) msg='DE '//c2(:n2)//' '//grid
      if(iv2.eq.8) msg=' '
      go to 100
   else

   endif

   grid6=grid//'ma'
   call grid2k(grid6,k)
   if(k.ge.1 .and. k.le.450)   call getpfx2(k,c1)
   if(k.ge.451 .and. k.le.900) call getpfx2(k,c2)

   i=index(c1,char(0))
   if(i.ge.3) c1=c1(1:i-1)//'            '
   i=index(c2,char(0))
   if(i.ge.3) c2=c2(1:i-1)//'            '

   msg='                      '
   j=0
   if(cqnnn) then
      msg=c1//'          '
      j=7                                  !### ??? ###
      go to 10
   endif

   do i=1,12
      j=j+1
      msg(j:j)=c1(i:i)
      if(c1(i:i).eq.' ') go to 10
   enddo
   j=j+1
   msg(j:j)=' '

 10 do i=1,12
      if(j.le.21) j=j+1
      msg(j:j)=c2(i:i)
      if(c2(i:i).eq.' ') go to 20
   enddo
   if(j.le.21) j=j+1
   msg(j:j)=' '

 20 if(k.eq.0) then
      do i=1,4
         if(j.le.21) j=j+1
         msg(j:j)=grid(i:i)
      enddo
      if(j.le.21) j=j+1
      msg(j:j)=' '
   endif

 100 continue
   if(msg(1:6).eq.'CQ9DX ') msg(3:3)=' '
   if(msg(1:2).eq.'E9' .and.                                          &
        msg(3:3).ge.'A' .and. msg(3:3).le.'Z' .and.                   &
        msg(4:4).ge.'A' .and. msg(4:4).le.'Z' .and.                   &
        msg(5:5).eq.' ') msg='CQ '//msg(3:)

   return
 end subroutine unpackmsg144


code:
subroutine update_hasharray(recent_calls,nrecent,nhasharray)

  character*12 recent_calls(nrecent)
  character*22 hashmsg
  integer nhasharray(nrecent,nrecent)

  nhasharray=-1
  do i=1,nrecent
    do j=i+1,nrecent
      if( recent_calls(i)(1:1) .ne. ' ' .and. recent_calls(j)(1:1) .ne. ' ' ) then
        hashmsg=trim(recent_calls(i))//' '//trim(recent_calls(j))
        call fmtmsg(hashmsg,iz)
        call hash(hashmsg,22,ihash)
        ihash=iand(ihash,4095)
        nhasharray(i,j)=ihash
        hashmsg=trim(recent_calls(j))//' '//trim(recent_calls(i))
        call fmtmsg(hashmsg,iz)
        call hash(hashmsg,22,ihash)
        ihash=iand(ihash,4095)
        nhasharray(j,i)=ihash
      endif
    enddo
  enddo

end subroutine update_hasharray


code:
subroutine update_recent_calls(call,calls_hrd,nsize)
character*12 call,calls_hrd(nsize)

 new=1
 do ic=1,nsize
   if( calls_hrd(ic).eq.call ) then
     new=0
   endif
 enddo

 if( new.eq.1 ) then
   do ic=nsize-1,1,-1
     calls_hrd(ic+1)(1:12)=calls_hrd(ic)(1:12)
   enddo
   calls_hrd(1)(1:12)=call(1:12)
 endif

 return
 end subroutine update_recent_calls


code:
#include <unistd.h>
#include "sleep.h"

/* usleep(3) */
void usleep_(unsigned long *microsec)
{
  usleep(*microsec);
}


code:
/* Viterbi decoder for arbitrary convolutional code
 * viterbi27 and viterbi37 for the r=1/2 and r=1/3 K=7 codes are faster
 * Copyright 1999 Phil Karn, KA9Q
 * Modifications by Joe Taylor, K1JT
 * May be used under the terms of the GNU Public License
 */

#include <stdlib.h>

/* Select code here */

#define V213

#ifdef V213
#define	K 13			/* Constraint length */
#define N 2			/* Number of symbols per data bit */
#define Polys	Poly213		/* Select polynomials here */
#endif

/* Rate 1/2 codes */
unsigned int Poly213[] = {012767,016461};	/* k = 13  */

#include <memory.h>

#define LONGBITS 32
#define LOGLONGBITS 5

#undef max
#define max(x,y) ((x) > (y) ? (x) : (y))
#define D       (1 << max(0,K-LOGLONGBITS-1))
#define MAXNBITS 200            /* Maximum frame size (user bits) */

extern unsigned char Partab[];	/* Parity lookup table */

int Syms[1 << K];

int parity(int x)
{
  x ^= (x >> 16);
  x ^= (x >> 8);
  return Partab[x & 0xff];
}

/* Convolutionally encode data into binary symbols */
int enc213(unsigned char symbols[], unsigned char data[],
       unsigned int nbytes, unsigned int startstate,
       unsigned int endstate)
{
  unsigned int i,j,k;
  int l,n=-1;
  unsigned int encstate = startstate;

  for(k=0; k<nbytes; k++) {
    for(l=7;l>=0;l--){
      encstate = (encstate + encstate) + ((data[k] >> l) & 1);
      for(j=0;j<N;j++) {
	n=n+1;
	symbols[n] = parity(encstate & Polys[j]);
      }
    }
  }
  // Flush out with zero tail.  (No need, if tail-biting code.)
  for(i=0; i<K-1;i++){
    encstate = (encstate << 1) | ((endstate >> i) & 1);
    for(j=0;j<N;j++) {
      n=n+1;
      symbols[n] = parity(encstate & Polys[j]);
    }
  }
  return 0;
}

/* Viterbi decoder */
int vit213(
int *metric,           /* Final path metric (returned value) */
unsigned char *data,	/* Decoded output data */
unsigned char *symbols,	/* Raw deinterleaved input symbols */
unsigned int nbits,	/* Number of output bits */
int mettab[2][256],	/* Metric table, [sent sym][rx symbol] */
unsigned int startstate,         /* Encoder starting state */
unsigned int endstate            /* Encoder ending state */
){
  int bitcnt = -(K-1);
  int m0,m1;
  int i,j,sym,ipp;
  int mets[1 << N];
  unsigned int paths[(MAXNBITS+K-1)*D];
  unsigned int *pp,mask;
  int cmetric[1 << (K-1)],nmetric[1 << (K-1)];
  static int VDInit = 0;

  memset(paths,0,sizeof(paths));

  // Initialize on first time through:
  if(!VDInit){
    for(i=0;i<(1<<K);i++){
      sym = 0;
      for(j=0;j<N;j++)
	sym = (sym << 1) + parity(i & Polys[j]);
      Syms[i] = sym;
    }
    VDInit++;
  }

  // Keep only lower K-1 bits of specified startstate and endstate
  startstate &= ~((1<<(K-1)) - 1);
  endstate &=   ~((1<<(K-1)) - 1);

  /* Initialize starting metrics */
  for(i=0;i< 1<<(K-1);i++)
    cmetric[i] = -999999;
  cmetric[startstate] = 0;

  pp = paths;
  ipp=0;
  for(;;){ /* For each data bit */
    /* Read input symbols and compute branch metrics */
    for(i=0;i< 1<<N;i++){
      mets[i] = 0;
      for(j=0;j<N;j++){
	mets[i] += mettab[(i >> (N-j-1)) & 1][symbols[j]];
      }
    }
    symbols += N;
    /* Run the add-compare-select operations */
    mask = 1;
    for(i=0;i< 1 << (K-1);i+=2){
      int b1,b2;

      b1 = mets[Syms[i]];
      nmetric[i] = m0 = cmetric[i/2] + b1;
      b2 = mets[Syms[i+1]];
      b1 -= b2;
      m1 = cmetric[(i/2) + (1<<(K-2))] + b2;

      if(m1 > m0){
	nmetric[i] = m1;
	*pp |= mask;
      }

      m0 -= b1;
      nmetric[i+1] = m0;
      m1 += b1;

      if(m1 > m0){
	nmetric[i+1] = m1;
	*pp |= mask << 1;
      }

      mask <<= 2;
      if(mask == 0){
	mask = 1;
	pp++;
	ipp++;
      }
    }
    if(mask != 1){
      pp++;
      ipp++;
    }
    if(++bitcnt == (int)nbits){
      *metric = nmetric[endstate];
      break;
    }
    memcpy(cmetric,nmetric,sizeof(cmetric));
  }

  /* Chain back from terminal state to produce decoded data */
  if(data == NULL)
    return 0;/* Discard output */
  memset(data,0,(nbits+7)/8); /* round up in case nbits % 8 != 0 */

  for(i=nbits-1;i >= 0;i--){
    //    int a0,a1;
    pp -= D;
    ipp -= D;
    m0=endstate >> LOGLONGBITS;
    m1=1L << (endstate & (LONGBITS-1));
    if(pp[m0] & m1) {
      //      a0=nmetric[endstate];
      endstate |= (1 << (K-1));
      //      a1=nmetric[endstate];
      data[i>>3] |= 0x80 >> (i&7);
      //      printf("B  %d  %d  %d  %d\n",*metric,i,a0,a1);
    }
    endstate >>= 1;
  }
  return 0;
}

// Wrapper for calling "encode" from Fortran:
void enc213_(
unsigned char data[],           // User data, 8 bits per byte
int *nbits,                     // Number of user bits
unsigned char symbols[],        // Encoded one-bit symbols, 8 per byte
int *nsymbols,                  // Number of symbols
int *kk,                        // K
int *nn)                        // N
{
  int nbytes;
  nbytes=(*nbits+7)/8;          // Always encode multiple of 8 information bits
  enc213(symbols,data,nbytes,0,0); // Do the encoding
  *nsymbols=(*nbits+K-1)*N;        // Return number of encoded symbols
  *kk=K;
  *nn=N;
}

// Wrapper for calling "viterbi" from Fortran:
void vit213_(
unsigned char symbols[],  /* Raw deinterleaved input symbols */
unsigned int *Nbits,	  /* Number of decoded information bits */
int mettab[2][256],	  /* Metric table, [sent sym][rx symbol] */
unsigned char ddec[],	  /* Decoded output data */
int *Metric              /* Final path metric (bigger is better) */
){
  int metric;
  vit213(&metric,ddec,symbols,*Nbits,mettab,0,0);
  *Metric=metric;
}


code:
/* Viterbi decoder for arbitrary convolutional code
 * viterbi27 and viterbi37 for the r=1/2 and r=1/3 K=7 codes are faster
 * Copyright 1999 Phil Karn, KA9Q
 * May be used under the terms of the GNU Public License
 */

/* Select code here */

#define V216

#ifdef V216
#define	K 16			/* Constraint length */
#define N 2			/* Number of symbols per data bit */
#define Polys	Poly216		/* Select polynomials here */
#endif

/* Rate 1/2 codes */
unsigned int Poly216[] = {0126723, 0152711};	/* k = 16  */

#include <stdlib.h>
#include <memory.h>

#define LONGBITS 32
#define LOGLONGBITS 5

#undef max
#define max(x,y) ((x) > (y) ? (x) : (y))
#define D       (1 << max(0,K-LOGLONGBITS-1))
#define MAXNBITS 200            /* Maximum frame size (user bits) */

extern unsigned char Partab[];	/* Parity lookup table */

int Syms[1 << K];
int VDInit = 0;

int parity(int x)
{
  x ^= (x >> 16);
  x ^= (x >> 8);
  return Partab[x & 0xff];
}

// Wrapper for calling "encode" from Fortran:
//void __stdcall ENCODE(
void enc216_(
unsigned char data[],           // User data, 8 bits per byte
int *nbits,                     // Number of user bits
unsigned char symbols[],        // Encoded one-bit symbols, 8 per byte
int *nsymbols,                  // Number of symbols
int *kk,                        // K
int *nn)                        // N
{
  int nbytes;
  nbytes=(*nbits+7)/8;          // Always encode multiple of 8 information bits
  enc216(symbols,data,nbytes,0,0); // Do the encoding
  *nsymbols=(*nbits+K-1)*N;        // Return number of encoded symbols
  *kk=K;
  *nn=N;
}

/* Convolutionally encode data into binary symbols */
  enc216(unsigned char symbols[], unsigned char data[],
       unsigned int nbytes, unsigned int startstate,
       unsigned int endstate)
{
  int i,j,k,n=-1;
  unsigned int encstate = startstate;

  for(k=0; k<nbytes; k++) {
    for(i=7;i>=0;i--){
      encstate = (encstate + encstate) + ((data[k] >> i) & 1);
      for(j=0;j<N;j++) {
	n=n+1;
	symbols[n] = parity(encstate & Polys[j]);
      }
    }
  }
  // Flush out with zero tail.  (No need, if tail-biting code.)
  for(i=0; i<K-1;i++){
    encstate = (encstate << 1) | ((endstate >> i) & 1);
    for(j=0;j<N;j++) {
      n=n+1;
      symbols[n] = parity(encstate & Polys[j]);
    }
  }
  return 0;
}

// Wrapper for calling "viterbi" from Fortran:
//void __stdcall VITERBI(
void vit216_(
unsigned char symbols[],  /* Raw deinterleaved input symbols */
unsigned int *Nbits,	  /* Number of decoded information bits */
int mettab[2][256],	  /* Metric table, [sent sym][rx symbol] */
unsigned char ddec[],	  /* Decoded output data */
long *Metric              /* Final path metric (bigger is better) */
){
  long metric;
  vit216(&metric,ddec,symbols,*Nbits,mettab,0,0);
  *Metric=metric;
}

/* Viterbi decoder */
int vit216(
long *metric,           /* Final path metric (returned value) */
unsigned char *data,	/* Decoded output data */
unsigned char *symbols,	/* Raw deinterleaved input symbols */
unsigned int nbits,	/* Number of output bits */
int mettab[2][256],	/* Metric table, [sent sym][rx symbol] */
unsigned int startstate,         /* Encoder starting state */
unsigned int endstate            /* Encoder ending state */
){
  int bitcnt = -(K-1);
  long m0,m1;
  int i,j,sym,ipp;
  int mets[1 << N];
  unsigned long paths[(MAXNBITS+K-1)*D];
  unsigned long *pp,mask;
  long cmetric[1 << (K-1)],nmetric[1 << (K-1)];

  memset(paths,0,sizeof(paths));

  // Initialize on first time through:
  if(!VDInit){
    for(i=0;i<(1<<K);i++){
      sym = 0;
      for(j=0;j<N;j++)
	sym = (sym << 1) + parity(i & Polys[j]);
      Syms[i] = sym;
    }
    VDInit++;
  }

  // Keep only lower K-1 bits of specified startstate and endstate
  startstate &= ~((1<<(K-1)) - 1);
  endstate &=   ~((1<<(K-1)) - 1);

  /* Initialize starting metrics */
  for(i=0;i< 1<<(K-1);i++)
    cmetric[i] = -999999;
  cmetric[startstate] = 0;

  pp = paths;
  ipp=0;
  for(;;){ /* For each data bit */
    /* Read input symbols and compute branch metrics */
    for(i=0;i< 1<<N;i++){
      mets[i] = 0;
      for(j=0;j<N;j++){
	mets[i] += mettab[(i >> (N-j-1)) & 1][symbols[j]];
      }
    }
    symbols += N;
    /* Run the add-compare-select operations */
    mask = 1;
    for(i=0;i< 1 << (K-1);i+=2){
      int b1,b2;

      b1 = mets[Syms[i]];
      nmetric[i] = m0 = cmetric[i/2] + b1;
      b2 = mets[Syms[i+1]];
      b1 -= b2;
      m1 = cmetric[(i/2) + (1<<(K-2))] + b2;

      if(m1 > m0){
	nmetric[i] = m1;
	*pp |= mask;
      }

      m0 -= b1;
      nmetric[i+1] = m0;
      m1 += b1;

      if(m1 > m0){
	nmetric[i+1] = m1;
	*pp |= mask << 1;
      }

      mask <<= 2;
      if(mask == 0){
	mask = 1;
	pp++;
	ipp++;
      }
    }
    if(mask != 1){
      pp++;
      ipp++;
    }
    if(++bitcnt == nbits){
      *metric = nmetric[endstate];
      break;
    }
    memcpy(cmetric,nmetric,sizeof(cmetric));
  }

  /* Chain back from terminal state to produce decoded data */
  if(data == NULL)
    return 0;/* Discard output */
  memset(data,0,(nbits+7)/8); /* round up in case nbits % 8 != 0 */

  for(i=nbits-1;i >= 0;i--){
    //    int a0,a1;
    pp -= D;
    ipp -= D;
    m0=endstate >> LOGLONGBITS;
    m1=1L << (endstate & (LONGBITS-1));
    if(pp[m0] & m1) {
      //      a0=nmetric[endstate];
      endstate |= (1 << (K-1));
      //      a1=nmetric[endstate];
      data[i>>3] |= 0x80 >> (i&7);
      //      printf("B  %d  %d  %d  %d\n",*metric,i,a0,a1);
    }
    endstate >>= 1;
  }
  return 0;
}


code:
subroutine wav11(d2,npts,dd)

! Convert i*2 data sampled at 12000 Hz to r*4 sampled at 11025 Hz.

  parameter (NZ11=60*11025,NZ12=60*12000)
  parameter (NFFT1=64*12000,NFFT2=64*11025)
  integer*2 d2(NZ12)
  real*4 dd(NZ11)
  real x(NFFT2)
  complex cx(0:NFFT1/2)
  equivalence (x,cx)
  save x,cx

  jz=min(NZ12,npts)
  x(1:jz)=d2(1:jz)
  x(jz+1:)=0.0
  call four2a(x,nfft1,1,-1,0)                    !Forward FFT, r2c
  df=12000.0/NFFT1
  ia=5000.0/df
  cx(ia:)=0.0
  call four2a(cx,nfft2,1,1,-1)                   !Inverse FFT, c2r
  npts=jz*11025.0/12000.0
  fac=1.e-6
  dd(1:npts)=fac*x(1:npts)

  return
end subroutine wav11


code:
subroutine wav12(d2,d1,npts,nbitsam2)

! Convert i*2 or i*1 data at 11025 Hz (from WSJT *.wav files)
! to i*2 data at 12000 Hz.

! Input:  i*2 d2(npts) or i*1 d1(npts)
!         i*2 nbitsam2 = 8 or 16 (bits per sample)

! Output: npts = (12000*npts)/11025
!         i*2 d2(npts)

  parameter (NZ11=60*11025,NZ12=60*12000)
  parameter (NFFT1=64*11025,NFFT2=64*12000)
  integer*1 d1(NZ11)
  integer*1 d1a(NZ11)
  integer*1 i1
  integer*2 i2
  integer*2 d2(NZ12)
  real x(NFFT2)
  complex cx(0:NFFT2/2)
  integer*2 nbitsam2
  equivalence (x,cx),(i1,i2)

  jz=min(NZ11,npts)
  if(nbitsam2.eq.8) then
     jz=min(NZ11,npts)
     d1a(1:jz)=d1(1:jz)            !d1 and d2 may be same array in calling prog
     do i=1,jz                     !Move data from d1a into d2
        i2=0
        i1=d1a(i)
        d2(i)=10*(i2-128)
     enddo
  endif

  x(1:jz)=d2(1:jz)
  x(jz+1:)=0.0
  call four2a(x,nfft1,1,-1,0)                    !Forwarxd FFT, r2c
  cx(nfft1/2:)=0.0
  call four2a(cx,nfft2,1,1,-1)                   !Inverse FFT, c2r

  npts=jz*12000.0/11025.0
  fac=1.e-6
!  if(nbitsam2.eq.16) fac=3.e-6
  x=fac*x
  d2(1:npts)=nint(x(1:npts))
  if(npts.lt.NZ12) d2(npts+1:NZ12)=0

  return
end subroutine wav12


code:
module wavhdr
  type hdr
     character*4 ariff
     integer*4 lenfile
     character*4 awave
     character*4 afmt
     integer*4 lenfmt
     integer*2 nfmt2
     integer*2 nchan2
     integer*4 nsamrate
     integer*4 nbytesec
     integer*2 nbytesam2
     integer*2 nbitsam2
     character*4 adata
     integer*4 ndata
  end type hdr

  contains

    function default_header(nsamrate,npts)
      type(hdr) default_header,h
      h%ariff='RIFF'
      h%awave='WAVE'
      h%afmt='fmt '
      h%lenfmt=16
      h%nfmt2=1
      h%nchan2=1
      h%nsamrate=nsamrate
      h%nbitsam2=16
      h%nbytesam2=h%nbitsam2 * h%nchan2 / 8
      h%adata='data'
      h%nbytesec=h%nsamrate * h%nbitsam2 * h%nchan2 / 8
      h%ndata=2*npts
      h%lenfile=h%ndata + 44 - 8
      default_header=h
    end function default_header

    subroutine set_wsjtx_wav_params(fMHz,mode,nsubmode,ntrperiod,id2)

      parameter (NBANDS=23,NMODES=11)
      character*8 mode,modes(NMODES)
      integer*2 id2(4)
      integer iperiod(7)
      real fband(NBANDS)
      data fband/0.137,0.474,1.8,3.5,5.1,7.0,10.14,14.0,18.1,21.0,24.9,  &
           28.0,50.0,144.0,222.0,432.0,902.0,1296.0,2304.0,3400.0,       &
           5760.0,10368.0,24048.0/
      data modes/'Echo','FSK441','ISCAT','JT4','JT65','JT6M','JT9',      &
           'JT9+JT65','JTMS','JTMSK','WSPR'/
      data iperiod/5,10,15,30,60,120,900/

      dmin=1.e30
      iband=0
      do i=1,NBANDS
         if(abs(fMHz-fband(i)).lt.dmin) then
            dmin=abs(fMHz-fband(i))
            iband=i
         endif
      enddo

      imode=0
      do i=1,NMODES
         if(mode.eq.modes(i)) imode=i
      enddo

      ip=0
      do i=1,7
         if(ntrperiod.eq.iperiod(i)) ip=i
      enddo

      id2(1)=iband
      id2(2)=imode
      id2(3)=nsubmode
      id2(4)=ip

      return
    end subroutine set_wsjtx_wav_params

    subroutine get_wsjtx_wav_params(id2,band,mode,nsubmode,ntrperiod,ok)

      parameter (NBANDS=23,NMODES=11)
      character*8 mode,modes(NMODES)
      character*6 band,bands(NBANDS)
      integer*2 id2(4)
      integer iperiod(7)
      logical ok
      data modes/'Echo','FSK441','ISCAT','JT4','JT65','JT6M','JT9',    &
           'JT9+JT65','JTMS','JTMSK','WSPR'/
      data iperiod/5,10,15,30,60,120,900/
      data bands/'2190m','630m','160m','80m','60m','40m','30m','20m',  &
           '17m','15m','12m','10m','6m','2m','1.25m','70cm','33cm',    &
           '23cm','13cm','9cm','6cm','3cm','1.25cm'/

      ok=.true.
      if(id2(1).lt.1 .or. id2(1).gt.NBANDS) ok=.false.
      if(id2(2).lt.1 .or. id2(2).gt.NMODES) ok=.false.
      if(id2(3).lt.1 .or. id2(3).gt.8) ok=.false.
      if(id2(4).lt.1 .or. id2(4).gt.7) ok=.false.

      if(ok) then
         band=bands(id2(1))
         mode=modes(id2(2))
         nsubmode=id2(3)
         ntrperiod=iperiod(id2(4))
      endif

      return
    end subroutine get_wsjtx_wav_params

end module wavhdr


code:
void export_wisdom_(char fname[], int len)
{
  int fftwf_export_wisdom_to_filename(const char *);
  fname[len-1]=0;
  fftwf_export_wisdom_to_filename(fname);
}

void import_wisdom_(char fname[], int *success, int len)
{
  int fftwf_import_wisdom_from_filename(const char *);
  fname[len-1]=0;
  *success = fftwf_import_wisdom_from_filename(fname);
}


code: C:\JTSDK-QT\fftw3f\fftwf-wisdom -o fftwf_wisdom.dat -t 1 -m rif672000 cif77175 cib77175 rif16384 rif884736 cib2048 rif8192 rif512 rib512 cib512


code:
!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name:    wqdecode.f90
! Description:
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------
subroutine wqdecode(data0,message,ntype)

  parameter (N15=32768)
  integer*1 data0(11)
  character*22 message
  character*12 callsign
  character*3 cdbm
  character grid4*4,grid6*6
  logical first
  character*12 dcall(0:N15-1)
  data first/.true./
  save first,dcall

! May want to have a timeout (say, one hour?) on calls fetched
! from the hash table.

  if(first) then
     dcall='            '
     first=.false.
  endif

  message='                      '
  call unpack50(data0,n1,n2)
!  print*,data0,n1,n2
  call unpackcall(n1,callsign)
  i1=index(callsign,' ')
  call unpackgrid(n2/128,grid4)
  ntype=iand(n2,127) -64

! Standard WSPR message (types 0 3 7 10 13 17 ... 60)
  if(ntype.ge.0 .and. ntype.le.62) then
     nu=mod(ntype,10)
     if(nu.eq.0 .or. nu.eq.3 .or. nu.eq.7) then
        write(cdbm,'(i3)'),ntype
        if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
        if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
        message=callsign(1:i1)//grid4//' '//cdbm
        call hash(callsign,i1-1,ih)
        dcall(ih)=callsign(:i1)
     else
        nadd=nu
        if(nu.gt.3) nadd=nu-3
        if(nu.gt.7) nadd=nu-7
        ng=n2/128 + 32768*(nadd-1)
        call unpackpfx(ng,callsign)
        ndbm=ntype-nadd
        write(cdbm,'(i3)'),ndbm
        if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
        if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
        i2=index(callsign,' ')
        message=callsign(:i2)//cdbm
        call hash(callsign,i2-1,ih)
        dcall(ih)=callsign(:i2)
     endif
  else if(ntype.lt.0) then
     ndbm=-(ntype+1)
     grid6=callsign(6:6)//callsign(1:5)
     ih=(n2-ntype-64)/128
     callsign=dcall(ih)
     write(cdbm,'(i3)'),ndbm
     if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
     if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
     i2=index(callsign,' ')
     if(dcall(ih)(1:1).ne.' ') then
        message='<'//callsign(:i2-1)//'> '//grid6//' '//cdbm
     else
        message='<...> '//grid6//' '//cdbm
     endif
  endif

  return
end subroutine wqdecode

!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name:    unpack50.f90
! Description:
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------
subroutine unpack50(dat,n1,n2)

  integer*1 dat(11)

  i=dat(1)
  i4=iand(i,255)
  n1=ishft(i4,20)
  i=dat(2)
  i4=iand(i,255)
  n1=n1 + ishft(i4,12)
  i=dat(3)
  i4=iand(i,255)
  n1=n1 + ishft(i4,4)
  i=dat(4)
  i4=iand(i,255)
  n1=n1 + iand(ishft(i4,-4),15)
  n2=ishft(iand(i4,15),18)
  i=dat(5)
  i4=iand(i,255)
  n2=n2 + ishft(i4,10)
  i=dat(6)
  i4=iand(i,255)
  n2=n2 + ishft(i4,2)
  i=dat(7)
  i4=iand(i,255)
  n2=n2 + iand(ishft(i4,-6),3)

  return
end subroutine unpack50

!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name:    unpackcall.f90
! Description:
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------
subroutine unpackcall(ncall,word)

  character word*12,c*37

  data c/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ '/

  n=ncall
  word='......'
  if(n.ge.262177560) go to 999            !Plain text message ...
  i=mod(n,27)+11
  word(6:6)=c(i:i)
  n=n/27
  i=mod(n,27)+11
  word(5:5)=c(i:i)
  n=n/27
  i=mod(n,27)+11
  word(4:4)=c(i:i)
  n=n/27
  i=mod(n,10)+1
  word(3:3)=c(i:i)
  n=n/10
  i=mod(n,36)+1
  word(2:2)=c(i:i)
  n=n/36
  i=n+1
  word(1:1)=c(i:i)
  do i=1,4
     if(word(i:i).ne.' ') go to 10
  enddo
  go to 999
10 word=word(i:)

999 if(word(1:3).eq.'3D0') word='3DA0'//word(4:)
  return
end subroutine unpackcall

!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name:    unpackgrid.f90
! Description:
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------
subroutine unpackgrid(ng,grid)

  parameter (NGBASE=180*180)
  character grid*4,grid6*6,digit*10
  data digit/'0123456789'/

  grid='    '
  if(ng.ge.32400) go to 10
  dlat=mod(ng,180)-90
  dlong=(ng/180)*2 - 180 + 2
  call deg2grid(dlong,dlat,grid6)
  grid=grid6(1:4) !XXX explicitly truncate this -db
  go to 100

10 n=ng-NGBASE-1
  if(n.ge.1 .and.n.le.30) then
     grid(1:1)='-'
     grid(2:2)=char(48+n/10)
     grid(3:3)=char(48+mod(n,10))
  else if(n.ge.31 .and.n.le.60) then
     n=n-30
     grid(1:2)='R-'
     grid(3:3)=char(48+n/10)
     grid(4:4)=char(48+mod(n,10))
  else if(n.eq.61) then
     grid='RO'
  else if(n.eq.62) then
     grid='RRR'
  else if(n.eq.63) then
     grid='73'
  endif

100 return
end subroutine unpackgrid

!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name:    unpackpfx.f90
! Description:
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------
subroutine unpackpfx(ng,call1)

  character*12 call1
  character*3 pfx

  if(ng.lt.60000) then
! Add-on prefix of 1 to 3 characters
     n=ng
     do i=3,1,-1
        nc=mod(n,37)
        if(nc.ge.0 .and. nc.le.9) then
           pfx(i:i)=char(nc+48)
        else if(nc.ge.10 .and. nc.le.35) then
           pfx(i:i)=char(nc+55)
        else
           pfx(i:i)=' '
        endif
        n=n/37
     enddo
     call1=pfx//'/'//call1
     if(call1(1:1).eq.' ') call1=call1(2:)
     if(call1(1:1).eq.' ') call1=call1(2:)
  else
! Add-on suffix, one or teo characters
     i1=index(call1,' ')
     nc=ng-60000
     if(nc.ge.0 .and. nc.le.9) then
        call1=call1(:i1-1)//'/'//char(nc+48)
     else if(nc.ge.10 .and. nc.le.35) then
        call1=call1(:i1-1)//'/'//char(nc+55)
     else if(nc.ge.36 .and. nc.le.125) then
        nc1=(nc-26)/10
        nc2=mod(nc-26,10)
        call1=call1(:i1-1)//'/'//char(nc1+48)//char(nc2+48)
     endif
  endif

  return
end subroutine unpackpfx


code:
subroutine wqencode(msg,ntype,data0)

!  Parse and encode a WSPR message.

  use packjt
  parameter (MASK15=32767)
  character*22 msg
  character*12 call1,call2
  character grid4*4
  logical lbad1,lbad2
  integer*1 data0(11)
  integer nu(0:9)
  data nu/0,-1,1,0,-1,2,1,0,-1,1/

! Standard WSPR message (types 0 3 7 10 13 17 ... 60)
  i1=index(msg,' ')
  i2=index(msg,'/')
  i3=index(msg,'<')
  call1=msg(:i1-1)
  if(i1.lt.3 .or. i1.gt.7 .or. i2.gt.0 .or. i3.gt.0) go to 10
  grid4=msg(i1+1:i1+4)
  call packcall(call1,n1,lbad1)
  call packgrid(grid4,ng,lbad2)
  if(lbad1 .or. lbad2) go to 10
  ndbm=0
  read(msg(i1+5:),*) ndbm
  if(ndbm.lt.0) ndbm=0
  if(ndbm.gt.60) ndbm=60
  ndbm=ndbm+nu(mod(ndbm,10))
  n2=128*ng + (ndbm+64)
  call pack50(n1,n2,data0)
  ntype=ndbm
  go to 900

10 if(i2.ge.2 .and. i3.lt.1) then
     call packpfx(call1,n1,ng,nadd)
     ndbm=0
     read(msg(i1+1:),*) ndbm
     if(ndbm.lt.0) ndbm=0
     if(ndbm.gt.60) ndbm=60
     ndbm=ndbm+nu(mod(ndbm,10))
     ntype=ndbm + 1 + nadd
     n2=128*ng + ntype + 64
     call pack50(n1,n2,data0)
  else if(i3.eq.1) then
     i4=index(msg,'>')
     call1=msg(2:i4-1)
     call hash(call1,i4-2,ih)
     i5=index(trim(msg(i1+1:)),' ')
! Convert grid to valid callsign format - first character moved to end
     call2=msg(i1+2:i1+i5-1)//msg(i1+1:i1+1)//'        '
     call packcall(call2,n1,lbad1)
     ndbm=0
     read(msg(i1+i5+1:),*) ndbm
     if(ndbm.lt.0) ndbm=0
     if(ndbm.gt.60) ndbm=60
     ndbm=ndbm+nu(mod(ndbm,10))
     ntype=-(ndbm+1)
     n2=128*ih + ntype + 64
     call pack50(n1,n2,data0)
  endif
  go to 900

900 continue
  return
end subroutine wqencode


code:
#include <math.h>
#include <stdio.h>
#include <float.h>
#include <limits.h>
#include <stdlib.h>
#include "rs.h"

static void *rs;
static int first=1;

void rs_encode_(int *dgen, int *sent)
// Encode JT65 data dgen[12], producing sent[63].
{
  int dat1[12];
  int b[51];
  int i;

  if(first) {
    // Initialize the JT65 codec
    rs=init_rs_int(6,0x43,3,1,51,0);
    first=0;
  }

  // Reverse data order for the Karn codec.
  for(i=0; i<12; i++) {
    dat1[i]=dgen[11-i];
  }
  // Compute the parity symbols
  encode_rs_int(rs,dat1,b);

  // Move parity symbols and data into sent[] array, in reverse order.
  for (i = 0; i < 51; i++) sent[50-i] = b[i];
  for (i = 0; i < 12; i++) sent[i+51] = dat1[11-i];
}

void rs_decode_(int *recd0, int *era0, int *numera0, int *decoded, int *nerr)
// Decode JT65 received data recd0[63], producing decoded[12].
// Erasures are indicated in era0[numera].  The number of corrected
// errors is *nerr.  If the data are uncorrectable, *nerr=-1 is returned.
{
  int numera;
  int i;
  int era_pos[50];
  int recd[63];

  if(first) {
    rs=init_rs_int(6,0x43,3,1,51,0);
    first=0;
  }

  numera=*numera0;
  for(i=0; i<12; i++) recd[i]=recd0[62-i];
  for(i=0; i<51; i++) recd[12+i]=recd0[50-i];
  if(numera)
    for(i=0; i<numera; i++) era_pos[i]=era0[i];
  *nerr=decode_rs_int(rs,recd,era_pos,numera);
  for(i=0; i<12; i++) decoded[i]=recd[11-i];
}

void rs_encode__(int *dgen, int *sent)
{
  rs_encode_(dgen, sent);
}

void rs_decode__(int *recd0, int *era0, int *numera0, int *decoded, int *nerr)
{
  rs_decode_(recd0, era0, numera0, decoded, nerr);
}


code:
subroutine wspr_downsample(id2,k)

! Input:
!  id2       raw 16-bit integer data, 12000 Hz sample rate
!  k         pointer to the most recent new data

! Output (in common/c0com)
!  c0        complex data downsampled to 1500 Hz

  parameter (NMAX=120*12000)         !Total sample intervals per 30 minutes
  parameter (NDMAX=120*1500)         !Sample intervals at 1500 Hz rate
  parameter (NSMAX=1366)             !Max length of saved spectra
  parameter (NFFT1=1024)
  parameter (MAXFFT3=32768)
  real*4 w3(MAXFFT3)
  real*4 x0(NFFT1)
  real*4 x2(NFFT1+105)
  real*4 ssum(NSMAX)
  integer*2 id2(NMAX)
  complex c0
  common/c0com/c0(NDMAX)
  data rms/999.0/,k0/99999999/,nfft3z/0/,nsps/8192/,nbfo/1500/
  save

  nfft3=nsps/4
  jstep=nsps/16
  if(k.gt.NMAX) go to 999
  if(k.lt.nfft3) go to 999                   !Wait for enough samples to start
  if(nfft3.ne.nfft3z) then
     pi=4.0*atan(1.0)
     do i=1,nfft3
        w3(i)=2.0*(sin(i*pi/nfft3))**2             !Window for nfft3
     enddo
     nfft3z=nfft3
  endif

  if(k.lt.k0) then
     ja=0
     ssum=0.
     k1=0
     k8=0
     x2=0.
!     if(ndiskdat.eq.0) then
!        id2(k+1:)=0
!        c0=0.          !This is necessary to prevent "ghosts".  Not sure why.
!     endif
  endif
  k0=k

  nzap=0
  nbslider=0
  sigmas=1.0*(10.0**(0.01*nbslider)) + 0.7
  peaklimit=sigmas*max(10.0,rms)
  px=0.

  nwindow=2
  kstep1=NFFT1
  if(nwindow.ne.0) kstep1=NFFT1/2
  fac=2.0/NFFT1
  nblks=(k-k1)/kstep1
  gain=1.0
  nb=0
  do nblk=1,nblks
     do i=1,NFFT1
        x0(i)=gain*id2(k1+i)
     enddo
!     call timf2(x0,k,NFFT1,nwindow,nb,peaklimit,x1,   &
!          slimit,lstrong,px,nzap)
! Mix at nbfo Hz, lowpass at +/-750 Hz, and downsample to 1500 Hz complex.
     call mixlpf(x0,nbfo,c0(k8+1))
     k1=k1+kstep1
     k8=k8+kstep1/8
  enddo

999 return
end subroutine wspr_downsample


code: gfortran -o wsprcode -fbounds-check wsprcode.f90 nhash.c


code:
/*
 *-------------------------------------------------------------------------------
 *
 * This file is part of the WSPR application, Weak Signal Propagation Reporter
 *
 * File Name:   nhash.c
 * Description: Functions to produce 32-bit hashes for hash table lookup
 *
 * Copyright (C) 2008-2014 Joseph Taylor, K1JT
 * License: GPL-3
 *
 * This program is free software; you can redistribute it and/or modify it under
 * the terms of the GNU General Public License as published by the Free Software
 * Foundation; either version 3 of the License, or (at your option) any later
 * version.
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
 * details.
 *
 * You should have received a copy of the GNU General Public License along with
 * this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
 * Street, Fifth Floor, Boston, MA 02110-1301, USA.
 *
 * Files: lookup3.c
 * Copyright: Copyright (C) 2006 Bob Jenkins <bob_jenkins@burtleburtle.net>
 * License: public-domain
 *  You may use this code any way you wish, private, educational, or commercial.
 *  It's free.
 *
 *-------------------------------------------------------------------------------
*/

/*
These are functions for producing 32-bit hashes for hash table lookup.
hashword(), hashlittle(), hashlittle2(), hashbig(), mix(), and final()
are externally useful functions.  Routines to test the hash are included
if SELF_TEST is defined.  You can use this free for any purpose.  It's in
the public domain.  It has no warranty.

You probably want to use hashlittle().  hashlittle() and hashbig()
hash byte arrays.  hashlittle() is is faster than hashbig() on
little-endian machines.  Intel and AMD are little-endian machines.
On second thought, you probably want hashlittle2(), which is identical to
hashlittle() except it returns two 32-bit hashes for the price of one.
You could implement hashbig2() if you wanted but I haven't bothered here.

If you want to find a hash of, say, exactly 7 integers, do
  a = i1;  b = i2;  c = i3;
  mix(a,b,c);
  a += i4; b += i5; c += i6;
  mix(a,b,c);
  a += i7;
  final(a,b,c);
then use c as the hash value.  If you have a variable length array of
4-byte integers to hash, use hashword().  If you have a byte array (like
a character string), use hashlittle().  If you have several byte arrays, or
a mix of things, see the comments above hashlittle().

Why is this so big?  I read 12 bytes at a time into 3 4-byte integers,
then mix those integers.  This is fast (you can do a lot more thorough
mixing with 12*3 instructions on 3 integers than you can with 3 instructions
on 1 byte), but shoehorning those bytes into integers efficiently is messy.
*/

#define SELF_TEST 1

#include <stdio.h>      /* defines printf for tests */
#include <time.h>       /* defines time_t for timings in the test */
#ifdef Win32
#include "win_stdint.h"	/* defines uint32_t etc */
#else
#include <stdint.h>	/* defines uint32_t etc */
#endif
//#include <sys/param.h>  /* attempt to define endianness */
//#ifdef linux
//# include <endian.h>    /* attempt to define endianness */
//#endif

#define HASH_LITTLE_ENDIAN 1

#define hashsize(n) ((uint32_t)1<<(n))
#define hashmask(n) (hashsize(n)-1)
#define rot(x,k) (((x)<<(k)) | ((x)>>(32-(k))))

/*
-------------------------------------------------------------------------------
mix -- mix 3 32-bit values reversibly.

This is reversible, so any information in (a,b,c) before mix() is
still in (a,b,c) after mix().

If four pairs of (a,b,c) inputs are run through mix(), or through
mix() in reverse, there are at least 32 bits of the output that
are sometimes the same for one pair and different for another pair.
This was tested for:
* pairs that differed by one bit, by two bits, in any combination
  of top bits of (a,b,c), or in any combination of bottom bits of
  (a,b,c).
* "differ" is defined as +, -, ^, or ~^.  For + and -, I transformed
  the output delta to a Gray code (a^(a>>1)) so a string of 1's (as
  is commonly produced by subtraction) look like a single 1-bit
  difference.
* the base values were pseudorandom, all zero but one bit set, or
  all zero plus a counter that starts at zero.

Some k values for my "a-=c; a^=rot(c,k); c+=b;" arrangement that
satisfy this are
    4  6  8 16 19  4
    9 15  3 18 27 15
   14  9  3  7 17  3
Well, "9 15 3 18 27 15" didn't quite get 32 bits diffing
for "differ" defined as + with a one-bit base and a two-bit delta.  I
used http://burtleburtle.net/bob/hash/avalanche.html to choose
the operations, constants, and arrangements of the variables.

This does not achieve avalanche.  There are input bits of (a,b,c)
that fail to affect some output bits of (a,b,c), especially of a.  The
most thoroughly mixed value is c, but it doesn't really even achieve
avalanche in c.

This allows some parallelism.  Read-after-writes are good at doubling
the number of bits affected, so the goal of mixing pulls in the opposite
direction as the goal of parallelism.  I did what I could.  Rotates
seem to cost as much as shifts on every machine I could lay my hands
on, and rotates are much kinder to the top and bottom bits, so I used
rotates.
-------------------------------------------------------------------------------
*/
#define mix(a,b,c) \
{ \
  a -= c;  a ^= rot(c, 4);  c += b; \
  b -= a;  b ^= rot(a, 6);  a += c; \
  c -= b;  c ^= rot(b, 8);  b += a; \
  a -= c;  a ^= rot(c,16);  c += b; \
  b -= a;  b ^= rot(a,19);  a += c; \
  c -= b;  c ^= rot(b, 4);  b += a; \
}

/*
-------------------------------------------------------------------------------
final -- final mixing of 3 32-bit values (a,b,c) into c

Pairs of (a,b,c) values differing in only a few bits will usually
produce values of c that look totally different.  This was tested for
* pairs that differed by one bit, by two bits, in any combination
  of top bits of (a,b,c), or in any combination of bottom bits of
  (a,b,c).
* "differ" is defined as +, -, ^, or ~^.  For + and -, I transformed
  the output delta to a Gray code (a^(a>>1)) so a string of 1's (as
  is commonly produced by subtraction) look like a single 1-bit
  difference.
* the base values were pseudorandom, all zero but one bit set, or
  all zero plus a counter that starts at zero.

These constants passed:
 14 11 25 16 4 14 24
 12 14 25 16 4 14 24
and these came close:
  4  8 15 26 3 22 24
 10  8 15 26 3 22 24
 11  8 15 26 3 22 24
-------------------------------------------------------------------------------
*/
#define final(a,b,c) \
{ \
  c ^= b; c -= rot(b,14); \
  a ^= c; a -= rot(c,11); \
  b ^= a; b -= rot(a,25); \
  c ^= b; c -= rot(b,16); \
  a ^= c; a -= rot(c,4);  \
  b ^= a; b -= rot(a,14); \
  c ^= b; c -= rot(b,24); \
}

/*
-------------------------------------------------------------------------------
hashlittle() -- hash a variable-length key into a 32-bit value
  k       : the key (the unaligned variable-length array of bytes)
  length  : the length of the key, counting by bytes
  initval : can be any 4-byte value
Returns a 32-bit value.  Every bit of the key affects every bit of
the return value.  Two keys differing by one or two bits will have
totally different hash values.

The best hash table sizes are powers of 2.  There is no need to do
mod a prime (mod is sooo slow!).  If you need less than 32 bits,
use a bitmask.  For example, if you need only 10 bits, do
  h = (h & hashmask(10));
In which case, the hash table should have hashsize(10) elements.

If you are hashing n strings (uint8_t **)k, do it like this:
  for (i=0, h=0; i<n; ++i) h = hashlittle( k[i], len[i], h);

By Bob Jenkins, 2006.  bob_jenkins@burtleburtle.net.  You may use this
code any way you wish, private, educational, or commercial.  It's free.

Use for hash table lookup, or anything where one collision in 2^^32 is
acceptable.  Do NOT use for cryptographic purposes.
-------------------------------------------------------------------------------
*/

//uint32_t hashlittle( const void *key, size_t length, uint32_t initval)
#ifdef STDCALL
uint32_t __stdcall NHASH( const void *key, size_t *length0, uint32_t *initval0)
#else
uint32_t nhash_( const void *key, int *length0, uint32_t *initval0)
#endif
{
  uint32_t a,b,c;                                          /* internal state */
  size_t length;
  uint32_t initval;
  union { const void *ptr; size_t i; } u;     /* needed for Mac Powerbook G4 */

  length=*length0;
  initval=*initval0;

  /* Set up the internal state */
  a = b = c = 0xdeadbeef + ((uint32_t)length) + initval;

  u.ptr = key;
  if (HASH_LITTLE_ENDIAN && ((u.i & 0x3) == 0)) {
    const uint32_t *k = (const uint32_t *)key;         /* read 32-bit chunks */
    const uint8_t  *k8;

    k8=0;                                     //Silence compiler warning
    /*------ all but last block: aligned reads and affect 32 bits of (a,b,c) */
    while (length > 12)
    {
      a += k[0];
      b += k[1];
      c += k[2];
      mix(a,b,c);
      length -= 12;
      k += 3;
    }

    /*----------------------------- handle the last (probably partial) block */
    /*
     * "k[2]&0xffffff" actually reads beyond the end of the string, but
     * then masks off the part it's not allowed to read.  Because the
     * string is aligned, the masked-off tail is in the same word as the
     * rest of the string.  Every machine with memory protection I've seen
     * does it on word boundaries, so is OK with this.  But VALGRIND will
     * still catch it and complain.  The masking trick does make the hash
     * noticably faster for short strings (like English words).
     */
#ifndef VALGRIND

    switch(length)
    {
    case 12: c+=k[2]; b+=k[1]; a+=k[0]; break;
    case 11: c+=k[2]&0xffffff; b+=k[1]; a+=k[0]; break;
    case 10: c+=k[2]&0xffff; b+=k[1]; a+=k[0]; break;
    case 9 : c+=k[2]&0xff; b+=k[1]; a+=k[0]; break;
    case 8 : b+=k[1]; a+=k[0]; break;
    case 7 : b+=k[1]&0xffffff; a+=k[0]; break;
    case 6 : b+=k[1]&0xffff; a+=k[0]; break;
    case 5 : b+=k[1]&0xff; a+=k[0]; break;
    case 4 : a+=k[0]; break;
    case 3 : a+=k[0]&0xffffff; break;
    case 2 : a+=k[0]&0xffff; break;
    case 1 : a+=k[0]&0xff; break;
    case 0 : return c;              /* zero length strings require no mixing */
    }

#else /* make valgrind happy */

    k8 = (const uint8_t *)k;
    switch(length)
    {
    case 12: c+=k[2]; b+=k[1]; a+=k[0]; break;
    case 11: c+=((uint32_t)k8[10])<<16;  /* fall through */
    case 10: c+=((uint32_t)k8[9])<<8;    /* fall through */
    case 9 : c+=k8[8];                   /* fall through */
    case 8 : b+=k[1]; a+=k[0]; break;
    case 7 : b+=((uint32_t)k8[6])<<16;   /* fall through */
    case 6 : b+=((uint32_t)k8[5])<<8;    /* fall through */
    case 5 : b+=k8[4];                   /* fall through */
    case 4 : a+=k[0]; break;
    case 3 : a+=((uint32_t)k8[2])<<16;   /* fall through */
    case 2 : a+=((uint32_t)k8[1])<<8;    /* fall through */
    case 1 : a+=k8[0]; break;
    case 0 : return c;
    }

#endif /* !valgrind */

  } else if (HASH_LITTLE_ENDIAN && ((u.i & 0x1) == 0)) {
    const uint16_t *k = (const uint16_t *)key;         /* read 16-bit chunks */
    const uint8_t  *k8;

    /*--------------- all but last block: aligned reads and different mixing */
    while (length > 12)
    {
      a += k[0] + (((uint32_t)k[1])<<16);
      b += k[2] + (((uint32_t)k[3])<<16);
      c += k[4] + (((uint32_t)k[5])<<16);
      mix(a,b,c);
      length -= 12;
      k += 6;
    }

    /*----------------------------- handle the last (probably partial) block */
    k8 = (const uint8_t *)k;
    switch(length)
    {
    case 12: c+=k[4]+(((uint32_t)k[5])<<16);
             b+=k[2]+(((uint32_t)k[3])<<16);
             a+=k[0]+(((uint32_t)k[1])<<16);
             break;
    case 11: c+=((uint32_t)k8[10])<<16;     /* fall through */
    case 10: c+=k[4];
             b+=k[2]+(((uint32_t)k[3])<<16);
             a+=k[0]+(((uint32_t)k[1])<<16);
             break;
    case 9 : c+=k8[8];                      /* fall through */
    case 8 : b+=k[2]+(((uint32_t)k[3])<<16);
             a+=k[0]+(((uint32_t)k[1])<<16);
             break;
    case 7 : b+=((uint32_t)k8[6])<<16;      /* fall through */
    case 6 : b+=k[2];
             a+=k[0]+(((uint32_t)k[1])<<16);
             break;
    case 5 : b+=k8[4];                      /* fall through */
    case 4 : a+=k[0]+(((uint32_t)k[1])<<16);
             break;
    case 3 : a+=((uint32_t)k8[2])<<16;      /* fall through */
    case 2 : a+=k[0];
             break;
    case 1 : a+=k8[0];
             break;
    case 0 : return c;                     /* zero length requires no mixing */
    }

  } else {                        /* need to read the key one byte at a time */
    const uint8_t *k = (const uint8_t *)key;

    /*--------------- all but the last block: affect some 32 bits of (a,b,c) */
    while (length > 12)
    {
      a += k[0];
      a += ((uint32_t)k[1])<<8;
      a += ((uint32_t)k[2])<<16;
      a += ((uint32_t)k[3])<<24;
      b += k[4];
      b += ((uint32_t)k[5])<<8;
      b += ((uint32_t)k[6])<<16;
      b += ((uint32_t)k[7])<<24;
      c += k[8];
      c += ((uint32_t)k[9])<<8;
      c += ((uint32_t)k[10])<<16;
      c += ((uint32_t)k[11])<<24;
      mix(a,b,c);
      length -= 12;
      k += 12;
    }

    /*-------------------------------- last block: affect all 32 bits of (c) */
    switch(length)                   /* all the case statements fall through */
    {
    case 12: c+=((uint32_t)k[11])<<24;
    case 11: c+=((uint32_t)k[10])<<16;
    case 10: c+=((uint32_t)k[9])<<8;
    case 9 : c+=k[8];
    case 8 : b+=((uint32_t)k[7])<<24;
    case 7 : b+=((uint32_t)k[6])<<16;
    case 6 : b+=((uint32_t)k[5])<<8;
    case 5 : b+=k[4];
    case 4 : a+=((uint32_t)k[3])<<24;
    case 3 : a+=((uint32_t)k[2])<<16;
    case 2 : a+=((uint32_t)k[1])<<8;
    case 1 : a+=k[0];
             break;
    case 0 : return c;
    }
  }

  final(a,b,c);
  return c;
}

//uint32_t __stdcall NHASH(const void *key, size_t length, uint32_t initval)


code:
!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name:    wspr_old_subs.f90
! Description:  Utility subroutines from WSPR 2.0
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------

subroutine deg2grid(dlong0,dlat,grid)

  real dlong                        !West longitude (deg)
  real dlat                         !Latitude (deg)
  character grid*6

  dlong=dlong0
  if(dlong.lt.-180.0) dlong=dlong+360.0
  if(dlong.gt.180.0) dlong=dlong-360.0

! Convert to units of 5 min of longitude, working east from 180 deg.
  nlong=60.0*(180.0-dlong)/5.0
  n1=nlong/240                      !20-degree field
  n2=(nlong-240*n1)/24              !2 degree square
  n3=nlong-240*n1-24*n2             !5 minute subsquare
  grid(1:1)=char(ichar('A')+n1)
  grid(3:3)=char(ichar('0')+n2)
  grid(5:5)=char(ichar('a')+n3)

! Convert to units of 2.5 min of latitude, working north from -90 deg.
  nlat=60.0*(dlat+90)/2.5
  n1=nlat/240                       !10-degree field
  n2=(nlat-240*n1)/24               !1 degree square
  n3=nlat-240*n1-24*n2              !2.5 minuts subsquare
  grid(2:2)=char(ichar('A')+n1)
  grid(4:4)=char(ichar('0')+n2)
  grid(6:6)=char(ichar('a')+n3)

  return
end subroutine deg2grid

subroutine encode232(dat,nbytes,symbol,maxsym)

! Convolutional encoder for a K=32, r=1/2 code.

  integer*1 dat(nbytes)             !User data, packed 8 bits per byte
  integer*1 symbol(maxsym)          !Channel symbols, one bit per byte
  integer*1 i1

! Layland-Lushbaugh polynomials for a K=32, r=1/2 convolutional code,
! and 8-bit parity lookup table.

  data npoly1/-221228207/,npoly2/-463389625/
  integer*1 partab(0:255)
  data partab/                 &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0/

  nstate=0
  k=0
  do j=1,nbytes
     do i=7,0,-1
        i1=dat(j)
        i4=i1
        if (i4.lt.0) i4=i4+256
        nstate=ior(ishft(nstate,1),iand(ishft(i4,-i),1))
        n=iand(nstate,npoly1)
        n=ieor(n,ishft(n,-16))
        k=k+1
        symbol(k)=partab(iand(ieor(n,ishft(n,-8)),255))
        n=iand(nstate,npoly2)
        n=ieor(n,ishft(n,-16))
        k=k+1
        symbol(k)=partab(iand(ieor(n,ishft(n,-8)),255))
     enddo
  enddo

  return
end subroutine encode232

subroutine fano232(symbol,nbits,mettab,ndelta,maxcycles,dat,ncycles,metric,ierr)

! Sequential decoder for K=32, r=1/2 convolutional code using
! the Fano algorithm.  Translated from C routine for same purpose
! written by Phil Karn, KA9Q.

  parameter (MAXBITS=103)
  parameter (MAXDAT=(MAXBITS+7)/8)
  integer*1 symbol(0:2*MAXBITS-1)
  integer*1 dat(MAXDAT)               !Decoded user data, 8 bits per byte
  integer mettab(0:255,0:1)           !Metric table

! These were the "node" structure in Karn's C code:
  integer nstate(0:MAXBITS-1)      !Encoder state of next node
  integer gamma(0:MAXBITS-1)       !Cumulative metric to this node
  integer metrics(0:3,0:MAXBITS-1) !Metrics indexed by all possible Tx syms
  integer tm(0:1,0:MAXBITS-1)      !Sorted metrics for current hypotheses
  integer ii(0:MAXBITS-1)          !Current branch being tested

  logical noback

! Layland-Lushbaugh polynomials for a K=32, r=1/2 convolutional code,
! and 8-bit parity lookup table.

  data npoly1/-221228207/,npoly2/-463389625/
  integer*1 partab(0:255)
  data partab/                 &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       0, 1, 1, 0, 1, 0, 0, 1, &
       1, 0, 0, 1, 0, 1, 1, 0/

  ntail=nbits-31

! Compute all possible branch metrics for each symbol pair.
! This is the only place we actually look at the raw input symbols
  i4a=0
  i4b=0
  do np=0,nbits-1
     j=2*np
     i4a=symbol(j)
     i4b=symbol(j+1)
     if (i4a.lt.0) i4a=i4a+256
     if (i4b.lt.0) i4b=i4b+256
     metrics(0,np) = mettab(i4a,0) + mettab(i4b,0)
     metrics(1,np) = mettab(i4a,0) + mettab(i4b,1)
     metrics(2,np) = mettab(i4a,1) + mettab(i4b,0)
     metrics(3,np) = mettab(i4a,1) + mettab(i4b,1)
  enddo

  np=0
  nstate(np)=0

! Compute and sort branch metrics from the root node
  n=iand(nstate(np),npoly1)
  n=ieor(n,ishft(n,-16))
  lsym=partab(iand(ieor(n,ishft(n,-8)),255))
  n=iand(nstate(np),npoly2)
  n=ieor(n,ishft(n,-16))
  lsym=lsym+lsym+partab(iand(ieor(n,ishft(n,-8)),255))
  m0=metrics(lsym,np)
  m1=metrics(ieor(3,lsym),np)
  if(m0.gt.m1) then
     tm(0,np)=m0                      !0-branch has better metric
     tm(1,np)=m1
  else
     tm(0,np)=m1                      !1-branch is better
     tm(1,np)=m0
     nstate(np)=nstate(np) + 1        !Set low bit
  endif

! Start with best branch
  ii(np)=0
  gamma(np)=0
  nt=0

! Start the Fano decoder
  do i=1,nbits*maxcycles
! Look forward
     ngamma=gamma(np) + tm(ii(np),np)
     if(ngamma.ge.nt) then

! Node is acceptable.  If first time visiting this node, tighten threshold:
        if(gamma(np).lt.(nt+ndelta)) nt=nt +                     &
             ndelta * ((ngamma-nt)/ndelta)

! Move forward
        gamma(np+1)=ngamma
        nstate(np+1)=ishft(nstate(np),1)
        np=np+1
        if(np.eq.nbits-1) go to 100     !We're done!

        n=iand(nstate(np),npoly1)
        n=ieor(n,ishft(n,-16))
        lsym=partab(iand(ieor(n,ishft(n,-8)),255))
        n=iand(nstate(np),npoly2)
        n=ieor(n,ishft(n,-16))
        lsym=lsym+lsym+partab(iand(ieor(n,ishft(n,-8)),255))

        if(np.ge.ntail) then
           tm(0,np)=metrics(lsym,np)      !We're in the tail, all zeros
        else
           m0=metrics(lsym,np)
           m1=metrics(ieor(3,lsym),np)
           if(m0.gt.m1) then
              tm(0,np)=m0                 !0-branch has better metric
              tm(1,np)=m1
           else
              tm(0,np)=m1                 !1-branch is better
              tm(1,np)=m0
              nstate(np)=nstate(np) + 1   !Set low bit
           endif
        endif

        ii(np)=0                          !Start with best branch
        go to 99
     endif

! Threshold violated, can't go forward
10   noback=.false.
     if(np.eq.0) noback=.true.
     if(np.gt.0) then
        if(gamma(np-1).lt.nt) noback=.true.
     endif

     if(noback) then
! Can't back up, either.  Relax threshold and look forward again
! to a better branch.
        nt=nt-ndelta
        if(ii(np).ne.0) then
           ii(np)=0
           nstate(np)=ieor(nstate(np),1)
        endif
        go to 99
     endif

! Back up
     np=np-1
     if(np.lt.ntail .and. ii(np).ne.1) then
! Search the next best branch
        ii(np)=ii(np)+1
        nstate(np)=ieor(nstate(np),1)
        go to 99
     endif
     go to 10
99   continue
  enddo
  i=nbits*maxcycles

100 metric=gamma(np)                       !Final path metric

! Copy decoded data to user's buffer
  nbytes=(nbits+7)/8
  np=7
  do j=1,nbytes-1
     i4a=nstate(np)
     dat(j)=i4a
     np=np+8
  enddo
  dat(nbytes)=0

  ncycles=i+1
  ierr=0
  if(i.ge.maxcycles*nbits) ierr=-1

  return
end subroutine fano232

subroutine grid2deg(grid0,dlong,dlat)

! Converts Maidenhead grid locator to degrees of West longitude
! and North latitude.

  character*6 grid0,grid
  character*1 g1,g2,g3,g4,g5,g6

  grid=grid0
  i=ichar(grid(5:5))
  if(grid(5:5).eq.' ' .or. i.le.64 .or. i.ge.128) grid(5:6)='mm'

  if(grid(1:1).ge.'a' .and. grid(1:1).le.'z') grid(1:1)=            &
       char(ichar(grid(1:1))+ichar('A')-ichar('a'))
  if(grid(2:2).ge.'a' .and. grid(2:2).le.'z') grid(2:2)=            &
       char(ichar(grid(2:2))+ichar('A')-ichar('a'))
  if(grid(5:5).ge.'A' .and. grid(5:5).le.'Z') grid(5:5)=            &
       char(ichar(grid(5:5))-ichar('A')+ichar('a'))
  if(grid(6:6).ge.'A' .and. grid(6:6).le.'Z') grid(6:6)=            &
       char(ichar(grid(6:6))-ichar('A')+ichar('a'))

  g1=grid(1:1)
  g2=grid(2:2)
  g3=grid(3:3)
  g4=grid(4:4)
  g5=grid(5:5)
  g6=grid(6:6)

  nlong = 180 - 20*(ichar(g1)-ichar('A'))
  n20d = 2*(ichar(g3)-ichar('0'))
  xminlong = 5*(ichar(g5)-ichar('a')+0.5)
  dlong = nlong - n20d - xminlong/60.0
  nlat = -90+10*(ichar(g2)-ichar('A')) + ichar(g4)-ichar('0')
  xminlat = 2.5*(ichar(g6)-ichar('a')+0.5)
  dlat = nlat + xminlat/60.0

  return
end subroutine grid2deg

subroutine hash(string,len,ihash)

  parameter (MASK15=32767)
  character*(*) string
  integer*1 ic(12)

     do i=1,len
        ic(i)=ichar(string(i:i))
     enddo
     i=nhash(ic,len,146)
     ihash=iand(i,MASK15)

!     print*,'C',ihash,len,string
  return
end subroutine hash

subroutine inter_mept(id,ndir)

! Interleave (ndir=1) or de-interleave (ndir=-1) the array id.

  integer*1 id(0:161),itmp(0:161)
  integer j0(0:161)
  logical first
  data first/.true./
  save

  if(first) then
! Compute the interleave table using bit reversal.
     k=-1
     do i=0,255
        n=0
        ii=i
        do j=0,7
           n=n+n
           if(iand(ii,1).ne.0) n=n+1
           ii=ii/2
        enddo
        if(n.le.161) then
           k=k+1
           j0(k)=n
        endif
     enddo
     first=.false.
  endif

  if(ndir.eq.1) then
     do i=0,161
        itmp(j0(i))=id(i)
     enddo
  else
     do i=0,161
        itmp(i)=id(j0(i))
     enddo
  endif

  do i=0,161
     id(i)=itmp(i)
  enddo

  return
end subroutine inter_mept

function nchar(c)

! Convert ASCII number, letter, or space to 0-36 for callsign packing.

  character c*1
  data n/0/                            !Silence compiler warning

  if(c.ge.'0' .and. c.le.'9') then
     n=ichar(c)-ichar('0')
  else if(c.ge.'A' .and. c.le.'Z') then
     n=ichar(c)-ichar('A') + 10
  else if(c.ge.'a' .and. c.le.'z') then
     n=ichar(c)-ichar('a') + 10
  else if(c.ge.' ') then
     n=36
  else
     Print*,'Invalid character in callsign ',c,' ',ichar(c)
     stop
  endif
  nchar=n

  return
end function nchar

subroutine pack50(n1,n2,dat)

  integer*1 dat(11),i1

  i1=iand(ishft(n1,-20),255)                !8 bits
  dat(1)=i1
  i1=iand(ishft(n1,-12),255)                 !8 bits
  dat(2)=i1
  i1=iand(ishft(n1, -4),255)                 !8 bits
  dat(3)=i1
  i1=16*iand(n1,15)+iand(ishft(n2,-18),15)   !4+4 bits
  dat(4)=i1
  i1=iand(ishft(n2,-10),255)                 !8 bits
  dat(5)=i1
  i1=iand(ishft(n2, -2),255)                 !8 bits
  dat(6)=i1
  i1=64*iand(n2,3)                           !2 bits
  dat(7)=i1
  dat(8)=0
  dat(9)=0
  dat(10)=0
  dat(11)=0

  return
end subroutine pack50

subroutine packcall(callsign,ncall,text)

! Pack a valid callsign into a 28-bit integer.

  parameter (NBASE=37*36*10*27*27*27)
  character callsign*6,c*1,tmp*6,digit*10
  logical text
  data digit/'0123456789'/

  text=.false.

! Work-around for Swaziland prefix:
  if(callsign(1:4).eq.'3DA0') callsign='3D0'//callsign(5:6)

  if(callsign(1:3).eq.'CQ ') then
     ncall=NBASE + 1
     if(callsign(4:4).ge.'0' .and. callsign(4:4).le.'9' .and.       &
          callsign(5:5).ge.'0' .and. callsign(5:5).le.'9' .and.     &
          callsign(6:6).ge.'0' .and. callsign(6:6).le.'9') then
        nfreq=100*(ichar(callsign(4:4))-48) +                       &
             10*(ichar(callsign(5:5))-48) +                         &
             ichar(callsign(6:6))-48
        ncall=NBASE + 3 + nfreq
     endif
     return
  else if(callsign(1:4).eq.'QRZ ') then
     ncall=NBASE + 2
     return
  endif

  tmp='      '
  if(callsign(3:3).ge.'0' .and. callsign(3:3).le.'9') then
     tmp=callsign
  else if(callsign(2:2).ge.'0' .and. callsign(2:2).le.'9') then
     if(callsign(6:6).ne.' ') then
        text=.true.
        return
     endif
     tmp=' '//callsign
  else
     text=.true.
     return
  endif

  do i=1,6
     c=tmp(i:i)
     if(c.ge.'a' .and. c.le.'z')                             &
          tmp(i:i)=char(ichar(c)-ichar('a')+ichar('A'))
  enddo

  n1=0
  if((tmp(1:1).ge.'A'.and.tmp(1:1).le.'Z').or.tmp(1:1).eq.' ') n1=1
  if(tmp(1:1).ge.'0' .and. tmp(1:1).le.'9') n1=1
  n2=0
  if(tmp(2:2).ge.'A' .and. tmp(2:2).le.'Z') n2=1
  if(tmp(2:2).ge.'0' .and. tmp(2:2).le.'9') n2=1
  n3=0
  if(tmp(3:3).ge.'0' .and. tmp(3:3).le.'9') n3=1
  n4=0
  if((tmp(4:4).ge.'A'.and.tmp(4:4).le.'Z').or.tmp(4:4).eq.' ') n4=1
  n5=0
  if((tmp(5:5).ge.'A'.and.tmp(5:5).le.'Z').or.tmp(5:5).eq.' ') n5=1
  n6=0
  if((tmp(6:6).ge.'A'.and.tmp(6:6).le.'Z').or.tmp(6:6).eq.' ') n6=1

  if(n1+n2+n3+n4+n5+n6 .ne. 6) then
     text=.true.
     return
  endif

  ncall=nchar(tmp(1:1))
  ncall=36*ncall+nchar(tmp(2:2))
  ncall=10*ncall+nchar(tmp(3:3))
  ncall=27*ncall+nchar(tmp(4:4))-10
  ncall=27*ncall+nchar(tmp(5:5))-10
  ncall=27*ncall+nchar(tmp(6:6))-10

  return
end subroutine packcall

subroutine packgrid(grid,ng,text)

  parameter (NGBASE=180*180)
  character*4 grid
  logical text

  text=.false.
  if(grid.eq.'    ') go to 90                 !Blank grid is OK

! Test for numerical signal report, etc.
  if(grid(1:1).eq.'-') then
     n=10*(ichar(grid(2:2))-48) + ichar(grid(3:3)) - 48
     ng=NGBASE+1+n
     go to 100
  else if(grid(1:2).eq.'R-') then
     n=10*(ichar(grid(3:3))-48) + ichar(grid(4:4)) - 48
     if(n.eq.0) go to 90
     ng=NGBASE+31+n
     go to 100
  else if(grid(1:2).eq.'RO') then
     ng=NGBASE+62
     go to 100
  else if(grid(1:3).eq.'RRR') then
     ng=NGBASE+63
     go to 100
  else if(grid(1:2).eq.'73') then
     ng=NGBASE+64
     go to 100
  endif

  if(grid(1:1).lt.'A' .or. grid(1:1).gt.'R') text=.true.
  if(grid(2:2).lt.'A' .or. grid(2:2).gt.'R') text=.true.
  if(grid(3:3).lt.'0' .or. grid(3:3).gt.'9') text=.true.
  if(grid(4:4).lt.'0' .or. grid(4:4).gt.'9') text=.true.
  if(text) go to 100

  call grid2deg(grid//'mm',dlong,dlat)
  long=dlong
  lat=dlat+ 90.0
  ng=((long+180)/2)*180 + lat
  go to 100

90 ng=NGBASE + 1

100 return
end subroutine packgrid

subroutine packpfx(call1,n1,ng,nadd)

  character*12 call1,call0
  character*3 pfx
  logical text

  i1=index(call1,'/')
  if(call1(i1+2:i1+2).eq.' ') then
! Single-character add-on suffix (maybe also fourth suffix letter?)
     call0=call1(:i1-1)
     call packcall(call0,n1,text)
     nadd=1
     nc=ichar(call1(i1+1:i1+1))
     if(nc.ge.48 .and. nc.le.57) then
        n=nc-48
     else if(nc.ge.65 .and. nc.le.90) then
        n=nc-65+10
     else
        n=38
     endif
     nadd=1
     ng=60000-32768+n
  else if(call1(i1+3:i1+3).eq.' ') then
! Two-character numerical suffix, /10 to /99
     call0=call1(:i1-1)
     call packcall(call0,n1,text)
     nadd=1
     n=10*(ichar(call1(i1+1:i1+1))-48) + ichar(call1(i1+2:i1+2)) - 48
     nadd=1
     ng=60000 + 26 + n
  else
! Prefix of 1 to 3 characters
     pfx=call1(:i1-1)
     if(pfx(3:3).eq.' ') pfx=' '//pfx
     if(pfx(3:3).eq.' ') pfx=' '//pfx
     call0=call1(i1+1:)
     call packcall(call0,n1,text)

     ng=0
     do i=1,3
        nc=ichar(pfx(i:i))
        if(nc.ge.48 .and. nc.le.57) then
           n=nc-48
        else if(nc.ge.65 .and. nc.le.90) then
           n=nc-65+10
        else
           n=36
        endif
        ng=37*ng + n
     enddo
     nadd=0
     if(ng.ge.32768) then
        ng=ng-32768
        nadd=1
     endif
  endif

  return
end subroutine packpfx

subroutine unpack50(dat,n1,n2)

  integer*1 dat(11)

  i=dat(1)
  i4=iand(i,255)
  n1=ishft(i4,20)
  i=dat(2)
  i4=iand(i,255)
  n1=n1 + ishft(i4,12)
  i=dat(3)
  i4=iand(i,255)
  n1=n1 + ishft(i4,4)
  i=dat(4)
  i4=iand(i,255)
  n1=n1 + iand(ishft(i4,-4),15)
  n2=ishft(iand(i4,15),18)
  i=dat(5)
  i4=iand(i,255)
  n2=n2 + ishft(i4,10)
  i=dat(6)
  i4=iand(i,255)
  n2=n2 + ishft(i4,2)
  i=dat(7)
  i4=iand(i,255)
  n2=n2 + iand(ishft(i4,-6),3)

  return
end subroutine unpack50

subroutine unpackcall(ncall,word)

  character word*12,c*37

  data c/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ '/

  n=ncall
  word='......'
  if(n.ge.262177560) go to 999            !Plain text message ...
  i=mod(n,27)+11
  word(6:6)=c(i:i)
  n=n/27
  i=mod(n,27)+11
  word(5:5)=c(i:i)
  n=n/27
  i=mod(n,27)+11
  word(4:4)=c(i:i)
  n=n/27
  i=mod(n,10)+1
  word(3:3)=c(i:i)
  n=n/10
  i=mod(n,36)+1
  word(2:2)=c(i:i)
  n=n/36
  i=n+1
  word(1:1)=c(i:i)
  do i=1,4
     if(word(i:i).ne.' ') go to 10
  enddo
  go to 999
10 word=word(i:)

999 if(word(1:3).eq.'3D0') word='3DA0'//word(4:)
  return
end subroutine unpackcall

subroutine unpackgrid(ng,grid)

  parameter (NGBASE=180*180)
  character grid*4,grid6*6,digit*10
  data digit/'0123456789'/

  grid='    '
  if(ng.ge.32400) go to 10
  dlat=mod(ng,180)-90
  dlong=(ng/180)*2 - 180 + 2
  call deg2grid(dlong,dlat,grid6)
  grid=grid6(1:4) !XXX explicitly truncate this -db
  go to 100

10 n=ng-NGBASE-1
  if(n.ge.1 .and.n.le.30) then
     grid(1:1)='-'
     grid(2:2)=char(48+n/10)
     grid(3:3)=char(48+mod(n,10))
  else if(n.ge.31 .and.n.le.60) then
     n=n-30
     grid(1:2)='R-'
     grid(3:3)=char(48+n/10)
     grid(4:4)=char(48+mod(n,10))
  else if(n.eq.61) then
     grid='RO'
  else if(n.eq.62) then
     grid='RRR'
  else if(n.eq.63) then
     grid='73'
  endif

100 return
end subroutine unpackgrid

subroutine unpackpfx(ng,call1)

  character*12 call1
  character*3 pfx

  if(ng.lt.60000) then
! Add-on prefix of 1 to 3 characters
     n=ng
     do i=3,1,-1
        nc=mod(n,37)
        if(nc.ge.0 .and. nc.le.9) then
           pfx(i:i)=char(nc+48)
        else if(nc.ge.10 .and. nc.le.35) then
           pfx(i:i)=char(nc+55)
        else
           pfx(i:i)=' '
        endif
        n=n/37
     enddo
     call1=pfx//'/'//call1
     if(call1(1:1).eq.' ') call1=call1(2:)
     if(call1(1:1).eq.' ') call1=call1(2:)
  else
! Add-on suffix, one or teo characters
     i1=index(call1,' ')
     nc=ng-60000
     if(nc.ge.0 .and. nc.le.9) then
        call1=call1(:i1-1)//'/'//char(nc+48)
     else if(nc.ge.10 .and. nc.le.35) then
        call1=call1(:i1-1)//'/'//char(nc+55)
     else if(nc.ge.36 .and. nc.le.125) then
        nc1=(nc-26)/10
        nc2=mod(nc-26,10)
        call1=call1(:i1-1)//'/'//char(nc1+48)//char(nc2+48)
     endif
  endif

  return
end subroutine unpackpfx

subroutine wqdecode(data0,message,ntype)

  parameter (N15=32768)
  integer*1 data0(11)
  character*22 message
  character*12 callsign
  character*3 cdbm
  character grid4*4,grid6*6
  logical first
  character*12 dcall(0:N15-1)
  data first/.true./
  save first,dcall

! May want to have a timeout (say, one hour?) on calls fetched
! from the hash table.

  if(first) then
     dcall='            '
     first=.false.
  endif

  message='                      '
  call unpack50(data0,n1,n2)
  call unpackcall(n1,callsign)
  i1=index(callsign,' ')
  call unpackgrid(n2/128,grid4)
  ntype=iand(n2,127) -64

! Standard WSPR message (types 0 3 7 10 13 17 ... 60)
  if(ntype.ge.0 .and. ntype.le.62) then
     nu=mod(ntype,10)
     if(nu.eq.0 .or. nu.eq.3 .or. nu.eq.7) then
        write(cdbm,'(i3)'),ntype
        if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
        if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
        message=callsign(1:i1)//grid4//' '//cdbm
        call hash(callsign,i1-1,ih)
        dcall(ih)=callsign(:i1)
     else
        nadd=nu
        if(nu.gt.3) nadd=nu-3
        if(nu.gt.7) nadd=nu-7
        ng=n2/128 + 32768*(nadd-1)
        call unpackpfx(ng,callsign)
        ndbm=ntype-nadd
        write(cdbm,'(i3)'),ndbm
        if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
        if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
        i2=index(callsign,' ')
        message=callsign(:i2)//cdbm
        call hash(callsign,i2-1,ih)
        dcall(ih)=callsign(:i2)
     endif
  else if(ntype.lt.0) then
     ndbm=-(ntype+1)
     grid6=callsign(6:6)//callsign(1:5)
     ih=(n2-ntype-64)/128
     callsign=dcall(ih)
     write(cdbm,'(i3)'),ndbm
     if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
     if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
     i2=index(callsign,' ')
     if(dcall(ih)(1:1).ne.' ') then
        message='<'//callsign(:i2-1)//'> '//grid6//' '//cdbm
     else
        message='<...> '//grid6//' '//cdbm
     endif
  endif

  return
end subroutine wqdecode

subroutine wqencode(msg,ntype,data0)

!  Parse and encode a WSPR message.

  parameter (MASK15=32767)
  character*22 msg
  character*12 call1,call2
  character grid4*4,grid6*6
  logical lbad1,lbad2
  integer*1 data0(11)
  integer nu(0:9)
  data nu/0,-1,1,0,-1,2,1,0,-1,1/

! Standard WSPR message (types 0 3 7 10 13 17 ... 60)
  i1=index(msg,' ')
  i2=index(msg,'/')
  i3=index(msg,'<')
  call1=msg(:i1-1)
  if(i1.lt.3 .or. i1.gt.7 .or. i2.gt.0 .or. i3.gt.0) go to 10
  grid4=msg(i1+1:i1+4)
  call packcall(call1,n1,lbad1)
  call packgrid(grid4,ng,lbad2)
  if(lbad1 .or. lbad2) go to 10
  ndbm=0
  read(msg(i1+5:),*) ndbm
  if(ndbm.lt.0) ndbm=0
  if(ndbm.gt.60) ndbm=60
  ndbm=ndbm+nu(mod(ndbm,10))
  n2=128*ng + (ndbm+64)
  call pack50(n1,n2,data0)
  ntype=ndbm
  go to 900

10 if(i2.ge.2 .and. i3.lt.1) then
     call packpfx(call1,n1,ng,nadd)
     ndbm=0
     read(msg(i1+1:),*) ndbm
     if(ndbm.lt.0) ndbm=0
     if(ndbm.gt.60) ndbm=60
     ndbm=ndbm+nu(mod(ndbm,10))
     ntype=ndbm + 1 + nadd
     n2=128*ng + ntype + 64
     call pack50(n1,n2,data0)
  else if(i3.eq.1) then
     i4=index(msg,'>')
     call1=msg(2:i4-1)
     call hash(call1,i4-2,ih)
     grid6=msg(i1+1:i1+6)
     call2=grid6(2:6)//grid6(1:1)//'      '
     call packcall(call2,n1,lbad1)
     ndbm=0
     read(msg(i1+8:),*) ndbm
     if(ndbm.lt.0) ndbm=0
     if(ndbm.gt.60) ndbm=60
     ndbm=ndbm+nu(mod(ndbm,10))
     ntype=-(ndbm+1)
     n2=128*ih + ntype + 64
     call pack50(n1,n2,data0)
  endif
  go to 900

900 continue
  return
end subroutine wqencode


code:
!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name:    wsprcode.f90
! Description:  This program provides examples of the source encoding,
!               convulsional error-control coding, bit and symbol ordering,
!               and synchronizing information contained in WSPR messages.
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------

program wsprcode

  parameter (NSYM=162)
  parameter (MAXSYM=176)
  character*22 msg,msg2
  integer*1 data0(13)
  integer*1 data1(13)
  integer*1 dat(206)
  integer*1 softsym(206)

! Define the sync vector:
  integer*1 sync(NSYM)
  data sync/                                      &
       1,1,0,0,0,0,0,0,1,0,0,0,1,1,1,0,0,0,1,0,   &
       0,1,0,1,1,1,1,0,0,0,0,0,0,0,1,0,0,1,0,1,   &
       0,0,0,0,0,0,1,0,1,1,0,0,1,1,0,1,0,0,0,1,   &
       1,0,1,0,0,0,0,1,1,0,1,0,1,0,1,0,1,0,0,1,   &
       0,0,1,0,1,1,0,0,0,1,1,0,1,0,1,0,0,0,1,0,   &
       0,0,0,0,1,0,0,1,0,0,1,1,1,0,1,1,0,0,1,1,   &
       0,1,0,0,0,1,1,1,0,0,0,0,0,1,0,1,0,0,1,1,   &
       0,0,0,0,0,0,0,1,1,0,1,0,1,1,0,0,0,1,1,0,   &
       0,0/

! Metric table for decoding from soft symbols
  integer mettab(0:255,0:1)
  data mettab/                                            &
         5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   &
         5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   &
         5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   &
         5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   &
         5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   &
         5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   &
         5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   &
         5,   5,   5,   5,   5,   5,   5,   5,   5,   4,   &
         4,   4,   4,   4,   4,   4,   4,   4,   4,   4,   &
         4,   4,   4,   4,   4,   4,   4,   4,   4,   4,   &
         3,   3,   3,   3,   3,   3,   3,   3,   3,   2,   &
         2,   2,   2,   2,   1,   1,   1,   1,   0,   0,   &
        -1,  -1,  -1,  -2,  -2,  -3,  -4,  -4,  -5,  -6,   &
        -7,  -7,  -8,  -9, -10, -11, -12, -12, -13, -14,   &
       -15, -16, -17, -17, -18, -19, -20, -21, -22, -22,   &
       -23, -24, -25, -26, -26, -27, -28, -29, -30, -30,   &
       -31, -32, -33, -33, -34, -35, -36, -36, -37, -38,   &
       -38, -39, -40, -41, -41, -42, -43, -43, -44, -45,   &
       -45, -46, -47, -47, -48, -49, -49, -50, -51, -51,   &
       -52, -53, -53, -54, -54, -55, -56, -56, -57, -57,   &
       -58, -59, -59, -60, -60, -61, -62, -62, -62, -63,   &
       -64, -64, -65, -65, -66, -67, -67, -67, -68, -69,   &
       -69, -70, -70, -71, -72, -72, -72, -72, -73, -74,   &
       -75, -75, -75, -77, -76, -76, -78, -78, -80, -81,   &
       -80, -79, -83, -82, -81, -82, -82, -83, -84, -84,   &
       -84, -87, -86, -87, -88, -89, -89, -89, -88, -87,   &
       -86, -87, -84, -84, -84, -83, -82, -82, -81, -82,   &
       -83, -79, -80, -81, -80, -78, -78, -76, -76, -77,   &
       -75, -75, -75, -74, -73, -72, -72, -72, -72, -71,   &
       -70, -70, -69, -69, -68, -67, -67, -67, -66, -65,   &
       -65, -64, -64, -63, -62, -62, -62, -61, -60, -60,   &
       -59, -59, -58, -57, -57, -56, -56, -55, -54, -54,   &
       -53, -53, -52, -51, -51, -50, -49, -49, -48, -47,   &
       -47, -46, -45, -45, -44, -43, -43, -42, -41, -41,   &
       -40, -39, -38, -38, -37, -36, -36, -35, -34, -33,   &
       -33, -32, -31, -30, -30, -29, -28, -27, -26, -26,   &
       -25, -24, -23, -22, -22, -21, -20, -19, -18, -17,   &
       -17, -16, -15, -14, -13, -12, -12, -11, -10,  -9,   &
        -8,  -7,  -7,  -6,  -5,  -4,  -4,  -3,  -2,  -2,   &
        -1,  -1,  -1,   0,   0,   1,   1,   1,   1,   2,   &
         2,   2,   2,   2,   3,   3,   3,   3,   3,   3,   &
         3,   3,   3,   4,   4,   4,   4,   4,   4,   4,   &
         4,   4,   4,   4,   4,   4,   4,   4,   4,   4,   &
         4,   4,   4,   4,   5,   5,   5,   5,   5,   5,   &
         5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   &
         5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   &
         5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   &
         5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   &
         5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   &
         5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   &
         5,   5,   5,   5,   5,   5,   5,   5,   5,   5,   &
         5,   5/

! Get command-line argument(s)
  nargs=iargc()
  if(nargs.ne.1) then
     print*,'Usage: WSPRcode "message"'
     go to 999
  endif
  call getarg(1,msg)                             !Get message from command line
  write(*,1000) msg
1000 format('Message: ',a22)

  nbits=50+31               !User bits=50, constraint length=32
  nbytes=(nbits+7)/8
  ndelta=50
  limit=20000

  data0=0
  call wqencode(msg,ntype0,data0)             !Source encoding
  write(*,1002) data0(1:7),data0(1:6),data0(7)/64
1002 format(/'Source-encoded message, 50 bits:'/'Hex:    ',7z3.2/    &
          'Binary: ',6b9.8,b3.2)

  call encode232(data0,nbytes,dat,MAXSYM)     !Convolutional encoding
  call inter_mept(dat,1)                      !Interleaving

  write(*,1004)
1004 format(/'Data symbols:')
  write(*,1006) (dat(i),i=1,NSYM)
1006 format(5x,30i2)

  write(*,1008)
1008 format(/'Sync symbols:')
  write(*,1006) (sync(i),i=1,NSYM)

  write(*,1010)
1010 format(/'Channel symbols:')
  write(*,1006) (2*dat(i)+sync(i),i=1,NSYM)

  call inter_mept(dat,-1)                     !Remove interleaving
  softsym=-dat                                !Simulate soft symbols

! Call the sequential (Fano algorithm) decoder
  call fano232(softsym,nbits,mettab,ndelta,limit,data1,ncycles,metric,nerr)
  call wqdecode(data1,msg2,ntype1)

  write(*,1020) msg2,ntype1
1020 format(/'Decoded message: ',a22,'   ntype:',i3)

999 end program wsprcode

include 'wspr_old_subs.f90'


code:
CC = gcc
#CC = clang-3.5
FC = gfortran

CFLAGS= -I/usr/include -Wall -Wno-missing-braces -O3 -ffast-math
LDFLAGS = -L/usr/lib
FFLAGS =  -O2 -Wall -Wno-conversion
LIBS = -lfftw3f -lm

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

all:    wsprd wsprsim wsprd_exp

DEPS =  wsprsim_utils.h wsprd_utils.h fano.h jelinek.h nhash.h
OBJS1 = wsprd.o wsprsim_utils.o wsprd_utils.o tab.o fano.o jelinek.o nhash.o
wsprd: $(OBJS1)
	$(CC) -o $@ $^ $(CFLAGS) $(LDFLAGS) $(LIBS)

OBJS2 = wsprsim.o wsprsim_utils.o wsprd_utils.o tab.o fano.o nhash.o
wsprsim: $(OBJS2)
	$(CC) -o $@ $^ $(CFLAGS) $(LDFLAGS) $(LIBS)

clean:
	rm *.o wsprd wsprsim

path: lib/wsprd/Makefile


code:
CC = gcc
#CC = clang
FC = gfortran

FFLAGS = -O2 -Wall -Wno-conversion
CFLAGS= -Wall -Wno-missing-braces -O2
#LDFLAGS = -L/JTSDK/fftw3f
LIBS = c:/JTSDK/fftw3f/libfftw3-3.dll -lm

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

all:    wsprd

DEPS = fano.h
OBJS1 = wsprd.o wsprd_utils.o fano.o tab.o nhash.o
wsprd: $(OBJS1)
	$(CC) -o $@ $^ $(CFLAGS) $(LDFLAGS) $(LIBS)

clean:
	rm *.o wsprd


code:
CC = gcc
FC = gfortran

FFLAGS = -O2 -Wall -Wno-conversion
CFLAGS= -I/JTSDK/fftw3f -Wall -Wno-missing-braces -O2
#LDFLAGS = -L/JTSDK/fftw3f -Wl,--stack,4000000
LDFLAGS = -L/JTSDK/fftw3f
LIBS = c:/JTSDK/fftw3f/libfftw3-3.dll -lm

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

all:    wsprd.exe wsprsim.exe wsprd_exp.exe

DEPS =  wsprsim_utils.h wsprd_utils.h fano.h jelinek.h nhash.h

OBJS1 = wsprd.o wsprsim_utils.o wsprd_utils.o tab.o fano.o jelinek.o nhash.o
wsprd.exe: $(OBJS1)
	$(CC) -o $@ $^ $(CFLAGS) $(LDFLAGS) $(LIBS)

OBJS2 = wsprsim.o wsprsim_utils.o wsprd_utils.o tab.o fano.o nhash.o
wsprsim.exe: $(OBJS2)
	$(CC) -o $@ $^ $(CFLAGS) $(LDFLAGS) $(LIBS)

OBJS3 = wsprd_exp.o wsprsim_utils.o wsprd_utils.o tab.o fano.o jelinek.o \
      nhash.o
wsprd_exp.exe: $(OBJS3)
	$(CC) -o $@ $^ $(CFLAGS) $(LDFLAGS) $(LIBS)

clean:
	rm *.o wsprd.exe wsprsim.exe wsprd_exp.exe


code:
wsprd is a decoder for K1JT's Weak Signal Propagation Reporter (WSPR) mode.

The program is written in C and is a command-line program that reads from a
.c2 file or .wav file and writes output to the console. It is used by WSJT-X
for wspr-mode decoding.

USAGE:
       wsprd [options...] infile

OPTIONS:
       -a <path> path to writeable data files, default="."
       -c write .c2 file at the end of the first pass
       -e x (x is transceiver dial frequency error in Hz)
       -f x (x is transceiver dial frequency in MHz)
       -H do not use (or update) the hash table
       -m decode wspr-15 .wav file
       -q quick mode - doesn't dig deep for weak signals
       -s single pass mode, no subtraction (same as original wsprd)
       -v verbose mode (shows dupes)
       -w wideband mode - decode signals within +/- 150 Hz of center
       -z x (x is fano metric table bias, default is 0.42)

infile can be either .wav or .c2

e.g.
./wsprd -wf 14.0956 140709_2258.wav

Note that for .c2 files, the frequency within the file overrides the command
line value.

FEATURES:
By default, wsprd reports signals that are within +/- 110 Hz of the
subband center frequency. The wideband option (-w) extends this to +/- 150 Hz.

wsprd maintains a hashtable and will decode all three types of wspr
messages. An option (-H) is available to turn off use of the hashtable.

The symbols are decoded using Phil Karn's sequential decoder routine,
fano.c.

NOTES:
This program attempts to maximize the number of successful decodes per transmit
interval by trying to decode virtually every peak in the averaged spectrum.
The program also implements two-pass decoding, whereby signals that are successfully
decoded are subtracted one-by-one during the first decoding pass. Then, the
decoder is run again. In many cases the subtraction process will uncover signals
that can then be successfully decoded on the second pass.

There will be occasional duplicate decodes when two closely spaced
peaks come from the same signal. The program removes dupes based on callsign
and frequency. Two decodes that have the same callsign and estimated frequencies
that are within 1 Hz will be treated as decodes of the same signal. This
dupechecking is turned off with the -v flag.

path: lib/wsprd/README


code:
/*
 This file is part of wsprd.

 File name: fano.c

 Description: Soft decision Fano sequential decoder for K=32 r=1/2
 convolutional code.

 Copyright 1994, Phil Karn, KA9Q
 Minor modifications by Joe Taylor, K1JT
*/

#define	LL 1	                // Select Layland-Lushbaugh code
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "fano.h"

struct node {
  unsigned long encstate;	// Encoder state of next node
  long gamma;		        // Cumulative metric to this node
  int metrics[4];		// Metrics indexed by all possible tx syms
  int tm[2];		        // Sorted metrics for current hypotheses
  int i;			// Current branch being tested
};

// Convolutional coding polynomials. All are rate 1/2, K=32
#ifdef	NASA_STANDARD
/* "NASA standard" code by Massey & Costello
 * Nonsystematic, quick look-in, dmin=11, dfree=23
 * used on Pioneer 10-12, Helios A,B
 */
#define	POLY1	0xbbef6bb7
#define	POLY2	0xbbef6bb5
#endif

#ifdef	MJ
/* Massey-Johannesson code
 * Nonsystematic, quick look-in, dmin=13, dfree>=23
 * Purported to be more computationally efficient than Massey-Costello
 */
#define	POLY1	0xb840a20f
#define POLY2	0xb840a20d
#endif

#ifdef	LL
/* Layland-Lushbaugh code
 * Nonsystematic, non-quick look-in, dmin=?, dfree=?
 */
#define	POLY1	0xf2d05351
#define	POLY2	0xe4613c47
#endif

/* Convolutionally encode a packet. The input data bytes are read
 * high bit first and the encoded packet is written into 'symbols',
 * one symbol per byte. The first symbol is generated from POLY1,
 * the second from POLY2.
 *
 * Storing only one symbol per byte uses more space, but it is faster
 * and easier than trying to pack them more compactly.
 */
int encode(
	   unsigned char *symbols,	// Output buffer, 2*nbytes*8
	   unsigned char *data,		// Input buffer, nbytes
	   unsigned int nbytes)		// Number of bytes in data
{
  unsigned long encstate;
  int sym;
  int i;

  encstate = 0;
  while(nbytes-- != 0) {
    for(i=7;i>=0;i--) {
      encstate = (encstate << 1) | ((*data >> i) & 1);
      ENCODE(sym,encstate);
      *symbols++ = sym >> 1;
      *symbols++ = sym & 1;
    }
    data++;
  }
  return 0;
}

/* Decode packet with the Fano algorithm.
 * Return 0 on success, -1 on timeout
 */
int fano(
	 unsigned int  *metric,	   // Final path metric (returned value)
	 unsigned int  *cycles,	   // Cycle count (returned value)
	 unsigned int  *maxnp,     // Progress before timeout (returned value)
	 unsigned char *data,	   // Decoded output data
	 unsigned char *symbols,   // Raw deinterleaved input symbols
	 unsigned int nbits,	   // Number of output bits
	 int mettab[2][256],	   // Metric table, [sent sym][rx symbol]
	 int delta,		   // Threshold adjust parameter
	 unsigned int maxcycles)   // Decoding timeout in cycles per bit
{
  struct node *nodes;		   // First node
  struct node *np;	           // Current node
  struct node *lastnode;	   // Last node
  struct node *tail;		   // First node of tail
  int t;			   // Threshold
  int  m0,m1;
  int ngamma;
  unsigned int lsym;
  unsigned int i;

  if((nodes = (struct node *)malloc((nbits+1)*sizeof(struct node))) == NULL) {
    printf("malloc failed\n");
    return 0;
  }
  lastnode = &nodes[nbits-1];
  tail = &nodes[nbits-31];
  *maxnp = 0;

/* Compute all possible branch metrics for each symbol pair
 * This is the only place we actually look at the raw input symbols
 */
  for(np=nodes;np <= lastnode;np++) {
    np->metrics[0] = mettab[0][symbols[0]] + mettab[0][symbols[1]];
    np->metrics[1] = mettab[0][symbols[0]] + mettab[1][symbols[1]];
    np->metrics[2] = mettab[1][symbols[0]] + mettab[0][symbols[1]];
    np->metrics[3] = mettab[1][symbols[0]] + mettab[1][symbols[1]];
    symbols += 2;
  }
  np = nodes;
  np->encstate = 0;

// Compute and sort branch metrics from root node */
  ENCODE(lsym,np->encstate);	// 0-branch (LSB is 0)
  m0 = np->metrics[lsym];

/* Now do the 1-branch. To save another ENCODE call here and
 * inside the loop, we assume that both polynomials are odd,
 * providing complementary pairs of branch symbols.

 * This code should be modified if a systematic code were used.
 */

  m1 = np->metrics[3^lsym];
  if(m0 > m1) {
    np->tm[0] = m0;                             // 0-branch has better metric
    np->tm[1] = m1;
  } else {
    np->tm[0] = m1;                             // 1-branch is better
    np->tm[1] = m0;
    np->encstate++;	                        // Set low bit
  }
  np->i = 0;	                                // Start with best branch
  maxcycles *= nbits;
  np->gamma = t = 0;

  // Start the Fano decoder
  for(i=1;i <= maxcycles;i++) {
    if((int)(np-nodes) > (int)*maxnp) *maxnp=(int)(np-nodes);
#ifdef	debug
    printf("k=%ld, g=%ld, t=%d, m[%d]=%d, maxnp=%d, encstate=%lx\n",
	   np-nodes,np->gamma,t,np->i,np->tm[np->i],*maxnp,np->encstate);
#endif
// Look forward */
    ngamma = np->gamma + np->tm[np->i];
    if(ngamma >= t) {
      if(np->gamma < t + delta) {               // Node is acceptable
	/* First time we've visited this node;
	 * Tighten threshold.
	 *
	 * This loop could be replaced with
	 *   t += delta * ((ngamma - t)/delta);
	 * but the multiply and divide are slower.
	 */
	while(ngamma >= t + delta) t += delta;
      }
      np[1].gamma = ngamma;                     // Move forward
      np[1].encstate = np->encstate << 1;
      if( ++np == (lastnode+1) ) {
	break;	                                // Done!
      }

      /* Compute and sort metrics, starting with the
       * zero branch
       */
      ENCODE(lsym,np->encstate);
      if(np >= tail) {
	/* The tail must be all zeroes, so don't
	 * bother computing the 1-branches here.
	 */
	np->tm[0] = np->metrics[lsym];
      } else {
	m0 = np->metrics[lsym];
	m1 = np->metrics[3^lsym];
	if(m0 > m1) {
	  np->tm[0] = m0;                       // 0-branch is better
	  np->tm[1] = m1;
	} else {
	  np->tm[0] = m1;                       // 1-branch is better
	  np->tm[1] = m0;
	  np->encstate++;	                // Set low bit
	}
      }
      np->i = 0;	                        // Start with best branch
      continue;
    }
    // Threshold violated, can't go forward
    for(;;) {                                   // Look backward
      if(np == nodes || np[-1].gamma < t) {
	/* Can't back up either.
	 * Relax threshold and and look
	 * forward again to better branch.
	 */
	t -= delta;
	if(np->i != 0) {
	  np->i = 0;
	  np->encstate ^= 1;
	}
	break;
      }
      // Back up
      if(--np < tail && np->i != 1) {
	np->i++;                          // Search next best branch
	np->encstate ^= 1;
	break;
      }                                   // else keep looking back
    }
  }
  *metric =  np->gamma;	                  // Return the final path metric

  // Copy decoded data to user's buffer
  nbits >>= 3;
  np = &nodes[7];
  while(nbits-- != 0) {
    *data++ = np->encstate;
    np += 8;
  }
  *cycles = i+1;

  free(nodes);
  if(i >= maxcycles) return -1;	          // Decoder timed out
  return 0;		                  // Successful completion
}


code:
/*
 This file is part of wsprd.

 File name: fano.h

 Description: Header file for sequential Fano decoder.

 Copyright 1994, Phil Karn, KA9Q
 Minor modifications by Joe Taylor, K1JT
*/

#ifndef FANO_H
#define FANO_H

int fano(unsigned int *metric, unsigned int *cycles, unsigned int *maxnp,
	unsigned char *data,unsigned char *symbols, unsigned int nbits,
	 int mettab[2][256],int delta,unsigned int maxcycles);

int encode(unsigned char *symbols,unsigned char *data,unsigned int nbytes);

extern unsigned char Partab[];

/* Convolutional encoder macro. Takes the encoder state, generates
 * a rate 1/2 symbol pair and stores it in 'sym'. The symbol generated from
 * POLY1 goes into the 2-bit of sym, and the symbol generated from POLY2
 * goes into the 1-bit.
 */
#define	ENCODE(sym,encstate){\
unsigned long _tmp;\
\
_tmp = (encstate) & POLY1;\
_tmp ^= _tmp >> 16;\
(sym) = Partab[(_tmp ^ (_tmp >> 8)) & 0xff] << 1;\
_tmp = (encstate) & POLY2;\
_tmp ^= _tmp >> 16;\
(sym) |= Partab[(_tmp ^ (_tmp >> 8)) & 0xff];\
}

#endif