diff --git a/CMakeLists.txt b/CMakeLists.txt index a7bae6402..c3dbf6f47 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -592,6 +592,7 @@ set (wsjt_FSRCS lib/superfox/sfox_demod.f90 lib/superfox/sfox_clo.f90 lib/superfox/sym_prob.f90 + lib/superfox/getpp3.f90 lib/superfox/ran1.f90 ) @@ -637,6 +638,7 @@ set (wsjt_CSRCS lib/superfox/encode_rs.c lib/superfox/decode_rs.c lib/superfox/rs_sf.c + lib/superfox/ftrsd3.c ${ldpc_CSRCS} ${qra_CSRCS} ) diff --git a/lib/superfox/ftrsd3.f90 b/lib/superfox/ftrsd3.f90 new file mode 100644 index 000000000..0f8dffc04 --- /dev/null +++ b/lib/superfox/ftrsd3.f90 @@ -0,0 +1,189 @@ +subroutine ftrsd3(rxdat,rxprob,rxdat2,rxprob2,ntrials0,correct,param,ntry) + +! Soft-decision decoder for Reed-Solomon codes. + +! 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 + + use sfox_mod + + integer, dimension(0:NN-1) :: rxdat,rxprob,rxdat2,rxprob2,workdat, & + correct,indexes + integer rxdat(0:NN-1) !Hard-decision symbol values + integer rxprob(0:NN-1) !Probabilities that rxdat values are correct + integer rxdat2(0:NN-1) !Second most probable symbol values + integer rxprob2(0:NN-1) !Probabilities that rxdat2 values are correct + integer workdat(0:NN-1) !Work array + integer correct(0:NN-1) !Corrected codeword + integer indexes(0:NN-1) !For sorting probabilities + integer probs(0:NN-1) !Temp array for sorting probabilities + integer thresh0(0:NN-1) !Temp array for thresholds + integer era_pos(0:NN-KK-1) !Index values for erasures + integer*8 nseed,ir !No unsigned int in Fortran + integer pass,tmp + + integer perr(0:7,0:7) + data perr/ 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/ + + ntrials=ntrials0 + nhard=0 + nhard_min=32768 + nsoft=0 + nsoft_min=32768 + ntotal=0 + ntotal_min=32768 + nera_best=0 + nsym=nn + + do i=0,NN-1 + indexes(i)=i + probs(i)=rxprob(i) + enddo + + do pass=1,nsym-1 + do k=0,nsym-pass-1 + if(probs(k).lt.probs(k+1)) then + probs(k)=probs(k+1) + probs(k+1)=tmp + tmp=indexes(k) + indexes(k)=indexes(k+1) + indexes(k+1)=tmp + enddo + enddo + enddo + + era_pos=0 + numera=0 + workdat=rxdat + call rs_decode_sf(workdat,era_pos,numera,nerr) !Call the decoder + + if(nerr.ge.0) then +! Hard-decision decoding succeeded. Save codeword and some parameters. + nhard=count(workdat.ne.rxdat) + correct=workdat + param(0)=0 + param(1)=nhard + param(2)=0 + param(3)=0 + param(4)=0 + param(5)=0 + param(7)=1000*1000 !??? + ntry=0 + return + endif + +! 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 + ncandidates=0 + nsum=0 + do i=0,NN-1 + nsum=nsum+rxprob(i) + j=indexes(NN-1-i) + ratio=(float)rxprob2(j)/((float)rxprob(j)+0.01) + ii=7.999*ratio + jj=(NN-1-i)/8 + thresh0(i)=1.3*perr(jj)(ii) + enddo + if(nsum.le.0) return + + pp1=0. + pp2=0. + do k=1,ntrials + era_pos=0 + workdat=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 + do i=0,NN-1 + j=indexes(126-i) + thresh=thresh0(i) + +! Generate a random number ir, 0 <= ir <= 100 (see POSIX.1-2001 example). + nseed=nseed*1103515245 + 12345 + ir=mod(nseed/65536),32768) + ir=(100*ir)/32768 + nseed=iand(ir,4294967295) + + if((ir.lt.thresh ) .and. numera.lt.(NN-KK)) then + era_pos(numera)=j + numera=numera+1 + endif + enddo + +! nerr=decode_rs_int(rs,workdat,era_pos,numera,0); + call rs_decode_sf(workdat,era_pos,numera,nerr) !Call the decoder + + if( nerr.ge.0) then + ! 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,127) of synchronized symbol spectra. + ncandidates=ncandidates+1 + nhard=0 + nsoft=0 + do i=0,NN-1 + if(workdat(i).ne. rxdat(i)) then + nhard=nhard+1; + if(workdat(i) .ne. rxdat2(i)) nsoft=nsoft+rxprob(i) + endif + enddo + nsoft=(NN-1)*nsoft/nsum + ntotal=nsoft+nhard + + pp=0. + call getpp3(s3,workdat,pp) + if(pp.gt.pp1) then + pp2=pp1 + pp1=pp + nsoft_min=nsoft + nhard_min=nhard + ntotal_min=ntotal + correct=workdat + nera_best=numera + ntry=k + else + if(pp.gt.pp2 .and. pp.ne.pp1) pp2=pp + endif + if(nhard_min <= 41 && ntotal_min <= 71) exit !### New values ### + enddo + if(k.eq.ntrials) ntry=k + enddo + + param(0)=ncandidates + param(1)=nhard_min + param(2)=nsoft_min + param(3)=nera_best +! param(4)= pp1 > 0 ? 1000.0*pp2/pp1 : 1000.0 + param(5)=ntotal_min + param(6)=ntry + param(7)=1000.0*pp2 + param(8)=1000.0*pp1 + if(param(0).eq.0) param(2)=-1 + + return +end subroutine ftrsd3 + diff --git a/lib/superfox/sfox_demod.f90 b/lib/superfox/sfox_demod.f90 index d7b993bf1..d81886304 100644 --- a/lib/superfox/sfox_demod.f90 +++ b/lib/superfox/sfox_demod.f90 @@ -4,15 +4,14 @@ subroutine sfox_demod(crcvd,f,t,s3,chansym) complex crcvd(NMAX) !Signal as received complex c(0:NSPS-1) !Work array, one symbol long real s(0:NQ-1) !Power spectrum - real s3(0:NQ-1,0:ND-1) !Symbol spectra + real s3(0:NQ-1,0:NN-1) !Symbol spectra integer chansym(NN) !Hard-decision symbol values integer ipk(1) i0=nint(12000.0*t) df=12000.0/NSPS j0=nint(f/df)-NQ/2 - chansym=0 - do n=1,ND !Loop over all symbols + do n=1,NN !Loop over all symbols ib=n*NSPS + i0 if(n.gt.ND1) ib=(NS+n)*NSPS + i0 ia=ib-NSPS+1 diff --git a/lib/superfox/sfox_gen.f90 b/lib/superfox/sfox_gen.f90 index c2cfea961..81e7f4e5a 100644 --- a/lib/superfox/sfox_gen.f90 +++ b/lib/superfox/sfox_gen.f90 @@ -4,7 +4,7 @@ subroutine sfox_gen(idat,f0,fsample,syncwidth,cdat) ! include "sfox_params.f90" complex cdat(NMAX) !Generated complex waveform complex w,wstep - integer idat(ND) + integer idat(NN) twopi=8.0*atan(1.0) tsync=NS*NSPS/fsample diff --git a/lib/superfox/sfox_mod.f90 b/lib/superfox/sfox_mod.f90 index 4c84a9ae4..d9abb1752 100644 --- a/lib/superfox/sfox_mod.f90 +++ b/lib/superfox/sfox_mod.f90 @@ -1,7 +1,7 @@ module sfox_mod parameter (NMAX=15*12000) !Samples in iwave (180,000) - integer MM,NQ,NN,KK,ND1,ND2,ND,NFZ,NSPS,NS,NSYNC,NZ,NFFT,NFFT1 + integer MM,NQ,NN,KK,ND1,ND2,NFZ,NSPS,NS,NSYNC,NZ,NFFT,NFFT1 contains subroutine sfox_init(mm0,nn0,kk0,itu,fspread,delay) @@ -17,21 +17,20 @@ contains MM=mm0 !Bits per symbol NQ=2**MM !Q, number of MFSK tones - NN=nn0 !Channel symbols, before puncture - KK=kk0 !Information symbols, before puncture + NN=nn0 !Number of channel symbols + KK=kk0 !Information symbols ND1=25 !Data symbols before sync ND2=NN-ND1 !Data symbols after sync - ND=NN !Total data symbols NFZ=3 !First zero tsync=2.0 - jsps=nint((12.8-tsync)*12000.0/ND) + jsps=nint((12.8-tsync)*12000.0/NN) iloc=minloc(abs(isps-jsps)) NSPS=isps(iloc(1)) !Samples per symbol NS=nint(tsync*12000.0/NSPS) if(mod(NS,2).eq.1) NS=NS+1 NSYNC=NS*NSPS !Samples in sync waveform - NZ=NSPS*(ND+NS) !Samples in full Tx waveform + NZ=NSPS*(NN+NS) !Samples in full Tx waveform NFFT=32768 !Length of FFT for sync waveform NFFT1=2*NSPS !Length of FFTs for symbol spectra diff --git a/lib/superfox/sfoxtest.f90 b/lib/superfox/sfoxtest.f90 index b624792c0..f32b58df9 100644 --- a/lib/superfox/sfoxtest.f90 +++ b/lib/superfox/sfoxtest.f90 @@ -6,6 +6,7 @@ program sfoxtest use sfox_mod type(hdr) h !Header for .wav file integer*2 iwave(NMAX) !Generated i*2 waveform + integer nparam(0:7) real*4 xnoise(NMAX) !Random noise real*4 dat(NMAX) !Generated real data complex cdat(NMAX) !Generated complex waveform @@ -13,7 +14,7 @@ program sfoxtest complex cnoise(NMAX) !Complex noise complex crcvd(NMAX) !Signal as received real a(3) - real, allocatable :: s3(:,:) !Symbol spectra: will be s3(NQ,ND) + real, allocatable :: s3(:,:) !Symbol spectra: will be s3(NQ,NN) integer, allocatable :: msg0(:) !Information symbols integer, allocatable :: parsym(:) !Parity symbols integer, allocatable :: chansym0(:) !Encoded data, 7-bit integers @@ -23,6 +24,8 @@ program sfoxtest integer, allocatable :: rxprob(:) integer, allocatable :: rxdat2(:) integer, allocatable :: rxprob2(:) + integer, allocatable :: correct(:) + character fname*17,arg*12,itu*2 nargs=iargc() @@ -75,16 +78,17 @@ program sfoxtest ' MaxErr:',i3,' tsync:',f4.1,' TxT:',f5.1/) ! Allocate storage for arrays that depend on code parameters. - allocate(s3(0:NQ-1,0:ND-1)) + allocate(s3(0:NQ-1,0:NN-1)) allocate(msg0(1:KK)) allocate(parsym(1:NN-KK)) allocate(chansym0(1:NN)) allocate(chansym(1:NN)) allocate(iera(1:NN)) - allocate(rxdat(0:ND-1)) - allocate(rxprob(0:ND-1)) - allocate(rxdat2(0:ND-1)) - allocate(rxprob2(0:ND-1)) + allocate(rxdat(0:NN-1)) + allocate(rxprob(0:NN-1)) + allocate(rxdat2(0:NN-1)) + allocate(rxprob2(0:NN-1)) + allocate(correct(0:NN-1)) rms=100. fsample=12000.0 !Sample rate (Hz) @@ -138,7 +142,8 @@ program sfoxtest f1=f0 if(f0.eq.0.0) then f1=1500.0 + 200.0*(ran1(idummy)-0.5) - xdt=0.6*(ran1(idummy)-0.5) +! xdt=0.6*(ran1(idummy)-0.5) + xdt=0.3*ran1(idummy) call sfox_gen(chansym0,f1,fsample,syncwidth,cdat,clo) endif @@ -157,10 +162,15 @@ program sfoxtest call sfox_sync(crcvd,clo,nv,f,t) ferr=f-f1 terr=t-xdt + igoodsync=0 if(abs(ferr).lt.baud/2.0 .and. abs(terr).lt.tsym/8.0) then + igoodsync=1 ngoodsync=ngoodsync+1 sqt=sqt + terr*terr sqf=sqf + ferr*ferr +! else +! write(*,3003) ferr,terr +!3003 format('Sync failed:',f8.1,f8.3) endif a=0. @@ -169,22 +179,25 @@ program sfoxtest f=1500.0 call sfox_demod(crcvd,f,t,s3,chansym) !Get s3 and hard symbol values call sym_prob(s3,rxdat,rxprob,rxdat2,rxprob2) - chansym(1:ND)=rxdat !### TEMPORARY ? ### -! do j=0,ND-1 -! do j=0,5 -! write(*,3001) j,chansym(1+j),rxdat(j),rxprob(j),rxdat2(j),rxprob2(j) -!3001 format('prob'i5,5i8) -! enddo + if(igoodsync.eq.1) then + do j=0,NN-1 + if(chansym(1+j).ne.rxdat(j)) write(*,3001) xdt,j,chansym(1+j), & + rxdat(j),rxprob(j),rxdat2(j),rxprob2(j) +3001 format(f7.3,i5,5i8) + enddo + endif nera=0 chansym=mod(chansym,nq) !Enforce 0 to nq-1 nharderr=count(chansym.ne.chansym0) !Count hard errors -! nhard2=count(rxdat.ne.chansym0(1:ND)) !Count hard errors +! nhard2=count(rxdat.ne.chansym0(1:NN)) !Count hard errors ! print*,'A',nharderr,nhard2 ntot=ntot+nharderr nworst=max(nworst,nharderr) - call rs_decode_sf(chansym,iera,nera,nfixed) !Call the decoder - + call rs_decode_sf(rxdat,iera,nera,nfixed) !Call the BM decoder + ntrials=1000 +! call ftrsd3(rxdat,rxprob,rxdat2,rxprob2,ntrials,nparam,correct,ntry) + if(iand(nv,1).ne.0) then fname='000000_000001.wav' write(fname(8:13),'(i6.6)') ifile