diff --git a/libm65/CMakeLists.txt b/libm65/CMakeLists.txt index b296c99f5..f80199df4 100644 --- a/libm65/CMakeLists.txt +++ b/libm65/CMakeLists.txt @@ -62,8 +62,10 @@ set (FSRCS decode0.f90 decode1a.f90 deep65.f90 + demod64a.f90 display.f90 dpol.f90 + extract.f90 four2a.f90 ftninit.f90 ftnquit.f90 @@ -94,10 +96,8 @@ set (FSRCS dcoord.f decode65b.f deg2grid.f - demod64a.f dot.f encode65.f - extract.F f77_wisdom.f fchisq.f fil6521.f diff --git a/libm65/demod64a.f b/libm65/demod64a.f deleted file mode 100644 index b5fbd6087..000000000 --- a/libm65/demod64a.f +++ /dev/null @@ -1,73 +0,0 @@ - subroutine demod64a(s3,nadd,mrsym,mrprob, - + mr2sym,mr2prob,ntest,nlow) - -C Demodulate the 64-bin spectra for each of 63 symbols in a frame. - -C Parameters -C nadd number of spectra already summed -C mrsym most reliable symbol value -C mr2sym second most likely symbol value -C mrprob probability that mrsym was the transmitted value -C mr2prob probability that mr2sym was the transmitted value - - implicit real*8 (a-h,o-z) - real*4 s3(64,63) - real*8 fs(64) - integer mrsym(63),mrprob(63),mr2sym(63),mr2prob(63) - common/mrscom/ mrs(63),mrs2(63) - - afac=1.1 * float(nadd)**0.64 - scale=255.999 - -C Compute average spectral value - sum=0. - do j=1,63 - do i=1,64 - sum=sum+s3(i,j) - enddo - enddo - ave=sum/(64.*63.) - i1=1 !Silence warning - i2=1 - -C Compute probabilities for most reliable symbol values - do j=1,63 - s1=-1.e30 - fsum=0. - do i=1,64 - x=min(afac*s3(i,j)/ave,50.d0) - fs(i)=exp(x) - fsum=fsum+fs(i) - if(s3(i,j).gt.s1) then - s1=s3(i,j) - i1=i !Most reliable - endif - enddo - - s2=-1.e30 - do i=1,64 - if(i.ne.i1 .and. s3(i,j).gt.s2) then - s2=s3(i,j) - i2=i !Second most reliable - endif - enddo - p1=fs(i1)/fsum !Normalized probabilities - p2=fs(i2)/fsum - mrsym(j)=i1-1 - mr2sym(j)=i2-1 - mrprob(j)=scale*p1 - mr2prob(j)=scale*p2 - mrs(j)=i1 - mrs2(j)=i2 - enddo - - sum=0. - nlow=0 - do j=1,63 - sum=sum+mrprob(j) - if(mrprob(j).le.5) nlow=nlow+1 - enddo - ntest=sum/63 - - return - end diff --git a/libm65/demod64a.f90 b/libm65/demod64a.f90 new file mode 100644 index 000000000..0f37095a5 --- /dev/null +++ b/libm65/demod64a.f90 @@ -0,0 +1,72 @@ +subroutine demod64a(s3,nadd,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow) + +! Demodulate the 64-bin spectra for each of 63 symbols in a frame. + +! Parameters +! nadd number of spectra already summed +! mrsym most reliable symbol value +! mr2sym second most likely symbol value +! mrprob probability that mrsym was the transmitted value +! mr2prob probability that mr2sym was the transmitted value + + implicit real*8 (a-h,o-z) + real*4 s3(64,63) + real*8 fs(64) + integer mrsym(63),mrprob(63),mr2sym(63),mr2prob(63) + common/mrscom/ mrs(63),mrs2(63) + + afac=1.1 * float(nadd)**0.64 + scale=255.999 + +! Compute average spectral value + sum=0. + do j=1,63 + do i=1,64 + sum=sum+s3(i,j) + enddo + enddo + ave=sum/(64.*63.) + i1=1 !Silence warning + i2=1 + +! Compute probabilities for most reliable symbol values + do j=1,63 + s1=-1.e30 + fsum=0. + do i=1,64 + x=min(afac*s3(i,j)/ave,50.d0) + fs(i)=exp(x) + fsum=fsum+fs(i) + if(s3(i,j).gt.s1) then + s1=s3(i,j) + i1=i !Most reliable + endif + enddo + + s2=-1.e30 + do i=1,64 + if(i.ne.i1 .and. s3(i,j).gt.s2) then + s2=s3(i,j) + i2=i !Second most reliable + endif + enddo + p1=fs(i1)/fsum !Normalized probabilities + p2=fs(i2)/fsum + mrsym(j)=i1-1 + mr2sym(j)=i2-1 + mrprob(j)=scale*p1 + mr2prob(j)=scale*p2 + mrs(j)=i1 + mrs2(j)=i2 + enddo + + sum=0. + nlow=0 + do j=1,63 + sum=sum+mrprob(j) + if(mrprob(j).le.5) nlow=nlow+1 + enddo + ntest=sum/63 + + return +end subroutine demod64a diff --git a/libm65/extract.F b/libm65/extract.F deleted file mode 100644 index 04313c241..000000000 --- a/libm65/extract.F +++ /dev/null @@ -1,109 +0,0 @@ - subroutine extract(s3,nadd,ncount,nhist,decoded,ltext) - - real s3(64,63) - real tmp(4032) - character decoded*22 - integer era(51),dat4(12),indx(64) - integer mrsym(63),mr2sym(63),mrprob(63),mr2prob(63) - logical first,ltext - data first/.true./,nsec1/0/ - save - - nfail=0 - 1 continue -! call timer('demod64a',0) - call demod64a(s3,nadd,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow) -! call timer('demod64a',1) - if(ntest.lt.50 .or. nlow.gt.20) then - ncount=-999 !Flag bad data - go to 900 - endif - call chkhist(mrsym,nhist,ipk) - - if(nhist.ge.20) then - nfail=nfail+1 - call pctile(s3,tmp,4032,50,base) ! ### or, use ave from demod64a - do j=1,63 - s3(ipk,j)=base - enddo - if(nfail.gt.30) then - decoded=' ' - ncount=-1 - go to 900 - endif - go to 1 - endif - - call graycode(mrsym,63,-1) - call interleave63(mrsym,-1) - call interleave63(mrprob,-1) - - ndec=1 - nemax=30 !Was 200 (30) - maxe=8 - xlambda=13.0 !Was 12 - - if(ndec.eq.1) then - call graycode(mr2sym,63,-1) - call interleave63(mr2sym,-1) - call interleave63(mr2prob,-1) - - nsec1=nsec1+1 - write(22,rec=1) nsec1,xlambda,maxe,200, - + mrsym,mrprob,mr2sym,mr2prob - call flush(22) -! call timer('kvasd ',0) -#ifdef UNIX - iret=system('./kvasd -q > dev_null') -#else - iret=system('kvasd -q > dev_null') -#endif -! call timer('kvasd ',1) - if(iret.ne.0) then - if(first) write(*,1000) iret - 1000 format('Error in KV decoder, or no KV decoder present.'/ - + 'Return code:',i8,'. Will use BM algorithm.') - ndec=0 - first=.false. - go to 20 - endif - - read(22,rec=2) nsec2,ncount,dat4 - j=nsec2 !Silence compiler warning - decoded=' ' - ltext=.false. - if(ncount.ge.0) then - call unpackmsg(dat4,decoded) !Unpack the user message - if(iand(dat4(10),8).ne.0) ltext=.true. - do i=2,12 - if(dat4(i).ne.dat4(1)) go to 20 - enddo - write(13,*) 'Bad decode?',nhist,nfail,ipk, - + ' ',dat4,decoded - ncount=-1 !Suppress supposedly bogus decodes - decoded=' ' - endif - endif - 20 if(ndec.eq.0) then - call indexx(63,mrprob,indx) - do i=1,nemax - j=indx(i) - if(mrprob(j).gt.120) then - ne2=i-1 - go to 2 - endif - era(i)=j-1 - enddo - ne2=nemax - 2 decoded=' ' - do nerase=0,ne2,2 - call rs_decode(mrsym,era,nerase,dat4,ncount) - if(ncount.ge.0) then - call unpackmsg(dat4,decoded) - go to 900 - endif - enddo - endif - - 900 return - end diff --git a/libm65/extract.f90 b/libm65/extract.f90 new file mode 100644 index 000000000..6e3457ee0 --- /dev/null +++ b/libm65/extract.f90 @@ -0,0 +1,107 @@ +subroutine extract(s3,nadd,ncount,nhist,decoded,ltext) + + real s3(64,63) + real tmp(4032) + character decoded*22 + integer era(51),dat4(12),indx(64) + integer mrsym(63),mr2sym(63),mrprob(63),mr2prob(63) + logical first,ltext + data first/.true./,nsec1/0/ + save + + nfail=0 +1 continue +! call timer('demod64a',0) + call demod64a(s3,nadd,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow) +! call timer('demod64a',1) + if(ntest.lt.50 .or. nlow.gt.20) then + ncount=-999 !Flag bad data + go to 900 + endif + call chkhist(mrsym,nhist,ipk) + + if(nhist.ge.20) then + nfail=nfail+1 + call pctile(s3,tmp,4032,50,base) ! ### or, use ave from demod64a + do j=1,63 + s3(ipk,j)=base + enddo + if(nfail.gt.30) then + decoded=' ' + ncount=-1 + go to 900 + endif + go to 1 + endif + + call graycode(mrsym,63,-1) + call interleave63(mrsym,-1) + call interleave63(mrprob,-1) + + ndec=1 + nemax=30 !Was 200 (30) + maxe=8 + xlambda=13.0 !Was 12 + + if(ndec.eq.1) then + call graycode(mr2sym,63,-1) + call interleave63(mr2sym,-1) + call interleave63(mr2prob,-1) + + nsec1=nsec1+1 + write(22,rec=1) nsec1,xlambda,maxe,200,mrsym,mrprob,mr2sym,mr2prob + call flush(22) +! call timer('kvasd ',0) +!#ifdef UNIX +! iret=system('./kvasd -q > dev_null') +!#else + iret=system('kvasd -q > dev_null') +!#endif +! call timer('kvasd ',1) + if(iret.ne.0) then + if(first) write(*,1000) iret +1000 format('Error in KV decoder, or no KV decoder present.'/ & + 'Return code:',i8,'. Will use BM algorithm.') + ndec=0 + first=.false. + go to 20 + endif + + read(22,rec=2) nsec2,ncount,dat4 + j=nsec2 !Silence compiler warning + decoded=' ' + ltext=.false. + if(ncount.ge.0) then + call unpackmsg(dat4,decoded) !Unpack the user message + if(iand(dat4(10),8).ne.0) ltext=.true. + do i=2,12 + if(dat4(i).ne.dat4(1)) go to 20 + enddo + write(13,*) 'Bad decode?',nhist,nfail,ipk,' ',dat4,decoded + ncount=-1 !Suppress supposedly bogus decodes + decoded=' ' + endif + endif +20 if(ndec.eq.0) then + call indexx(63,mrprob,indx) + do i=1,nemax + j=indx(i) + if(mrprob(j).gt.120) then + ne2=i-1 + go to 2 + endif + era(i)=j-1 + enddo + ne2=nemax +2 decoded=' ' + do nerase=0,ne2,2 + call rs_decode(mrsym,era,nerase,dat4,ncount) + if(ncount.ge.0) then + call unpackmsg(dat4,decoded) + go to 900 + endif + enddo + endif + +900 return +end subroutine extract