From e90cc846bbf0d36055462d6350a1d1ce5d7c2ee2 Mon Sep 17 00:00:00 2001 From: Steven Franke Date: Sat, 2 Dec 2017 16:04:51 +0000 Subject: [PATCH] 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 | 2 +- lib/decode65b.f90 | 3 +- lib/decoder.f90 | 2 +- lib/extract.f90 | 123 +++++++++++++++++++----- lib/ftrsd/ftrsdap.c | 227 ++++++++++++++++++++++++++++++++++++++++++++ lib/jt65_decode.f90 | 5 +- 6 files changed, 333 insertions(+), 29 deletions(-) create mode 100644 lib/ftrsd/ftrsdap.c diff --git a/CMakeLists.txt b/CMakeLists.txt index 66d19471b..5ac010fea 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/lib/decode65b.f90 b/lib/decode65b.f90 index 3ff6dc017..215b6bacd 100644 --- a/lib/decode65b.f90 +++ b/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 diff --git a/lib/decoder.f90 b/lib/decoder.f90 index d4ecbe5e0..5b2e12a8d 100644 --- a/lib/decoder.f90 +++ b/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='# ' diff --git a/lib/extract.f90 b/lib/extract.f90 index 2a0a97205..cfd2c45d6 100644 --- a/lib/extract.f90 +++ b/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, & - ltext,nft,qual) + mycall_12,hiscall_12,hisgrid,nQSOProgress,ljt65apon, & + 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) - call timer('ftrsd ',0) - param=0 - call ftrsd2(mrsym,mrprob,mr2sym,mr2prob,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 + npass=1 ! if ap decoding is disabled + if(ljt65apon .and. len_trim(mycall).gt.0) then + npass=1+nappasses(nQSOProgress) +!write(*,*) 'ap is on: ',npass-1,'ap passes of types ',naptypes(nQSOProgress,:) 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) diff --git a/lib/ftrsd/ftrsdap.c b/lib/ftrsd/ftrsdap.c new file mode 100644 index 000000000..b2b85ddd9 --- /dev/null +++ b/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 +#include +#include +#include +#include +#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=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= 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; +} diff --git a/lib/jt65_decode.f90 b/lib/jt65_decode.f90 index f0dcb62af..2ae90f9df 100644 --- a/lib/jt65_decode.f90 +++ b/lib/jt65_decode.f90 @@ -490,8 +490,9 @@ contains enddo nadd=nsum*ismo - call extract(s3,nadd,mode65,ntrials,naggressive,ndepth,nflip,mycall, & - hiscall,hisgrid,nexp_decode,ncount,nhist,decoded,ltext,nft,qual) + call extract(s3c,nadd,mode65,ntrials,naggressive,ndepth,nflip,mycall, & + hiscall,hisgrid,nQSOProgress,ljt65apon,nexp_decode,ncount,nhist, & + avemsg,ltext,nftt,qual) if(nftt.eq.1) then nsmo=ismo param(9)=nsmo