diff --git a/CMakeLists.txt b/CMakeLists.txt index 2eefc415e..71fd5f789 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/wsprd/osdwspr.f90 lib/pctile.f90 lib/peakdt9.f90 lib/peakup.f90 @@ -586,6 +587,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 ) @@ -1195,7 +1198,7 @@ add_executable (wsprcode lib/wsprcode/wsprcode.f90 lib/wsprcode/nhash.c wsjtx.rc) target_link_libraries (wsprcode wsjt_fort wsjt_cxx) -add_executable (wsprd ${wsprd_CSRCS}) +add_executable (wsprd ${wsprd_CSRCS} lib/indexx.f90 lib/wsprd/osdwspr.f90) target_include_directories (wsprd PRIVATE ${FFTW3_INCLUDE_DIRS}) target_link_libraries (wsprd ${FFTW3_LIBRARIES}) @@ -1237,6 +1240,9 @@ target_link_libraries (ft8sim 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..26069eb14 --- /dev/null +++ b/lib/fsk4hf/osdwspr.f90 @@ -0,0 +1,373 @@ +subroutine osdwspr(ss,apmask,ndeep,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) + +real ss(N) +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 indx(N),ndeep,nhardmin +real rx(N),absrx(N),dmin +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=ss/127.0 +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=60 + ntheta=12 +elseif(ndeep.eq.2) then + nord=1 + npre1=1 + npre2=0 + nt=60 + ntheta=12 +elseif(ndeep.eq.3) then + nord=1 + npre1=1 + npre2=1 + nt=60 + ntheta=22 + ntau=16 +elseif(ndeep.eq.4) then + nord=2 + npre1=1 + npre2=0 + nt=60 + ntheta=22 + ntau=16 +elseif(ndeep.eq.5) then + nord=3 + npre1=1 + npre2=1 + nt=60 + 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 as-received order. +cw(indices)=cw +hdec(indices)=hdec +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..644b87a39 --- /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=162) + character*12 arg + character*22 msg,msg2 + integer*1 data0(13) + integer*1 data1(13) + integer*1 dat(206) + integer*1 softsym(162) + integer*1 apmask(162),cw0(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 +write(*,*) 'data0 ',data0 +! 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 + + write(*,*) 'Source encoded message bits: ' + write(*,'(6(8i1,1x),2i1)') mbits + +! Encode message bits using the generator matrix, generating a 162-bit codeword. + cw0=0 + do i=1,50 + if(mbits(i).eq.1) cw0=mod(cw0+gen(i,:),2) + enddo + + write(*,*) 'Codeword from generator matrix: ' + write(*,'(162i1)') cw0 + +! call encode232(data0,nbytes,dat) !Convolutional encoding +! write(*,*) 'Codeword from encode232: ' +! write(*,'(162i2)') 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*cw0(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) +nerr=1 + 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(softsym,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 +!dat=0 +!dat(1:162)=cw +! call fano232(dat,nbits,mettab,ndelta,limit,data1,ncycles,metric,nerr) +! iflag=1 + endif +! call wqdecode(data1,msg2,ntype1) +! write(*,*) msg2,iflag,nhardmin,dmin + if(any(cw.ne.cw0)) nbad=nbad+1 + if(all(cw.eq.cw0)) 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' + diff --git a/lib/wsprd/osdwspr.f90 b/lib/wsprd/osdwspr.f90 new file mode 100644 index 000000000..42b8e949c --- /dev/null +++ b/lib/wsprd/osdwspr.f90 @@ -0,0 +1,359 @@ +subroutine osdwspr(ss,apmask,ndeep,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) + +real ss(N) +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 indx(N),ndeep,nhardmin +real rx(N),absrx(N),dmin +logical first,reset +data first/.true./ +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 + 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 + +rx=ss/127.0 +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.le.0) goto 998 ! norder=0 +if(ndeep.gt.5) ndeep=5 +if( ndeep.eq. 1) then + nord=1 + npre1=0 + npre2=0 + nt=66 + ntheta=16 +elseif(ndeep.eq.2) then + nord=1 + npre1=1 + npre2=0 + nt=66 + ntheta=16 +elseif(ndeep.eq.3) then + nord=2 + npre1=1 + npre2=0 + nt=66 + ntheta=22 + ntau=16 +elseif(ndeep.eq.4) then + nord=2 + npre1=1 + npre2=1 + nt=66 + ntheta=22 + ntau=16 +elseif(ndeep.eq.5) then + nord=3 + npre1=1 + npre2=0 + nt=66 + ntheta=22 + ntau=20 +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 as-received order. +cw(indices)=cw +hdec(indices)=hdec +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/wsprd/wsprd.c b/lib/wsprd/wsprd.c index 48ec1faeb..024afaa48 100644 --- a/lib/wsprd/wsprd.c +++ b/lib/wsprd/wsprd.c @@ -693,6 +693,7 @@ void usage(void) printf(" -H do not use (or update) the hash table\n"); printf(" -J use the stack decoder instead of Fano decoder\n"); printf(" -m decode wspr-15 .wav file\n"); + printf(" -o n (0<=n<=5), decoding depth for OSD, default is disabled\n"); printf(" -q quick mode - doesn't dig deep for weak signals\n"); printf(" -s single pass mode, no subtraction (same as original wsprd)\n"); printf(" -v verbose mode (shows dupes)\n"); @@ -708,7 +709,7 @@ int main(int argc, char *argv[]) extern char *optarg; extern int optind; int i,j,k; - unsigned char *symbols, *decdata, *channel_symbols; + unsigned char *symbols, *decdata, *channel_symbols, *apmask, *cw; signed char message[]={-9,13,-35,123,57,-39,64,0,0,0,0}; char *callsign, *call_loc_pow; char *ptr_to_infile,*ptr_to_infile_suffix; @@ -717,19 +718,22 @@ int main(int argc, char *argv[]) char timer_fname[200],hash_fname[200]; char uttime[5],date[7]; int c,delta,maxpts=65536,verbose=0,quickmode=0,more_candidates=0, stackdecoder=0; - int writenoise=0,usehashtable=1,wspr_type=2, ipass, nblocksize; + int writenoise=0,usehashtable=1,wspr_type=2, ipass, nblocksize,use_osd=0; + int nhardmin,ihash; int writec2=0,maxdrift; int shift1, lagmin, lagmax, lagstep, ifmin, ifmax, worth_a_try, not_decoded; unsigned int nbits=81, stacksize=200000; unsigned int npoints, metric, cycles, maxnp; float df=375.0/256.0/2; float freq0[200],snr0[200],drift0[200],sync0[200]; + float fsymbs[162]; int shift0[200]; float dt=1.0/375.0, dt_print; double dialfreq_cmdline=0.0, dialfreq, freq_print; double dialfreq_error=0.0; float fmin=-110, fmax=110; float f1, fstep, sync1, drift1; + float dmin; float psavg[512]; float *idat, *qdat; clock_t t0,t00; @@ -738,13 +742,16 @@ int main(int argc, char *argv[]) struct result { char date[7]; char time[5]; float sync; float snr; float dt; double freq; char message[23]; float drift; - unsigned int cycles; int jitter; int blocksize; unsigned int metric; }; + unsigned int cycles; int jitter; int blocksize; unsigned int metric; + unsigned char osd_decode }; struct result decodes[50]; char *hashtab; hashtab=calloc(32768*13,sizeof(char)); int nh; symbols=calloc(nbits*2,sizeof(char)); + apmask=calloc(162,sizeof(char)); + cw=calloc(162,sizeof(char)); decdata=calloc(11,sizeof(char)); channel_symbols=calloc(nbits*2,sizeof(char)); callsign=calloc(13,sizeof(char)); @@ -765,6 +772,7 @@ int main(int argc, char *argv[]) int block_demod=1; //Default is to use block demod on pass 2 int subtraction=1; int npasses=2; + int ndepth=-1; //Depth for OSD float minrms=52.0 * (symfac/64.0); //Final test for plausible decoding delta=60; //Fano threshold step @@ -779,7 +787,7 @@ int main(int argc, char *argv[]) idat=calloc(maxpts,sizeof(float)); qdat=calloc(maxpts,sizeof(float)); - while ( (c = getopt(argc, argv, "a:BcC:de:f:HJmqstwvz:")) !=-1 ) { + while ( (c = getopt(argc, argv, "a:BcC:de:f:HJmo:qstwvz:")) !=-1 ) { switch (c) { case 'a': data_dir = optarg; @@ -812,6 +820,9 @@ int main(int argc, char *argv[]) case 'm': //15-minute wspr mode wspr_type = 15; break; + case 'o': //use ordered-statistics-decoder + ndepth=(int) strtol(optarg,NULL,10); + break; case 'q': //no shift jittering quickmode = 1; break; @@ -1242,7 +1253,9 @@ int main(int argc, char *argv[]) int idt, ii, jittered_shift; float y,sq,rms; not_decoded=1; - int ib=1, blocksize; + int osd_decode=0; + int ib=1, blocksize; + int n1,n2; while( ib <= nblocksize && not_decoded ) { blocksize=ib; idt=0; ii=0; @@ -1278,6 +1291,34 @@ int main(int argc, char *argv[]) } tfano += (float)(clock()-t0)/CLOCKS_PER_SEC; + + if( (ndepth >= 0) && not_decoded ) { + for(i=0; i<162; i++) { + fsymbs[i]=symbols[i]-127; + } + osdwspr_(fsymbs,apmask,&ndepth,cw,&nhardmin,&dmin); + for(i=0; i<162; i++) { + symbols[i]=255*cw[i]; + } + fano(&metric,&cycles,&maxnp,decdata,symbols,nbits, + mettab,delta,maxcycles); + for(i=0; i<11; i++) { + if( decdata[i]>127 ) { + message[i]=decdata[i]-256; + } else { + message[i]=decdata[i]; + } + } + unpack50(message,&n1,&n2); + if( !unpackcall(n1,callsign) ) break; + callsign[12]=0; + ihash=nhash(callsign,strlen(callsign),(uint32_t)146); + if(strncmp(hashtab+ihash*13,callsign,13)==0) { + not_decoded=0; + osd_decode =1; + break; + } + } } idt++; @@ -1285,7 +1326,7 @@ int main(int argc, char *argv[]) } ib++; } - + if( worth_a_try && !not_decoded ) { ndecodes_pass++; @@ -1347,6 +1388,7 @@ int main(int argc, char *argv[]) decodes[uniques-1].jitter=ii; decodes[uniques-1].blocksize=blocksize; decodes[uniques-1].metric=metric; + decodes[uniques-1].osd_decode=osd_decode; } } } @@ -1378,11 +1420,11 @@ int main(int argc, char *argv[]) decodes[i].time, decodes[i].snr,decodes[i].dt, decodes[i].freq, (int)decodes[i].drift, decodes[i].message); fprintf(fall_wspr, - "%6s %4s %3d %3.0f %5.2f %11.7f %-22s %2d %5u %4d %4d %4d\n", + "%6s %4s %3d %3.0f %5.2f %11.7f %-22s %2d %5u %4d %4d %4d %2u\n", decodes[i].date, decodes[i].time, (int)(10*decodes[i].sync), decodes[i].snr, decodes[i].dt, decodes[i].freq, decodes[i].message, (int)decodes[i].drift, decodes[i].cycles/81, - decodes[i].jitter,decodes[i].blocksize,decodes[i].metric); + decodes[i].jitter,decodes[i].blocksize,decodes[i].metric,decodes[i].osd_decode); fprintf(fwsprd, "%6s %4s %3d %3.0f %4.1f %10.6f %-22s %2d %5u %4d\n", decodes[i].date, decodes[i].time, (int)(10*decodes[i].sync), diff --git a/mainwindow.cpp b/mainwindow.cpp index 37f120623..7fd49db9a 100644 --- a/mainwindow.cpp +++ b/mainwindow.cpp @@ -1379,7 +1379,7 @@ void MainWindow::dataSink(qint64 frames) t2.sprintf(" -f %.6f ",f0m1500); if((m_ndepth&7)==1) depth_string=" -qB "; //2 pass w subtract, no Block detection, no shift jittering if((m_ndepth&7)==2) depth_string=" -B "; //2 pass w subtract, no Block detection - if((m_ndepth&7)==3) depth_string=" "; //2 pass w subtract, Block detection + if((m_ndepth&7)==3) depth_string=" -C 5000 -o 4"; //2 pass w subtract, Block detection and OSD. QString degrade; degrade.sprintf("-d %4.1f ",m_config.degrade());