diff --git a/lib/detectmsk32.f90 b/lib/detectmsk32.f90 index 7041a25c4..076bbe49b 100644 --- a/lib/detectmsk32.f90 +++ b/lib/detectmsk32.f90 @@ -1,7 +1,7 @@ subroutine detectmsk32(cbig,n,mycall,partnercall,lines,nmessages,nutc,ntol,t00) use timer_module, only: timer - parameter (NSPM=192, NPTS=3*NSPM, MAXSTEPS=7500, NFFT=3*NSPM, MAXCAND=10) + parameter (NSPM=192, NPTS=3*NSPM, MAXSTEPS=7500, NFFT=3*NSPM, MAXCAND=5) character*4 rpt(0:63) character*6 mycall,partnercall character*22 msg,msgsent,msgreceived,allmessages(32) @@ -19,6 +19,8 @@ subroutine detectmsk32(cbig,n,mycall,partnercall,lines,nmessages,nutc,ntol,t00) real detmet(-2:MAXSTEPS+3) real detfer(MAXSTEPS) real ferrs(MAXCAND) + real f2(64) + real peak(64) real pp(12) real rcw(12) real snrs(MAXCAND) @@ -166,13 +168,13 @@ subroutine detectmsk32(cbig,n,mycall,partnercall,lines,nmessages,nutc,ntol,t00) nmessages=0 lines=char(0) - fbest=1e6 - pkbest=-1e6 + pkbest=-1.0 + ratbest=0.0 imsgbest=-1 - istartbest=-1 + fbest=0.0 ipbest=-1 - nsnrbest=-100 - t0best=-1e6 + nsnrbest=-1 + do ip=1,ndet !run through the candidates and try to sync/demod/decode imid=times(ip)*fs if( imid .lt. NPTS/2 ) imid=NPTS/2 @@ -180,12 +182,30 @@ subroutine detectmsk32(cbig,n,mycall,partnercall,lines,nmessages,nutc,ntol,t00) t0=times(ip) + t00 cdat=cbig(imid-NPTS/2+1:imid+NPTS/2) ferr=ferrs(ip) - nsnr=2*nint(snrs(ip)/2.0) - if( nsnr .lt. -4 ) nsnr=-4 - if( nsnr .gt. 24 ) nsnr=24 + nsnr=nint(snrs(ip)) + if( nsnr .lt. -5 ) nsnr=-5 + if( nsnr .gt. 25 ) nsnr=25 - do imsg=1,64 - do istart=NSPM-NSPM/2,NPTS-NSPM + ic0=NSPM + do i=1,6 +! if( ic0+11+NSPM .le. NPTS ) then + bb(i) = sum( ( cdat(ic0+i-1+6:ic0+i-1+6+NSPM:6) * conjg( cdat(ic0+i-1:ic0+i-1+NSPM:6) ) )**2 ) +! else +! bb(i) = sum( ( cdat(ic0+i-1+6:NPTS:6) * conjg( cdat(ic0+i-1:NPTS-6:6) ) )**2 ) +! endif +! write(*,*) ip,i,abs(bb(i)) + enddo + iloc=maxloc(abs(bb)) + ibb=iloc(1) + bba=abs(bb(ibb)) + bbp=atan2(-imag(bb(ibb)),-real(bb(ibb)))/(2*twopi*6*dt) + if( ibb .le. 3 ) ibb=ibb-1 + if( ibb .gt. 3 ) ibb=ibb-7 +! write(*,*) ibb + ic0=ic0+ibb+2 + + do istart=ic0,ic0+32*6-1,6 + do imsg=1,64 cft(1:144)=cdat(istart:istart+144-1)*conjg(cwaveforms(1:144,imsg)) cft(145:512)=0. df=12000.0/512.0 @@ -195,26 +215,36 @@ subroutine detectmsk32(cbig,n,mycall,partnercall,lines,nmessages,nutc,ntol,t00) pk=abs(cft(ipk)) fpk=(ipk-1)*df if( fpk.gt.12000.0 ) fpk=fpk-12000.0 - if( pk .gt. pkbest .and. abs(fpk-1500.0) .le. ntol) then - ipbest=ip - pkbest=pk - fbest=fpk - imsgbest=imsg - istartbest=istart - nsnrbest=nsnr - t0best=t0 - endif + f2(imsg)=fpk + peak(imsg)=pk enddo + iloc=maxloc(peak) + imsg1=iloc(1) + pk1=peak(imsg1) + peak(imsg1)=-1 + pk2=maxval(peak) + rat=pk1/pk2 + if( abs(f2(imsg1)-1500) .le. ntol .and. (pk1 .gt. pkbest) ) then + pkbest=pk1 + ratbest=rat + imsgbest=imsg1 + fbest=f2(imsg1) + ipbest=ip + t0best=t0 + nsnrbest=nsnr + istartbest=istart + endif enddo +! write(*,*) ip,imid,istart,imsgbest,pkbest,ratbest,nsnrbest +! if( pkbest .gt. 110.0 .and. ratbest .gt. 1.2 ) goto 999 enddo ! candidate loop 999 continue msgreceived=' ' - if( imsgbest .gt. 0 .and. pkbest .ge. 108.0) then + if( imsgbest .gt. 0 .and. pkbest .ge. 110.0 .and. ratbest .ge. 1.20) then nrxrpt=iand(imsgbest-1,63) - nrxhash=(imsgbest-1-nrxrpt)/64 -!write(*,*) ipbest,pkbest,fbest,imsgbest,istartbest,nsnrbest,t0best,nrxrpt,nrxhash +!write(*,*) ipbest,pkbest,fbest,imsgbest,istartbest,nsnrbest,t0best nmessages=1 write(msgreceived,'(a1,a,1x,a,a1,1x,a4)') "<",trim(mycall), & trim(partnercall),">",rpt(nrxrpt) diff --git a/lib/msk32d.f90 b/lib/msk32d.f90 index 24238ae09..ff2138308 100644 --- a/lib/msk32d.f90 +++ b/lib/msk32d.f90 @@ -4,12 +4,12 @@ program msk32d parameter (NSPM=6*32) complex c0(0:NZ0-1) complex c(0:NZ-1) - complex cmsg(0:NSPM-1,0:31) + complex cmsg(0:NSPM-1,0:63) complex z real a(3) real p0(0:NSPM-1) real p(0:NSPM-1) - real s0(0:31) + real s0(0:63) real dd(NZ) integer itone(144) integer ihdr(11) @@ -17,12 +17,7 @@ program msk32d integer*2 id2(NZ) character*22 msg,msgsent character mycall*8,hiscall*6,arg*12,infile*80,datetime*13 - character*4 rpt(0:31) - data rpt /'-04 ','-02 ','+00 ','+02 ','+04 ','+06 ','+08 ','+10 ','+12 ', & - '+14 ','+16 ','+18 ','+20 ','+22 ','+24 ', & - 'R-04','R-02','R+00','R+02','R+04','R+06','R+08','R+10','R+12', & - 'R+14','R+16','R+18','R+20','R+22','R+24', & - 'RRR ','73 '/ + character*4 rpt(0:63) equivalence (ipk0,ipk) nargs=iargc() @@ -40,13 +35,25 @@ program msk32d idf1=nf1-1500 idf2=nf2-1500 + do i=0,30 + if( i.lt.5 ) then + write(rpt(i),'(a1,i2.2,a1)') '-',abs(i-5) + write(rpt(i+31),'(a2,i2.2,a1)') 'R-',abs(i-5) + else + write(rpt(i),'(a1,i2.2,a1)') '+',i-5 + write(rpt(i+31),'(a2,i2.2,a1)') 'R+',i-5 + endif + enddo + rpt(62)='RRR ' + rpt(63)='73 ' + ! Generate the test messages twopi=8.0*atan(1.0) nsym=32 freq=1500.0 dphi0=twopi*(freq-500.0)/12000.0 dphi1=twopi*(freq+500.0)/12000.0 - do imsg=0,31 + do imsg=0,63 i=index(hiscall," ") msg="<"//mycall//" "//hiscall(1:i-1)//"> "//rpt(imsg) call fmtmsg(msg,iz) @@ -81,7 +88,7 @@ program msk32d nfft=min(2**n,1024*1024) call analytic(dd,npts,nfft,c0) !Convert to analytic signal sbest=0. - do imsg=0,31 !Try all short messages + do imsg=0, 63 !Try all short messages do idf=idf1,idf2,10 !Frequency dither a(1)=-idf a(2:3)=0. @@ -140,7 +147,7 @@ program msk32d s0=s0-ave s1=sbest-ave s2=0. - do i=0,31 + do i=0,63 if(i.ne.ibest .and. s0(i).gt.s2) s2=s0(i) write(15,1020) i,idf,jpk/12000.0,s0(i) 1020 format(2i6,2f10.2) diff --git a/lib/msk32d_ldpc.f90 b/lib/msk32d_ldpc.f90 new file mode 100644 index 000000000..ff2138308 --- /dev/null +++ b/lib/msk32d_ldpc.f90 @@ -0,0 +1,171 @@ +program msk32d + + parameter (NZ=15*12000,NZ0=262144) + parameter (NSPM=6*32) + complex c0(0:NZ0-1) + complex c(0:NZ-1) + complex cmsg(0:NSPM-1,0:63) + complex z + real a(3) + real p0(0:NSPM-1) + real p(0:NSPM-1) + real s0(0:63) + real dd(NZ) + integer itone(144) + integer ihdr(11) + integer ipk0(1) + integer*2 id2(NZ) + character*22 msg,msgsent + character mycall*8,hiscall*6,arg*12,infile*80,datetime*13 + character*4 rpt(0:63) + equivalence (ipk0,ipk) + + nargs=iargc() + if(nargs.lt.5) then + print*,'Usage: msk32d Call_1 Call_2 f1 f2 file1 [file2 ...]' + print*,'Example: msk32d K9AN K1JT 1500 1500 fort.61' + go to 999 + endif + call getarg(1,mycall) + call getarg(2,hiscall) + call getarg(3,arg) + read(arg,*) nf1 + call getarg(4,arg) + read(arg,*) nf2 + idf1=nf1-1500 + idf2=nf2-1500 + + do i=0,30 + if( i.lt.5 ) then + write(rpt(i),'(a1,i2.2,a1)') '-',abs(i-5) + write(rpt(i+31),'(a2,i2.2,a1)') 'R-',abs(i-5) + else + write(rpt(i),'(a1,i2.2,a1)') '+',i-5 + write(rpt(i+31),'(a2,i2.2,a1)') 'R+',i-5 + endif + enddo + rpt(62)='RRR ' + rpt(63)='73 ' + +! Generate the test messages + twopi=8.0*atan(1.0) + nsym=32 + freq=1500.0 + dphi0=twopi*(freq-500.0)/12000.0 + dphi1=twopi*(freq+500.0)/12000.0 + do imsg=0,63 + i=index(hiscall," ") + msg="<"//mycall//" "//hiscall(1:i-1)//"> "//rpt(imsg) + call fmtmsg(msg,iz) + ichk=0 + call genmsk32(msg,msgsent,ichk,itone,itype) + + phi=0.0 + k=0 + do i=1,nsym + dphi=dphi0 + if(itone(i).eq.1) dphi=dphi1 + do j=1,6 + x=cos(phi) + y=sin(phi) + cmsg(k,imsg)=cmplx(x,y) + k=k+1 + phi=phi+dphi + if(phi.gt.twopi) phi=phi-twopi + enddo + enddo + enddo + +! Process the specified files + nfiles=nargs-4 + do ifile=1,nfiles !Loop over all files + call getarg(ifile+4,infile) + open(10,file=infile,access='stream',status='old') + read(10) ihdr,id2 + dd=0.03*id2 + npts=NZ + n=log(float(npts))/log(2.0) + 1.0 + nfft=min(2**n,1024*1024) + call analytic(dd,npts,nfft,c0) !Convert to analytic signal + sbest=0. + do imsg=0, 63 !Try all short messages + do idf=idf1,idf2,10 !Frequency dither + a(1)=-idf + a(2:3)=0. + call twkfreq(c0,c,npts,12000.0,a) + smax=0. + p=0. + fac=1.0/192 + do j=0,npts-NSPM,2 + z=fac*dot_product(c(j:j+NSPM-1),cmsg(0:NSPM-1,imsg)) + s=real(z)**2 + aimag(z)**2 + k=mod(j,NSPM) + p(k)=p(k)+s +! if(imsg.eq.30) write(13,1010) j/12000.0,s,k +!1010 format(2f12.6,i5) + if(s.gt.smax) then + smax=s + jpk=j + f0=idf + if(smax.gt.sbest) then + sbest=smax + p0=p + ibest=imsg + endif + endif + enddo + s0(imsg)=smax + enddo + enddo + + ipk0=maxloc(p0) + ps=0. + sq=0. + ns=0 + pmax=0. + do i=0,NSPM-1,2 + j=ipk-i + if(j.gt.96) j=j-192 + if(j.lt.-96) j=j+192 + if(abs(j).gt.4) then + ps=ps+p0(i) + sq=sq+p0(i)**2 + ns=ns+1 + endif + enddo + avep=ps/ns + rmsp=sqrt(sq/ns - avep*avep) + p0=(p0-avep)/rmsp + p1=maxval(p0) + + do i=0,NSPM-1,2 + write(14,1030) i,i/12000.0,p0(i) +1030 format(i5,f10.6,f10.3) + enddo + + ave=(sum(s0)-sbest)/31 + s0=s0-ave + s1=sbest-ave + s2=0. + do i=0,63 + if(i.ne.ibest .and. s0(i).gt.s2) s2=s0(i) + write(15,1020) i,idf,jpk/12000.0,s0(i) +1020 format(2i6,2f10.2) + enddo + + i=index(infile,".wav") + datetime=infile(i-13:i-1) + r1=s1/s2 + r2=r1+p1 + msg=" " +! if(r1.gt.2.2 .or. p1.gt.7.0) then + if(r2.gt.10.0) then + i=index(hiscall," ") + msg="<"//mycall//" "//hiscall(1:i-1)//"> "//rpt(ibest) + call fmtmsg(msg,iz) + endif + write(*,1040) datetime,r1,p1,r2,msg +1040 format(a13,3f7.1,2x,a22) + enddo + +999 end program msk32d