diff --git a/CMakeLists.txt b/CMakeLists.txt index 6a7f0b58b..ad8bdb970 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -521,6 +521,7 @@ set (wsjt_FSRCS lib/fsk4hf/msksoftsymw.f90 lib/ft8/osd174.f90 lib/fsk4hf/osd300.f90 + lib/fsk4hf/osdwspr.f90 lib/pctile.f90 lib/peakdt9.f90 lib/peakup.f90 @@ -565,6 +566,7 @@ set (wsjt_FSRCS lib/sync9f.f90 lib/sync9w.f90 lib/synciscat.f90 + lib/fsk4hf/tccsim.f90 lib/timf2.f90 lib/to_contest_msg.f90 lib/tweak1.f90 @@ -586,6 +588,8 @@ set (wsjt_FSRCS lib/fsk4hf/wspr_fsk8_wav.f90 lib/fsk4hf/wspr_fsk8_downsample.f90 lib/fsk4hf/wsprlfsim.f90 + lib/fsk4hf/wsprsimf.f90 + lib/fsk4hf/wspr_wav.f90 lib/wspr_downsample.f90 lib/zplot9.f90 ) @@ -1232,9 +1236,15 @@ target_link_libraries (ft8code wsjt_fort wsjt_cxx) add_executable (ft8sim lib/ft8/ft8sim.f90 wsjtx.rc) target_link_libraries (ft8sim wsjt_fort wsjt_cxx) +add_executable (tccsim lib/fsk4hf/tccsim.f90 wsjtx.rc) +target_link_libraries (tccsim wsjt_fort wsjt_cxx) + add_executable (wsprlfsim lib/fsk4hf/wsprlfsim.f90 wsjtx.rc) target_link_libraries (wsprlfsim wsjt_fort wsjt_cxx) +add_executable (wsprsimf lib/fsk4hf/wsprsimf.f90 wsjtx.rc) +target_link_libraries (wsprsimf wsjt_fort wsjt_cxx) + add_executable (wspr_fsk8d lib/fsk4hf/wspr_fsk8d.f90 wsjtx.rc) target_link_libraries (wspr_fsk8d wsjt_fort wsjt_cxx) diff --git a/lib/fsk4hf/encode4K25A.f90 b/lib/fsk4hf/encode4K25A.f90 new file mode 100644 index 000000000..cd1ae762a --- /dev/null +++ b/lib/fsk4hf/encode4K25A.f90 @@ -0,0 +1,56 @@ +subroutine encode4K25A(message,codeword) +! A (280,70) rate 1/4 tailbiting convolutional code using +! the "4K25A" polynomials from EbNaut website. +! Code is transparent, has constraint length 25, and has dmin=58 +character*10 g1,g2,g3,g4 +integer*1 codeword(280) +!integer*1 p1(25),p2(25),p3(25),p4(25) +integer*1 p1(16),p2(16),p3(16),p4(16) +integer*1 gg(100) +integer*1 gen(280,70) +integer*1 itmp(280) +integer*1 message(70) +logical first +data first/.true./ +data g1/"106042635"/ +data g2/"125445117"/ +data g3/"152646773"/ +data g4/"167561761"/ +!data p1/1,0,0,0,1,1,0,0,0,0,1,0,0,0,1,0,1,1,0,0,1,1,1,0,1/ +!data p2/1,0,1,0,1,0,1,1,0,0,1,0,0,1,0,1,0,0,1,0,0,1,1,1,1/ +!data p3/1,1,0,1,0,1,0,1,1,0,1,0,0,1,1,0,1,1,1,1,1,1,0,1,1/ +!data p4/1,1,1,0,1,1,1,1,0,1,1,1,0,0,0,1,1,1,1,1,1,0,0,0,1/ +data p1/1,0,1,0,1,1,0,0,1,1,0,1,1,1,1,1/ +data p2/1,0,1,1,0,1,0,0,1,1,1,1,1,0,0,1/ +data p3/1,1,0,0,1,0,1,1,0,1,1,1,0,0,1,1/ +data p4/1,1,1,0,1,1,0,1,1,1,1,0,0,1,0,1/ + +save first,gen + +if( first ) then ! fill the generator matrix + gg=0 +! gg(1:25)=p1 +! gg(26:50)=p2 +! gg(51:75)=p3 +! gg(76:100)=p4 + gg(1:16)=p1 + gg(17:32)=p2 + gg(33:48)=p3 + gg(49:64)=p4 + gen=0 +! gen(1:100,1)=gg(1:100) + gen(1:64,1)=gg(1:64) + do i=2,70 + gen(:,i)=cshift(gen(:,i-1),-4,1) + enddo + first=.false. +endif + +codeword=0 +do i=1,70 + if(message(i).eq.1) codeword=codeword+gen(:,i) +enddo +codeword=mod(codeword,2) + +return +end subroutine encode4K25A diff --git a/lib/fsk4hf/osdtbcc.f90 b/lib/fsk4hf/osdtbcc.f90 new file mode 100644 index 000000000..462fcc792 --- /dev/null +++ b/lib/fsk4hf/osdtbcc.f90 @@ -0,0 +1,372 @@ +subroutine osdtbcc(llr,apmask,ndeep,decoded,cw,nhardmin,dmin) +! +use iso_c_binding +parameter (N=280, K=70, L=16) + +integer*1 p1(L),p2(L),p3(L),p4(L) +integer*1 gg(100) + +integer*1 apmask(N),apmaskr(N) +integer*1 gen(K,N) +integer*1 genmrb(K,N),g2(N,K) +integer*1 temp(K),m0(K),me(K),mi(K),misub(K),e2sub(N-K),e2(N-K),ui(N-K) +integer*1 r2pat(N-K) +integer indices(N),nxor(N) +integer*1 cw(N),ce(N),c0(N),hdec(N) +integer*1 decoded(K) +integer indx(N) +real llr(N),rx(N),absrx(N) +logical first,reset +data first/.true./ +!data p1/1,0,0,0,1,1,0,0,0,0,1,0,0,0,1,0,1,1,0,0,1,1,1,0,1/ +!data p2/1,0,1,0,1,0,1,1,0,0,1,0,0,1,0,1,0,0,1,0,0,1,1,1,1/ +!data p3/1,1,0,1,0,1,0,1,1,0,1,0,0,1,1,0,1,1,1,1,1,1,0,1,1/ +!data p4/1,1,1,0,1,1,1,1,0,1,1,1,0,0,0,1,1,1,1,1,1,0,0,0,1/ +data p1/1,0,1,0,1,1,0,0,1,1,0,1,1,1,1,1/ +data p2/1,0,1,1,0,1,0,0,1,1,1,1,1,0,0,1/ +data p3/1,1,0,0,1,0,1,1,0,1,1,1,0,0,1,1/ +data p4/1,1,1,0,1,1,0,1,1,1,1,0,0,1,0,1/ + +save first,gen + +if( first ) then ! fill the generator matrix + gg=0 + gg(1:L)=p1 + gg(L+1:2*L)=p2 + gg(2*L+1:3*L)=p3 + gg(3*L+1:4*L)=p4 + gen=0 + gen(1,1:4*L)=gg(1:4*L) + do i=2,K + gen(i,:)=cshift(gen(i-1,:),-4) + enddo + first=.false. +endif + +! Re-order received vector to place systematic msg bits at the end. +rx=llr +apmaskr=apmask + +! Hard decisions on the received word. +hdec=0 +where(rx .ge. 0) hdec=1 + +! Use magnitude of received symbols as a measure of reliability. +absrx=abs(rx) +call indexx(absrx,N,indx) + +! Re-order the columns of the generator matrix in order of decreasing reliability. +do i=1,N + genmrb(1:K,i)=gen(1:K,indx(N+1-i)) + indices(i)=indx(N+1-i) +enddo + +! Do gaussian elimination to create a generator matrix with the most reliable +! received bits in positions 1:K in order of decreasing reliability (more or less). +do id=1,K ! diagonal element indices + do icol=id,K+20 ! The 20 is ad hoc - beware + iflag=0 + if( genmrb(id,icol) .eq. 1 ) then + iflag=1 + if( icol .ne. id ) then ! reorder column + temp(1:K)=genmrb(1:K,id) + genmrb(1:K,id)=genmrb(1:K,icol) + genmrb(1:K,icol)=temp(1:K) + itmp=indices(id) + indices(id)=indices(icol) + indices(icol)=itmp + endif + do ii=1,K + if( ii .ne. id .and. genmrb(ii,id) .eq. 1 ) then + genmrb(ii,1:N)=ieor(genmrb(ii,1:N),genmrb(id,1:N)) + endif + enddo + exit + endif + enddo +enddo + +g2=transpose(genmrb) + +! The hard decisions for the K MRB bits define the order 0 message, m0. +! Encode m0 using the modified generator matrix to find the "order 0" codeword. +! Flip various combinations of bits in m0 and re-encode to generate a list of +! codewords. Return the member of the list that has the smallest Euclidean +! distance to the received word. + +hdec=hdec(indices) ! hard decisions from received symbols +m0=hdec(1:K) ! zero'th order message +absrx=absrx(indices) +rx=rx(indices) +apmaskr=apmaskr(indices) + +call mrbencode(m0,c0,g2,N,K) +nxor=ieor(c0,hdec) +nhardmin=sum(nxor) +dmin=sum(nxor*absrx) + +cw=c0 +ntotal=0 +nrejected=0 + +if(ndeep.eq.0) goto 998 ! norder=0 +if(ndeep.gt.5) ndeep=5 +if( ndeep.eq. 1) then + nord=1 + npre1=0 + npre2=0 + nt=120 + ntheta=12 +elseif(ndeep.eq.2) then + nord=1 + npre1=1 + npre2=0 + nt=120 + ntheta=12 +elseif(ndeep.eq.3) then + nord=1 + npre1=1 + npre2=1 + nt=120 + ntheta=12 + ntau=15 +elseif(ndeep.eq.4) then + nord=2 + npre1=1 + npre2=0 + nt=120 + ntheta=12 + ntau=15 +elseif(ndeep.eq.5) then + nord=3 + npre1=1 + npre2=1 + nt=80 + ntheta=22 + ntau=16 +endif + +do iorder=1,nord + misub(1:K-iorder)=0 + misub(K-iorder+1:K)=1 + iflag=K-iorder+1 + do while(iflag .ge.0) + if(iorder.eq.nord .and. npre1.eq.0) then + iend=iflag + else + iend=1 + endif + do n1=iflag,iend,-1 + mi=misub + mi(n1)=1 + if(any(iand(apmaskr(1:K),mi).eq.1)) cycle + ntotal=ntotal+1 + me=ieor(m0,mi) + if(n1.eq.iflag) then + call mrbencode(me,ce,g2,N,K) + e2sub=ieor(ce(K+1:N),hdec(K+1:N)) + e2=e2sub + nd1Kpt=sum(e2sub(1:nt))+1 + d1=sum(ieor(me(1:K),hdec(1:K))*absrx(1:K)) + else + e2=ieor(e2sub,g2(K+1:N,n1)) + nd1Kpt=sum(e2(1:nt))+2 + endif + if(nd1Kpt .le. ntheta) then + call mrbencode(me,ce,g2,N,K) + nxor=ieor(ce,hdec) + if(n1.eq.iflag) then + dd=d1+sum(e2sub*absrx(K+1:N)) + else + dd=d1+ieor(ce(n1),hdec(n1))*absrx(n1)+sum(e2*absrx(K+1:N)) + endif + if( dd .lt. dmin ) then + dmin=dd + cw=ce + nhardmin=sum(nxor) + nd1Kptbest=nd1Kpt + endif + else + nrejected=nrejected+1 + endif + enddo +! Get the next test error pattern, iflag will go negative +! when the last pattern with weight iorder has been generated. + call nextpat(misub,k,iorder,iflag) + enddo +enddo + +if(npre2.eq.1) then + reset=.true. + ntotal=0 + do i1=K,1,-1 + do i2=i1-1,1,-1 + ntotal=ntotal+1 + mi(1:ntau)=ieor(g2(K+1:K+ntau,i1),g2(K+1:K+ntau,i2)) + call boxit(reset,mi(1:ntau),ntau,ntotal,i1,i2) + enddo + enddo + + ncount2=0 + ntotal2=0 + reset=.true. +! Now run through again and do the second pre-processing rule + misub(1:K-nord)=0 + misub(K-nord+1:K)=1 + iflag=K-nord+1 + do while(iflag .ge.0) + me=ieor(m0,misub) + call mrbencode(me,ce,g2,N,K) + e2sub=ieor(ce(K+1:N),hdec(K+1:N)) + do i2=0,ntau + ntotal2=ntotal2+1 + ui=0 + if(i2.gt.0) ui(i2)=1 + r2pat=ieor(e2sub,ui) +778 continue + call fetchit(reset,r2pat(1:ntau),ntau,in1,in2) + if(in1.gt.0.and.in2.gt.0) then + ncount2=ncount2+1 + mi=misub + mi(in1)=1 + mi(in2)=1 + if(sum(mi).lt.nord+npre1+npre2.or.any(iand(apmaskr(1:K),mi).eq.1)) cycle + me=ieor(m0,mi) + call mrbencode(me,ce,g2,N,K) + nxor=ieor(ce,hdec) + dd=sum(nxor*absrx) + if( dd .lt. dmin ) then + dmin=dd + cw=ce + nhardmin=sum(nxor) + endif + goto 778 + endif + enddo + call nextpat(misub,K,nord,iflag) + enddo +endif + +998 continue +! Re-order the codeword to place message bits at the end. +cw(indices)=cw +hdec(indices)=hdec +decoded=0 +return +end subroutine osdtbcc + +subroutine mrbencode(me,codeword,g2,N,K) +integer*1 me(K),codeword(N),g2(N,K) +! fast encoding for low-weight test patterns + codeword=0 + do i=1,K + if( me(i) .eq. 1 ) then + codeword=ieor(codeword,g2(1:N,i)) + endif + enddo +return +end subroutine mrbencode + +subroutine nextpat(mi,k,iorder,iflag) + integer*1 mi(k),ms(k) +! generate the next test error pattern + ind=-1 + do i=1,k-1 + if( mi(i).eq.0 .and. mi(i+1).eq.1) ind=i + enddo + if( ind .lt. 0 ) then ! no more patterns of this order + iflag=ind + return + endif + ms=0 + ms(1:ind-1)=mi(1:ind-1) + ms(ind)=1 + ms(ind+1)=0 + if( ind+1 .lt. k ) then + nz=iorder-sum(ms) + ms(k-nz+1:k)=1 + endif + mi=ms + do i=1,k ! iflag will point to the lowest-index 1 in mi + if(mi(i).eq.1) then + iflag=i + exit + endif + enddo + return +end subroutine nextpat + +subroutine boxit(reset,e2,ntau,npindex,i1,i2) + integer*1 e2(1:ntau) + integer indexes(4000,2),fp(0:525000),np(4000) + logical reset + common/boxes/indexes,fp,np + + if(reset) then + patterns=-1 + fp=-1 + np=-1 + sc=-1 + indexes=-1 + reset=.false. + endif + + indexes(npindex,1)=i1 + indexes(npindex,2)=i2 + ipat=0 + do i=1,ntau + if(e2(i).eq.1) then + ipat=ipat+ishft(1,ntau-i) + endif + enddo + + ip=fp(ipat) ! see what's currently stored in fp(ipat) + if(ip.eq.-1) then + fp(ipat)=npindex + else + do while (np(ip).ne.-1) + ip=np(ip) + enddo + np(ip)=npindex + endif + return +end subroutine boxit + +subroutine fetchit(reset,e2,ntau,i1,i2) + integer indexes(4000,2),fp(0:525000),np(4000) + integer lastpat + integer*1 e2(ntau) + logical reset + common/boxes/indexes,fp,np + save lastpat,inext + + if(reset) then + lastpat=-1 + reset=.false. + endif + + ipat=0 + do i=1,ntau + if(e2(i).eq.1) then + ipat=ipat+ishft(1,ntau-i) + endif + enddo + index=fp(ipat) + + if(lastpat.ne.ipat .and. index.gt.0) then ! return first set of indices + i1=indexes(index,1) + i2=indexes(index,2) + inext=np(index) + elseif(lastpat.eq.ipat .and. inext.gt.0) then + i1=indexes(inext,1) + i2=indexes(inext,2) + inext=np(inext) + else + i1=-1 + i2=-1 + inext=-1 + endif + lastpat=ipat + return +end subroutine fetchit + diff --git a/lib/fsk4hf/osdwspr.f90 b/lib/fsk4hf/osdwspr.f90 new file mode 100644 index 000000000..70a89c17c --- /dev/null +++ b/lib/fsk4hf/osdwspr.f90 @@ -0,0 +1,374 @@ +subroutine osdwspr(llr,apmask,ndeep,decoded,cw,nhardmin,dmin) +! +use iso_c_binding +parameter (N=162, K=50, L=32) + +!integer*1 p1(L),p2(L),p3(L),p4(L) +integer*1 gg(64) + +integer*1 apmask(N),apmaskr(N) +integer*1 gen(K,N) +integer*1 genmrb(K,N),g2(N,K) +integer*1 temp(K),m0(K),me(K),mi(K),misub(K),e2sub(N-K),e2(N-K),ui(N-K) +integer*1 r2pat(N-K) +integer indices(N),nxor(N) +integer*1 cw(N),ce(N),c0(N),hdec(N) +integer*1 decoded(K) +integer indx(N) +real llr(N),rx(N),absrx(N) +logical first,reset +data first/.true./ +!data p1/1,0,0,0,1,1,0,0,0,0,1,0,0,0,1,0,1,1,0,0,1,1,1,0,1/ +!data p2/1,0,1,0,1,0,1,1,0,0,1,0,0,1,0,1,0,0,1,0,0,1,1,1,1/ +!data p3/1,1,0,1,0,1,0,1,1,0,1,0,0,1,1,0,1,1,1,1,1,1,0,1,1/ +!data p4/1,1,1,0,1,1,1,1,0,1,1,1,0,0,0,1,1,1,1,1,1,0,0,0,1/ +!data p1/1,0,1,0,1,1,0,0,1,1,0,1,1,1,1,1/ +!data p2/1,0,1,1,0,1,0,0,1,1,1,1,1,0,0,1/ +!data p3/1,1,0,0,1,0,1,1,0,1,1,1,0,0,1,1/ +!data p4/1,1,1,0,1,1,0,1,1,1,1,0,0,1,0,1/ +data gg/1,1,0,1,0,1,0,0,1,0,0,0,1,1,0,0,1,0,1,0,0,1,0,1,1,1,0,1,1,0,0,0, & + 0,1,0,0,0,0,0,0,1,0,0,1,1,1,1,0,0,0,1,0,0,1,0,0,1,0,1,1,1,1,1,1/ + +save first,gen + +if( first ) then ! fill the generator matrix +! gg=0 +! gg(1:L)=p1 +! gg(L+1:2*L)=p2 +! gg(2*L+1:3*L)=p3 +! gg(3*L+1:4*L)=p4 + gen=0 + gen(1,1:2*L)=gg(1:2*L) + do i=2,K + gen(i,:)=cshift(gen(i-1,:),-2) + enddo + first=.false. +endif + +! Re-order received vector to place systematic msg bits at the end. +rx=llr +apmaskr=apmask + +! Hard decisions on the received word. +hdec=0 +where(rx .ge. 0) hdec=1 + +! Use magnitude of received symbols as a measure of reliability. +absrx=abs(rx) +call indexx(absrx,N,indx) + +! Re-order the columns of the generator matrix in order of decreasing reliability. +do i=1,N + genmrb(1:K,i)=gen(1:K,indx(N+1-i)) + indices(i)=indx(N+1-i) +enddo + +! Do gaussian elimination to create a generator matrix with the most reliable +! received bits in positions 1:K in order of decreasing reliability (more or less). +do id=1,K ! diagonal element indices + do icol=id,K+20 ! The 20 is ad hoc - beware + iflag=0 + if( genmrb(id,icol) .eq. 1 ) then + iflag=1 + if( icol .ne. id ) then ! reorder column + temp(1:K)=genmrb(1:K,id) + genmrb(1:K,id)=genmrb(1:K,icol) + genmrb(1:K,icol)=temp(1:K) + itmp=indices(id) + indices(id)=indices(icol) + indices(icol)=itmp + endif + do ii=1,K + if( ii .ne. id .and. genmrb(ii,id) .eq. 1 ) then + genmrb(ii,1:N)=ieor(genmrb(ii,1:N),genmrb(id,1:N)) + endif + enddo + exit + endif + enddo +enddo + +g2=transpose(genmrb) + +! The hard decisions for the K MRB bits define the order 0 message, m0. +! Encode m0 using the modified generator matrix to find the "order 0" codeword. +! Flip various combinations of bits in m0 and re-encode to generate a list of +! codewords. Return the member of the list that has the smallest Euclidean +! distance to the received word. + +hdec=hdec(indices) ! hard decisions from received symbols +m0=hdec(1:K) ! zero'th order message +absrx=absrx(indices) +rx=rx(indices) +apmaskr=apmaskr(indices) + +call mrbencode(m0,c0,g2,N,K) +nxor=ieor(c0,hdec) +nhardmin=sum(nxor) +dmin=sum(nxor*absrx) + +cw=c0 +ntotal=0 +nrejected=0 + +if(ndeep.eq.0) goto 998 ! norder=0 +if(ndeep.gt.5) ndeep=5 +if( ndeep.eq. 1) then + nord=1 + npre1=0 + npre2=0 + nt=120 + ntheta=12 +elseif(ndeep.eq.2) then + nord=1 + npre1=1 + npre2=0 + nt=120 + ntheta=12 +elseif(ndeep.eq.3) then + nord=1 + npre1=1 + npre2=1 + nt=120 + ntheta=12 + ntau=15 +elseif(ndeep.eq.4) then + nord=2 + npre1=1 + npre2=0 + nt=120 + ntheta=12 + ntau=15 +elseif(ndeep.eq.5) then + nord=3 + npre1=1 + npre2=1 + nt=80 + ntheta=22 + ntau=16 +endif + +do iorder=1,nord + misub(1:K-iorder)=0 + misub(K-iorder+1:K)=1 + iflag=K-iorder+1 + do while(iflag .ge.0) + if(iorder.eq.nord .and. npre1.eq.0) then + iend=iflag + else + iend=1 + endif + do n1=iflag,iend,-1 + mi=misub + mi(n1)=1 + if(any(iand(apmaskr(1:K),mi).eq.1)) cycle + ntotal=ntotal+1 + me=ieor(m0,mi) + if(n1.eq.iflag) then + call mrbencode(me,ce,g2,N,K) + e2sub=ieor(ce(K+1:N),hdec(K+1:N)) + e2=e2sub + nd1Kpt=sum(e2sub(1:nt))+1 + d1=sum(ieor(me(1:K),hdec(1:K))*absrx(1:K)) + else + e2=ieor(e2sub,g2(K+1:N,n1)) + nd1Kpt=sum(e2(1:nt))+2 + endif + if(nd1Kpt .le. ntheta) then + call mrbencode(me,ce,g2,N,K) + nxor=ieor(ce,hdec) + if(n1.eq.iflag) then + dd=d1+sum(e2sub*absrx(K+1:N)) + else + dd=d1+ieor(ce(n1),hdec(n1))*absrx(n1)+sum(e2*absrx(K+1:N)) + endif + if( dd .lt. dmin ) then + dmin=dd + cw=ce + nhardmin=sum(nxor) + nd1Kptbest=nd1Kpt + endif + else + nrejected=nrejected+1 + endif + enddo +! Get the next test error pattern, iflag will go negative +! when the last pattern with weight iorder has been generated. + call nextpat(misub,k,iorder,iflag) + enddo +enddo + +if(npre2.eq.1) then + reset=.true. + ntotal=0 + do i1=K,1,-1 + do i2=i1-1,1,-1 + ntotal=ntotal+1 + mi(1:ntau)=ieor(g2(K+1:K+ntau,i1),g2(K+1:K+ntau,i2)) + call boxit(reset,mi(1:ntau),ntau,ntotal,i1,i2) + enddo + enddo + + ncount2=0 + ntotal2=0 + reset=.true. +! Now run through again and do the second pre-processing rule + misub(1:K-nord)=0 + misub(K-nord+1:K)=1 + iflag=K-nord+1 + do while(iflag .ge.0) + me=ieor(m0,misub) + call mrbencode(me,ce,g2,N,K) + e2sub=ieor(ce(K+1:N),hdec(K+1:N)) + do i2=0,ntau + ntotal2=ntotal2+1 + ui=0 + if(i2.gt.0) ui(i2)=1 + r2pat=ieor(e2sub,ui) +778 continue + call fetchit(reset,r2pat(1:ntau),ntau,in1,in2) + if(in1.gt.0.and.in2.gt.0) then + ncount2=ncount2+1 + mi=misub + mi(in1)=1 + mi(in2)=1 + if(sum(mi).lt.nord+npre1+npre2.or.any(iand(apmaskr(1:K),mi).eq.1)) cycle + me=ieor(m0,mi) + call mrbencode(me,ce,g2,N,K) + nxor=ieor(ce,hdec) + dd=sum(nxor*absrx) + if( dd .lt. dmin ) then + dmin=dd + cw=ce + nhardmin=sum(nxor) + endif + goto 778 + endif + enddo + call nextpat(misub,K,nord,iflag) + enddo +endif + +998 continue +! Re-order the codeword to place message bits at the end. +cw(indices)=cw +hdec(indices)=hdec +decoded=0 +return +end subroutine osdwspr + +subroutine mrbencode(me,codeword,g2,N,K) +integer*1 me(K),codeword(N),g2(N,K) +! fast encoding for low-weight test patterns + codeword=0 + do i=1,K + if( me(i) .eq. 1 ) then + codeword=ieor(codeword,g2(1:N,i)) + endif + enddo +return +end subroutine mrbencode + +subroutine nextpat(mi,k,iorder,iflag) + integer*1 mi(k),ms(k) +! generate the next test error pattern + ind=-1 + do i=1,k-1 + if( mi(i).eq.0 .and. mi(i+1).eq.1) ind=i + enddo + if( ind .lt. 0 ) then ! no more patterns of this order + iflag=ind + return + endif + ms=0 + ms(1:ind-1)=mi(1:ind-1) + ms(ind)=1 + ms(ind+1)=0 + if( ind+1 .lt. k ) then + nz=iorder-sum(ms) + ms(k-nz+1:k)=1 + endif + mi=ms + do i=1,k ! iflag will point to the lowest-index 1 in mi + if(mi(i).eq.1) then + iflag=i + exit + endif + enddo + return +end subroutine nextpat + +subroutine boxit(reset,e2,ntau,npindex,i1,i2) + integer*1 e2(1:ntau) + integer indexes(4000,2),fp(0:525000),np(4000) + logical reset + common/boxes/indexes,fp,np + + if(reset) then + patterns=-1 + fp=-1 + np=-1 + sc=-1 + indexes=-1 + reset=.false. + endif + + indexes(npindex,1)=i1 + indexes(npindex,2)=i2 + ipat=0 + do i=1,ntau + if(e2(i).eq.1) then + ipat=ipat+ishft(1,ntau-i) + endif + enddo + + ip=fp(ipat) ! see what's currently stored in fp(ipat) + if(ip.eq.-1) then + fp(ipat)=npindex + else + do while (np(ip).ne.-1) + ip=np(ip) + enddo + np(ip)=npindex + endif + return +end subroutine boxit + +subroutine fetchit(reset,e2,ntau,i1,i2) + integer indexes(4000,2),fp(0:525000),np(4000) + integer lastpat + integer*1 e2(ntau) + logical reset + common/boxes/indexes,fp,np + save lastpat,inext + + if(reset) then + lastpat=-1 + reset=.false. + endif + + ipat=0 + do i=1,ntau + if(e2(i).eq.1) then + ipat=ipat+ishft(1,ntau-i) + endif + enddo + index=fp(ipat) + + if(lastpat.ne.ipat .and. index.gt.0) then ! return first set of indices + i1=indexes(index,1) + i2=indexes(index,2) + inext=np(index) + elseif(lastpat.eq.ipat .and. inext.gt.0) then + i1=indexes(inext,1) + i2=indexes(inext,2) + inext=np(inext) + else + i1=-1 + i2=-1 + inext=-1 + endif + lastpat=ipat + return +end subroutine fetchit + diff --git a/lib/fsk4hf/tccsim.f90 b/lib/fsk4hf/tccsim.f90 new file mode 100644 index 000000000..1cbe8d172 --- /dev/null +++ b/lib/fsk4hf/tccsim.f90 @@ -0,0 +1,194 @@ +! +! Simulator for terminated convolutional codes (so, far, only rate 1/2) +! BPSK on AWGN Channel +! +! Hybrid decoder - Fano Sequential Decoder and Ordered Statistics Decoder (OSD)a +! +program tccsim + + parameter (N=162,K=50) + integer*1 gen(K,N) + integer*1 gg(64) + integer*1 mbits(50),mbits2(50) + integer*4 mettab(-128:127,0:1) + + parameter (NSYM=162) + parameter (MAXSYM=220) + character*12 arg + character*22 msg,msg2 + integer*1 data0(11) + integer*1 data1(11) + integer*1 dat(162) + integer*1 softsym(162) + integer*1 apmask(162),cw(162) + real*4 xx0(0:255) + real ss(162) + + character*64 g ! Interleaved polynomial coefficients + data g/'1101010010001100101001011101100001000000100111100010010010111111'/ + + data xx0/ & !Metric table + 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, & + 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, & + 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, & + 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, & + 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, & + 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, & + 0.988, 1.000, 0.991, 0.993, 1.000, 0.995, 1.000, 0.991, & + 1.000, 0.991, 0.992, 0.991, 0.990, 0.990, 0.992, 0.996, & + 0.990, 0.994, 0.993, 0.991, 0.992, 0.989, 0.991, 0.987, & + 0.985, 0.989, 0.984, 0.983, 0.979, 0.977, 0.971, 0.975, & + 0.974, 0.970, 0.970, 0.970, 0.967, 0.962, 0.960, 0.957, & + 0.956, 0.953, 0.942, 0.946, 0.937, 0.933, 0.929, 0.920, & + 0.917, 0.911, 0.903, 0.895, 0.884, 0.877, 0.869, 0.858, & + 0.846, 0.834, 0.821, 0.806, 0.790, 0.775, 0.755, 0.737, & + 0.713, 0.691, 0.667, 0.640, 0.612, 0.581, 0.548, 0.510, & + 0.472, 0.425, 0.378, 0.328, 0.274, 0.212, 0.146, 0.075, & + 0.000,-0.079,-0.163,-0.249,-0.338,-0.425,-0.514,-0.606, & + -0.706,-0.796,-0.895,-0.987,-1.084,-1.181,-1.280,-1.376, & + -1.473,-1.587,-1.678,-1.790,-1.882,-1.992,-2.096,-2.201, & + -2.301,-2.411,-2.531,-2.608,-2.690,-2.829,-2.939,-3.058, & + -3.164,-3.212,-3.377,-3.463,-3.550,-3.768,-3.677,-3.975, & + -4.062,-4.098,-4.186,-4.261,-4.472,-4.621,-4.623,-4.608, & + -4.822,-4.870,-4.652,-4.954,-5.108,-5.377,-5.544,-5.995, & + -5.632,-5.826,-6.304,-6.002,-6.559,-6.369,-6.658,-7.016, & + -6.184,-7.332,-6.534,-6.152,-6.113,-6.288,-6.426,-6.313, & + -9.966,-6.371,-9.966,-7.055,-9.966,-6.629,-6.313,-9.966, & + -5.858,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966, & + -9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966, & + -9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966, & + -9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966, & + -9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966, & + -9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966,-9.966/ + + bias=0.42 + scale=120 +! ndelta=nint(3.4*scale) + ndelta=100 + ib=150 + slope=2 + do i=0,255 + mettab(i-128,0)=nint(scale*(xx0(i)-bias)) + if(i.gt.ib) mettab(i-128,0)=mettab(ib-128,0)-slope*(i-ib) + if(i.ge.1) mettab(128-i,1)=mettab(i-128,0) + enddo + mettab(-128,1)=mettab(-127,1) + +! Get command-line argument(s) + nargs=iargc() + if(nargs.ne.3) then + print*,'Usage: tccsim "message" ntrials ndepth' + go to 999 + endif + call getarg(1,msg) !Get message from command line + write(*,1000) msg +1000 format('Message: ',a22) + call getarg(2,arg) + read(arg,*) ntrials + call getarg(3,arg) + read(arg,*) ndepth + + nbits=50+31 !User bits=99, constraint length=32 + nbytes=(nbits+7)/8 + limit=20000 + + data0=0 +! call wqencode(msg,ntype0,data0) !Source encoding +data0(1:7)=1 + write(*,1002) data0(1:7),data0(1:6),data0(7)/64 +1002 format(/'Source-encoded message, 50 bits:'/'Hex: ',7z3.2/ & + 'Binary: ',6b9.8,b3.2) + +! Demonstrates how to create the generator matrix from a string that contains the interleaved +! polynomial coefficients + gen=0 + do j=1,64 + read(g(j:j),"(i1)") gg(j) ! read polynomial coeffs from string + enddo + do i=1,K + gen(i,2*(i-1)+1:2*(i-1)+64)=gg ! fill the generator matrix with cyclic shifts of gg + enddo + +! get message bits from data0 + nbits=0 + do i=1,7 + do ib=7,0,-1 + nbits=nbits+1 + if(nbits .le. 50) then + mbits(nbits)=0 + if(btest(data0(i),ib)) mbits(nbits)=1 + endif + enddo + enddo + +! Encode message bits using the generator matrix, generating a 162-bit codeword. + cw=0 + do i=1,50 + if(mbits(i).eq.1) cw=mod(cw+gen(i,:),2) + enddo + + write(*,*) 'Codeword from generator matrix: ' + write(*,'(162i1)') cw + + call encode232(data0,nbytes,dat,MAXSYM) !Convolutional encoding + write(*,*) 'Codeword from encode232: ' + write(*,'(162i1)') dat + +! call inter_mept(dat,1) !Interleaving + +! Here, we have channel symbols. + +! call inter_mept(dat,-1) !Remove interleaving + + call init_random_seed() + call sgran() + + do isnr=10,-20,-1 + sigma=1/sqrt(2*(10**((isnr/2.0)/10.0))) + ngood=0 + nbad=0 + do i=1,ntrials + do j=1,162 + ss(j)=-(2*dat(j)-1)+sigma*gran() !Simulate soft symbols + enddo + + rms=sqrt(sum(ss**2)) + ss=100*ss/rms + where(ss>127.0) ss=127.0 + where(ss<-127.0) ss=-127.0 + softsym=ss + +! Call the sequential (Fano algorithm) decoder + nbits=50+31 + call fano232(softsym,nbits,mettab,ndelta,limit,data1,ncycles,metric,nerr) + iflag=0 + nhardmin=0 + dmin=0.0 + +! If Fano fails, call OSD + if(nerr.ne.0 .and. ndepth.ge.0) then + apmask=0 + cw=0 + call osdwspr(-ss/127,apmask,ndepth,cw,nhardmin,dmin) + +! OSD produces a codeword, but code is not systematic +! Use Fano with hard decisions to retrieve the message from the codeword + cw=-(2*cw-1)*64 + nbits=50+31 + call fano232(cw,nbits,mettab,ndelta,limit,data1,ncycles,metric,nerr) + iflag=1 + endif + +! call wqdecode(data1,msg2,ntype1) +! write(*,*) msg2,iflag,nhardmin,dmin + if(nerr.eq.0 .and. any(data1.ne.data0)) nbad=nbad+1 + if(nerr.eq.0 .and. all(data1.eq.data0)) ngood=ngood+1 + enddo + + write(*,'(f4.1,i8,i8)') isnr/2.0,ngood,nbad + enddo + +999 end program tccsim + +include '../wsprcode/wspr_old_subs.f90' +