diff --git a/lib/fsk4hf/ft8b.f90 b/lib/fsk4hf/ft8b.f90 index 6d762cecb..c26ec0520 100644 --- a/lib/fsk4hf/ft8b.f90 +++ b/lib/fsk4hf/ft8b.f90 @@ -10,16 +10,17 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, & character*6 mygrid6 logical bcontest real a(5) - real s1(0:7,ND),s2(0:7,NN) + real s1(0:7,ND),s2(0:7,NN),s1sort(8*ND) real ps(0:7) - real bmeta(3*ND),bmetap(3*ND) - real llr(3*ND),llra(3*ND),llr0(3*ND),llrap(3*ND) !Soft symbols + real bmeta(3*ND),bmetb(3*ND),bmetap(3*ND) + real llr(3*ND),llra(3*ND),llr0(3*ND),llr1(3*ND),llrap(3*ND) !Soft symbols real dd0(15*12000) integer*1 decoded(KK),apmask(3*ND),cw(3*ND) integer*1 msgbits(KK) integer apsym(KK) integer mcq(28),mde(28),mrrr(16),m73(16),mrr73(16) integer itone(NN) + integer indxs1(8*ND) integer icos7(0:6),ip(1) 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 @@ -27,6 +28,7 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, & complex ctwk(32) complex csymb(32) logical first,newdat,lsubtract,lapon,nagain + equivalence (s1,s1sort) data icos7/2,5,6,0,4,1,3/ data mcq/1,1,1,1,1,0,1,0,0,0,0,0,1,0,0,0,0,0,1,1,0,0,0,1,1,0,0,1/ data mrrr/0,1,1,1,1,1,1,0,1,1,0,0,1,1,1,1/ @@ -153,10 +155,15 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, & s1(0:7,j)=s2(0:7,k) enddo + call indexx(s1sort,8*ND,indxs1) + xmeds1=s1sort(indxs1(nint(0.5*8*ND))) + s1=s1/xmeds1 + do j=1,ND i4=3*j-2 i2=3*j-1 i1=3*j +! Max amplitude ps=s1(0:7,j) r1=max(ps(1),ps(3),ps(5),ps(7))-max(ps(0),ps(2),ps(4),ps(6)) r2=max(ps(2),ps(3),ps(6),ps(7))-max(ps(0),ps(1),ps(4),ps(5)) @@ -167,6 +174,34 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, & bmetap(i4)=r4 bmetap(i2)=r2 bmetap(i1)=r1 +! Max log + ps=log(ps) + r1=max(ps(1),ps(3),ps(5),ps(7))-max(ps(0),ps(2),ps(4),ps(6)) + r2=max(ps(2),ps(3),ps(6),ps(7))-max(ps(0),ps(1),ps(4),ps(5)) + r4=max(ps(4),ps(5),ps(6),ps(7))-max(ps(0),ps(1),ps(2),ps(3)) +! Metric for Cauchy noise +! r1=log(ps(1)**3+ps(3)**3+ps(5)**3+ps(7)**3)- & +! log(ps(0)**3+ps(2)**3+ps(4)**3+ps(6)**3) +! r2=log(ps(2)**3+ps(3)**3+ps(6)**3+ps(7)**3)- & +! log(ps(0)**3+ps(1)**3+ps(4)**3+ps(5)**3) +! r4=log(ps(4)**3+ps(5)**3+ps(6)**3+ps(7)**3)- & +! log(ps(0)**3+ps(1)**3+ps(2)**3+ps(3)**3) +! Metric for AWGN, no fading +! bscale=2.5 +! b0=bessi0(bscale*ps(0)) +! b1=bessi0(bscale*ps(1)) +! b2=bessi0(bscale*ps(2)) +! b3=bessi0(bscale*ps(3)) +! b4=bessi0(bscale*ps(4)) +! b5=bessi0(bscale*ps(5)) +! b6=bessi0(bscale*ps(6)) +! b7=bessi0(bscale*ps(7)) +! r1=log(b1+b3+b5+b7)-log(b0+b2+b4+b6) +! r2=log(b2+b3+b6+b7)-log(b0+b1+b4+b5) +! r4=log(b4+b5+b6+b7)-log(b0+b1+b2+b3) + bmetb(i4)=r4 + bmetb(i2)=r2 + bmetb(i1)=r1 if(nQSOProgress .eq. 0 .or. nQSOProgress .eq. 5) then ! When bits 88:115 are set as ap bits, bit 115 lives in symbol 39 along @@ -209,12 +244,14 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, & enddo call normalizebmet(bmeta,3*ND) + call normalizebmet(bmetb,3*ND) call normalizebmet(bmetap,3*ND) - ss=0.84 - llr0=2.0*bmeta/(ss*ss) - llra=2.0*bmetap/(ss*ss) ! llr's for use with ap - apmag=4.0 + scalefac=2.83 + llr0=scalefac*bmeta + llr1=scalefac*bmetb + llra=scalefac*bmetap ! llr's for use with ap + apmag=scalefac*(maxval(abs(bmetap))*1.01) ! pass # !------------------------------ @@ -227,35 +264,36 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, & ! 7 ap pass 4, etc. if(lapon) then - npasses=3+nappasses(nQSOProgress) + npasses=4+nappasses(nQSOProgress) else - npasses=3 + npasses=4 endif do ipass=1,npasses llr=llr0 - if(ipass.eq.2) llr(1:24)=0. - if(ipass.eq.3) llr(1:48)=0. - if(ipass.le.3) then + if(ipass.eq.2) llr=llr1 + if(ipass.eq.3) llr(1:24)=0. + if(ipass.eq.4) llr(1:48)=0. + if(ipass.le.4) then apmask=0 llrap=llr iaptype=0 endif - if(ipass .gt. 3) then - iaptype=naptypes(nQSOProgress,ipass-3) + if(ipass .gt. 4) then + iaptype=naptypes(nQSOProgress,ipass-4) if(iaptype.ge.3 .and. (abs(f1-nfqso).gt.napwid .and. abs(f1-nftx).gt.napwid) ) cycle if(iaptype.eq.1 .or. iaptype.eq.2 ) then ! AP,???,??? apmask=0 apmask(88:115)=1 ! first 28 bits are AP apmask(144)=1 ! not free text llrap=llr - if(iaptype.eq.1) llrap(88:115)=apmag*mcq/ss - if(iaptype.eq.2) llrap(88:115)=apmag*apsym(1:28)/ss + if(iaptype.eq.1) llrap(88:115)=apmag*mcq + if(iaptype.eq.2) llrap(88:115)=apmag*apsym(1:28) llrap(116:117)=llra(116:117) llrap(142:143)=llra(142:143) - llrap(144)=-apmag/ss + llrap(144)=-apmag endif if(iaptype.eq.3) then ! mycall, dxcall, ??? apmask=0 @@ -263,8 +301,8 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, & apmask(116:143)=1 ! hiscall apmask(144)=1 ! not free text llrap=llr - llrap(88:143)=apmag*apsym(1:56)/ss - llrap(144)=-apmag/ss + llrap(88:143)=apmag*apsym(1:56) + llrap(144)=-apmag endif if(iaptype.eq.4 .or. iaptype.eq.5 .or. iaptype.eq.6) then apmask=0 @@ -272,10 +310,10 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, & apmask(116:143)=1 ! hiscall apmask(144:159)=1 ! RRR or 73 or RR73 llrap=llr - llrap(88:143)=apmag*apsym(1:56)/ss - if(iaptype.eq.4) llrap(144:159)=apmag*mrrr/ss - if(iaptype.eq.5) llrap(144:159)=apmag*m73/ss - if(iaptype.eq.6) llrap(144:159)=apmag*mrr73/ss + llrap(88:143)=apmag*apsym(1:56) + if(iaptype.eq.4) llrap(144:159)=apmag*mrrr + if(iaptype.eq.5) llrap(144:159)=apmag*m73 + if(iaptype.eq.6) llrap(144:159)=apmag*mrr73 endif if(iaptype.eq.7) then ! ???, dxcall, ??? apmask=0 @@ -283,8 +321,8 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, & apmask(144)=1 ! not free text llrap=llr llrap(115)=llra(115) - llrap(116:143)=apmag*apsym(29:56)/ss - llrap(144)=-apmag/ss + llrap(116:143)=apmag*apsym(29:56) + llrap(144)=-apmag endif endif @@ -297,7 +335,7 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, & if(ndepth.eq.3 .and. nharderrors.lt.0) then ndeep=3 if(abs(nfqso-f1).le.napwid .or. abs(nftx-f1).le.napwid) then - if((ipass.eq.2 .or. ipass.eq.3) .and. .not.nagain) then + if((ipass.eq.3 .or. ipass.eq.4) .and. .not.nagain) then ndeep=3 else ndeep=4 @@ -315,8 +353,8 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid, & if(any(decoded(73:75).ne.0)) cycle !Reject if any of the 3 extra bits are nonzero if(nharderrors.ge.0 .and. nharderrors+dmin.lt.60.0 .and. & .not.(sync.lt.2.0 .and. nharderrors.gt.35) .and. & - .not.(ipass.gt.1 .and. nharderrors.gt.39) .and. & - .not.(ipass.eq.3 .and. nharderrors.gt.30) & + .not.(ipass.gt.2 .and. nharderrors.gt.39) .and. & + .not.(ipass.eq.4 .and. nharderrors.gt.30) & ) then call chkcrc12a(decoded,nbadcrc) else @@ -363,3 +401,28 @@ subroutine normalizebmet(bmet,n) bmet=bmet/bmetsig return end subroutine normalizebmet + + +FUNCTION bessi0(x) +! From Numerical Recipes + REAL bessi0,x + DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y + SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 + DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,3.5156229d0,3.0899424d0,1.2067492d0, & + 0.2659732d0,0.360768d-1,0.45813d-2/ + DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,0.1328592d-1, & + 0.225319d-2,-0.157565d-2,0.916281d-2,-0.2057706d-1, & + 0.2635537d-1,-0.1647633d-1,0.392377d-2/ + + if (abs(x).lt.3.75) then + y=(x/3.75)**2 + bessi0=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))) + else + ax=abs(x) + y=3.75/ax + bessi0=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4 & + +y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))) + endif + return +end function bessi0 +