Implement AP decoding for JT65 (VHF/UHF/Microwave only).

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@8277 ab8295b8-cf94-4d9e-aec4-7959e3be5d79

CMakeLists.txt

@@ -604,7 +604,7 @@ set (qra_CSRCS
 
 set (wsjt_CSRCS
   ${ka9q_CSRCS}
-  lib/ftrsd/ftrsd2.c
+  lib/ftrsd/ftrsdap.c
   lib/sgran.c
   lib/golay24_table.c
   lib/gran.c

lib/decode65b.f90

@@ -21,7 +21,8 @@ subroutine decode65b(s2,nflip,nadd,mode65,ntrials,naggressive,ndepth,      &
   enddo
 
   call extract(s3,nadd,mode65,ntrials,naggressive,ndepth,nflip,mycall,   &
-       hiscall,hisgrid,nexp_decode,ncount,nhist,decoded,ltext,nft,qual) 
+      hiscall,hisgrid,nQSOProgress,ljt65apon,nexp_decode,ncount,         &
+      nhist,decoded,ltext,nft,qual) 
 
 ! Suppress "birdie messages" and other garbage decodes:
   if(decoded(1:7).eq.'000AAA ') ncount=-1

lib/decoder.f90

@@ -366,7 +366,7 @@ contains
           endif
           nap=ishft(ft,-2)
           if(nap.ne.0) then
-            write(cflags(1:3),'(a1,i2.2)') 'a',nap 
+            write(cflags(1:3),'(a1,i1)') 'a',nap 
           endif
        endif
        csync='# '

lib/extract.f90

@@ -1,6 +1,7 @@
 subroutine extract(s3,nadd,mode65,ntrials,naggressive,ndepth,nflip,     &
-     mycall_12,hiscall_12,hisgrid,nexp_decode,ncount,nhist,decoded,     &
+     mycall_12,hiscall_12,hisgrid,nQSOProgress,ljt65apon,               &
-     ltext,nft,qual)
+     nexp_decode,ncount,      &
+     nhist,decoded,ltext,nft,qual)
 
 ! Input:
 !   s3       64-point spectra for each of 63 data symbols
@@ -20,19 +21,73 @@ subroutine extract(s3,nadd,mode65,ntrials,naggressive,ndepth,nflip,     &
   use timer_module, only: timer
 
   real s3(64,63)
-  character decoded*22
+  character decoded*22, apmessage*22
   character*12 mycall_12,hiscall_12
   character*6 mycall,hiscall,hisgrid
+  character*6 mycall0,hiscall0,hisgrid0
+  integer apsymbols(5,12),ap(12)
+  integer nappasses(0:5)  ! the number of decoding passes to use for each QSO state
+  integer naptypes(0:5,4) ! (nQSOProgress, decoding pass)  maximum of 4 passes for now 
   integer dat4(12)
   integer mrsym(63),mr2sym(63),mrprob(63),mr2prob(63)
   integer correct(63),tmp(63)
-  logical ltext
+  logical first,ltext,ljt65apon
   common/chansyms65/correct
+  data first/.true./
   save
-
+  
   if(mode65.eq.-99) stop                   !Silence compiler warning
+  if(first) then
+
+! aptype
+!------------------------
+!   1        CQ     ???    ???
+!   2        MyCall ???    ???
+!   3        MyCall DxCall ???
+!   4        MyCall DxCall RRR
+!   5        MyCall DxCall 73
+
+     apsymbols=-1
+     nappasses=(/2,2,2,3,3,3/)
+     naptypes(0,1:4)=(/1,2,0,0/)
+     naptypes(1,1:4)=(/2,3,0,0/)
+     naptypes(2,1:4)=(/2,3,0,0/)
+     naptypes(3,1:4)=(/3,4,5,0/)
+     naptypes(4,1:4)=(/3,4,5,0/)
+     naptypes(5,1:4)=(/3,1,2,0/)  
+     first=.false.
+  endif
+
   mycall=mycall_12(1:6)
   hiscall=hiscall_12(1:6)
+! Fill apsymbols array
+  if(ljt65apon .and. (mycall.ne.mycall0 .or. hiscall.ne.hiscall0)) then 
+!write(*,*) 'initializing apsymbols '
+     apsymbols=-1
+     mycall0=mycall
+     hiscall0=hiscall
+     ap=-1
+     apsymbols(1,1:4)=(/62,32,32,49/) ! CQ
+     if(len_trim(mycall).gt.0) then
+        apmessage=mycall//" "//mycall//" RRR" 
+        call packmsg(apmessage,ap,itype,.false.)
+        if(itype.ne.1) ap=-1
+        apsymbols(2,1:4)=ap(1:4)
+!write(*,*) 'mycall symbols ',ap(1:4)
+        if(len_trim(hiscall).gt.0) then
+           apmessage=mycall//" "//hiscall//" RRR" 
+           call packmsg(apmessage,ap,itype,.false.)
+           if(itype.ne.1) ap=-1
+           apsymbols(3,1:9)=ap(1:9)
+           apsymbols(4,:)=ap
+           apmessage=mycall//" "//hiscall//" 73" 
+           call packmsg(apmessage,ap,itype,.false.)
+           if(itype.ne.1) ap=-1
+           apsymbols(5,:)=ap
+        endif
+     endif
+  endif
+   
   qual=0.
   nbirdie=20
   npct=50
@@ -71,28 +126,48 @@ subroutine extract(s3,nadd,mode65,ntrials,naggressive,ndepth,nflip,     &
   call graycode65(mr2sym,63,-1)      !Remove gray code and interleaving
   call interleave63(mr2sym,-1)       !from second-most-reliable symbols
   call interleave63(mr2prob,-1)
-  ntry=0
 
-!call gf64_osd(mrsym,mrprob,mr2sym,mr2prob,correct)
+  npass=1  ! if ap decoding is disabled  
-  call timer('ftrsd   ',0)
+  if(ljt65apon .and. len_trim(mycall).gt.0) then 
-  param=0
+     npass=1+nappasses(nQSOProgress)
-  call ftrsd2(mrsym,mrprob,mr2sym,mr2prob,ntrials,correct,param,ntry)
+!write(*,*) 'ap is on: ',npass-1,'ap passes of types ',naptypes(nQSOProgress,:)
-  call timer('ftrsd   ',1)
-  ncandidates=param(0)
-  nhard=param(1)
-  nsoft=param(2)
-  nerased=param(3)
-  rtt=0.001*param(4)
-  ntotal=param(5)
-  qual=0.001*param(7)
-  nd0=81
-  r0=0.87
-  if(naggressive.eq.10) then
-     nd0=83
-     r0=0.90
   endif
-  if(ntotal.le.nd0 .and. rtt.le.r0) nft=1
+  do ipass=1,npass
+     ap=-1
+     ntype=0
+     if(ipass.gt.1) then
+       ntype=naptypes(nQSOProgress,ipass-1)
+!write(*,*) 'ap pass, type ',ntype
+       ap=apsymbols(ntype,:)
+       if(count(ap.ge.0).eq.0) cycle  ! don't bother if all ap symbols are -1
+!write(*,'(12i3)') ap
+     endif
+     ntry=0
+     call timer('ftrsd   ',0)
+     param=0
+     call ftrsdap(mrsym,mrprob,mr2sym,mr2prob,ap,ntrials,correct,param,ntry)
+     call timer('ftrsd   ',1)
+     ncandidates=param(0)
+     nhard=param(1)
+     nsoft=param(2)
+     nerased=param(3)
+     rtt=0.001*param(4)
+     ntotal=param(5)
+     qual=0.001*param(7)
+     nd0=81
+     r0=0.87
+     if(naggressive.eq.10) then
+        nd0=83
+        r0=0.90
+     endif
 
+     if(ntotal.le.nd0 .and. rtt.le.r0) then
+        nft=1+ishft(ntype,2)
+     endif 
+  
+     if(nft.gt.0) exit
+  enddo
+!write(*,*) nft
   if(nft.eq.0 .and. iand(ndepth,32).eq.32) then
      qmin=2.0 - 0.1*naggressive
      call timer('hint65  ',0)

lib/ftrsd/ftrsdap.c

@@ -0,0 +1,227 @@
+/*
+ ftrsdap.c
+ 
+ A soft-decision decoder for the JT65 (63,12) Reed-Solomon code.
+ 
+ This decoding scheme is built around Phil Karn's Berlekamp-Massey
+ errors and erasures decoder. The approach is inspired by a number of
+ publications, including the stochastic Chase decoder described
+ in "Stochastic Chase Decoding of Reed-Solomon Codes", by Leroux et al.,
+ IEEE Communications Letters, Vol. 14, No. 9, September 2010 and
+ "Soft-Decision Decoding of Reed-Solomon Codes Using Successive Error-
+ and-Erasure Decoding," by Soo-Woong Lee and B. V. K. Vijaya Kumar.
+ 
+ Steve Franke K9AN and Joe Taylor K1JT
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <time.h>
+#include <string.h>
+#include "../ftrsd/rs2.h"
+
+static void *rs;
+void getpp_(int workdat[], float *pp);
+
+void ftrsdap_(int mrsym[], int mrprob[], int mr2sym[], int mr2prob[], 
+	     int ap[], int* ntrials0, int correct[], int param[], int ntry[])
+{
+  int rxdat[63], rxprob[63], rxdat2[63], rxprob2[63];
+  int workdat[63];
+  int indexes[63];
+  int era_pos[51];
+  int i, j, numera, nerr, nn=63;
+  int ntrials = *ntrials0;
+  int nhard=0,nhard_min=32768,nsoft=0,nsoft_min=32768;
+  int ntotal=0,ntotal_min=32768,ncandidates;
+  int nera_best=0;
+  float pp,pp1,pp2;
+  static unsigned int nseed;
+  
+// Power-percentage symbol metrics - composite gnnf/hf 
+  int perr[8][8] = {
+    { 4,      9,     11,     13,     14,     14,     15,     15},
+    { 2,     20,     20,     30,     40,     50,     50,     50},
+    { 7,     24,     27,     40,     50,     50,     50,     50},
+    {13,     25,     35,     46,     52,     70,     50,     50},
+    {17,     30,     42,     54,     55,     64,     71,     70},
+    {25,     39,     48,     57,     64,     66,     77,     77},
+    {32,     45,     54,     63,     66,     75,     78,     83},
+    {51,     58,     57,     66,     72,     77,     82,     86}};
+
+    
+// Initialize the KA9Q Reed-Solomon encoder/decoder
+  unsigned int symsize=6, gfpoly=0x43, fcr=3, prim=1, nroots=51;
+  rs=init_rs_int(symsize, gfpoly, fcr, prim, nroots, 0);
+
+// Reverse the received symbol vectors for BM decoder
+  for (i=0; i<63; i++) {
+    rxdat[i]=mrsym[62-i];
+    rxprob[i]=mrprob[62-i];
+    rxdat2[i]=mr2sym[62-i];
+    rxprob2[i]=mr2prob[62-i];
+  }
+
+// Set ap symbols and ap mask
+  for (i=0; i<12; i++) {
+    if(ap[i]>=0) {
+      rxdat[11-i]=ap[i];
+      rxprob2[11-i]=-1;
+    }
+  }
+ 
+// Sort rxprob to find indexes of the least reliable symbols
+  int k, pass, tmp, nsym=63;
+  int probs[63];
+  for (i=0; i<63; i++) {
+    indexes[i]=i;
+    probs[i]=rxprob[i];
+  }
+  for (pass = 1; pass <= nsym-1; pass++) {
+    for (k = 0; k < nsym - pass; k++) {
+      if( probs[k] < probs[k+1] ) {
+        tmp = probs[k];
+        probs[k] = probs[k+1];
+        probs[k+1] = tmp;
+        tmp = indexes[k];
+        indexes[k] = indexes[k+1];
+        indexes[k+1] = tmp;
+      }
+    }
+  }
+  
+// See if we can decode using BM HDD, and calculate the syndrome vector.
+  memset(era_pos,0,51*sizeof(int));
+  numera=0;
+  memcpy(workdat,rxdat,sizeof(rxdat));
+  nerr=decode_rs_int(rs,workdat,era_pos,numera,1);
+  if( nerr >= 0 ) {
+    // Hard-decision decoding succeeded.  Save codeword and some parameters.
+    nhard=0;
+    for (i=0; i<63; i++) {
+      if( workdat[i] != rxdat[i] ) nhard=nhard+1;
+    }
+    memcpy(correct,workdat,63*sizeof(int));
+    param[0]=0;
+    param[1]=nhard;
+    param[2]=0;
+    param[3]=0;
+    param[4]=0;
+    param[5]=0;
+    param[7]=1000*1000;
+    ntry[0]=0;
+    return;
+  }
+
+/*
+Hard-decision decoding failed.  Try the FT soft-decision method.
+Generate random erasure-locator vectors and see if any of them
+decode. This will generate a list of "candidate" codewords.  The
+soft distance between each candidate codeword and the received 
+word is estimated by finding the largest (pp1) and second-largest 
+(pp2) outputs from a synchronized filter-bank operating on the 
+symbol spectra, and using these to decide which candidate 
+codeword is "best".
+*/
+
+  nseed=1;                                 //Seed for random numbers
+  float ratio;
+  int thresh, nsum;
+  int thresh0[63];
+  ncandidates=0;
+  nsum=0;
+  int ii,jj;
+  for (i=0; i<nn; i++) {
+    nsum=nsum+rxprob[i];
+    j = indexes[62-i];
+    if( rxprob2[j]>=0 ) {
+      ratio = (float)rxprob2[j]/((float)rxprob[j]+0.01);
+      ii = 7.999*ratio;
+      jj = (62-i)/8;
+      thresh0[i] = 1.3*perr[ii][jj];
+    } else {
+      thresh0[i] = 0.0;
+    }
+//printf("%d %d %d\n",i,j,rxdat[i]);
+  }
+
+  if(nsum<=0) return;
+
+  pp1=0.0;
+  pp2=0.0;
+  for (k=1; k<=ntrials; k++) {
+    memset(era_pos,0,51*sizeof(int));
+    memcpy(workdat,rxdat,sizeof(rxdat));
+
+/* 
+Mark a subset of the symbols as erasures.
+Run through the ranked symbols, starting with the worst, i=0.
+NB: j is the symbol-vector index of the symbol with rank i.
+*/
+    numera=0;
+    for (i=0; i<nn; i++) {
+      j = indexes[62-i];
+      thresh=thresh0[i];
+      long int ir;
+
+// Generate a random number ir, 0 <= ir < 100 (see POSIX.1-2001 example).
+      nseed = nseed * 1103515245 + 12345;
+      ir = (unsigned)(nseed/65536) % 32768;
+      ir = (100*ir)/32768;
+
+      if((ir < thresh ) && numera < 51) {
+        era_pos[numera]=j;
+        numera=numera+1;
+      }
+    }
+
+    nerr=decode_rs_int(rs,workdat,era_pos,numera,0);        
+    if( nerr >= 0 ) {
+      // We have a candidate codeword.  Find its hard and soft distance from
+      // the received word.  Also find pp1 and pp2 from the full array 
+      // s3(64,63) of synchronized symbol spectra.
+      ncandidates=ncandidates+1;
+      nhard=0;
+      nsoft=0;
+      for (i=0; i<63; i++) {
+        if(workdat[i] != rxdat[i]) {
+          nhard=nhard+1;
+          if(workdat[i] != rxdat2[i]) {
+            nsoft=nsoft+rxprob[i];
+          }
+        }
+      }
+      nsoft=63*nsoft/nsum;
+      ntotal=nsoft+nhard;
+
+      getpp_(workdat,&pp);
+      if(pp>pp1) {
+        pp2=pp1;
+        pp1=pp;
+        nsoft_min=nsoft;
+        nhard_min=nhard;
+        ntotal_min=ntotal;
+        memcpy(correct,workdat,63*sizeof(int));
+        nera_best=numera;
+        ntry[0]=k;
+      } else {
+        if(pp>pp2 && pp!=pp1) pp2=pp;
+      }
+      if(nhard_min <= 41 && ntotal_min <= 71) break;
+    }
+    if(k == ntrials) ntry[0]=k;
+  }
+  
+  param[0]=ncandidates;
+  param[1]=nhard_min;
+  param[2]=nsoft_min;
+  param[3]=nera_best;
+  param[4]=1000.0*pp2/pp1;
+  param[5]=ntotal_min;
+  param[6]=ntry[0];
+  param[7]=1000.0*pp2;
+  param[8]=1000.0*pp1;
+  if(param[0]==0) param[2]=-1;
+  return;
+}

lib/jt65_decode.f90

@@ -490,8 +490,9 @@ contains
        enddo
 
        nadd=nsum*ismo
-       call extract(s3,nadd,mode65,ntrials,naggressive,ndepth,nflip,mycall,   &
+       call extract(s3c,nadd,mode65,ntrials,naggressive,ndepth,nflip,mycall, &
-            hiscall,hisgrid,nexp_decode,ncount,nhist,decoded,ltext,nft,qual) 
+            hiscall,hisgrid,nQSOProgress,ljt65apon,nexp_decode,ncount,nhist, &
+            avemsg,ltext,nftt,qual)
        if(nftt.eq.1) then
           nsmo=ismo
           param(9)=nsmo