diff --git a/CMakeLists.txt b/CMakeLists.txt index ad128730f..e1aa25a3d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1029,6 +1029,9 @@ target_link_libraries (jt65sim wsjt_fort wsjt_cxx) add_executable (qra64sim lib/qra/qra64/qra64sim.f90 wsjtx.rc) target_link_libraries (qra64sim wsjt_fort wsjt_cxx) +add_executable (msk32d lib/msk32d.f90 wsjtx.rc) +target_link_libraries (msk32d wsjt_fort wsjt_cxx) + add_executable (jt9sim lib/jt9sim.f90 wsjtx.rc) target_link_libraries (jt9sim wsjt_fort wsjt_cxx) @@ -1202,7 +1205,7 @@ install (TARGETS udp_daemon message_aggregator BUNDLE DESTINATION ${WSJT_BIN_DESTINATION} COMPONENT runtime ) -install (TARGETS jt9 jt65code qra64code qra64sim jt9code jt4code wsprd +install (TARGETS jt9 jt65code qra64code qra64sim jt9code jt4code wsprd msk32d RUNTIME DESTINATION ${WSJT_BIN_DESTINATION} COMPONENT runtime BUNDLE DESTINATION ${WSJT_BIN_DESTINATION} COMPONENT runtime ) diff --git a/lib/msk32d.f90 b/lib/msk32d.f90 index b1f26bbec..24238ae09 100644 --- a/lib/msk32d.f90 +++ b/lib/msk32d.f90 @@ -1,44 +1,66 @@ program msk32d - parameter (NZ=15*12000) + parameter (NZ=15*12000,NZ0=262144) parameter (NSPM=6*32) - complex c0(0:NZ-1) + complex c0(0:NZ0-1) complex c(0:NZ-1) complex cmsg(0:NSPM-1,0:31) complex z - real*8 twopi,freq,phi,dphi0,dphi1,dphi real a(3) real p0(0:NSPM-1) real p(0:NSPM-1) + real s0(0:31) + 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: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 '/ + equivalence (ipk0,ipk) - twopi=8.d0*atan(1.d0) + 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 + +! Generate the test messages + twopi=8.0*atan(1.0) nsym=32 - freq=1500.d0 - dphi0=twopi*(freq-500.d0)/12000.d0 - dphi1=twopi*(freq+500.d0)/12000.d0 - + freq=1500.0 + dphi0=twopi*(freq-500.0)/12000.0 + dphi1=twopi*(freq+500.0)/12000.0 do imsg=0,31 - msg=" "//rpt(imsg) + 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.d0 + phi=0.0 k=0 do i=1,nsym dphi=dphi0 if(itone(i).eq.1) dphi=dphi1 do j=1,6 - xphi=phi - x=cos(xphi) - y=sin(xphi) + x=cos(phi) + y=sin(phi) cmsg(k,imsg)=cmplx(x,y) k=k+1 phi=phi+dphi @@ -47,46 +69,96 @@ program msk32d enddo enddo - read(61) npts,c0(0:npts-1) +! 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,31 !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 - sbest=0. - do imsg=0,31 -! do idf=-50,50,10 - ! imsg=30 - idf=-5 - 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 - endif + 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 - write(*,1020) imsg,idf,jpk/12000.0,smax - write(15,1020) imsg,idf,jpk/12000.0,smax -1020 format(2i6,2f10.2) + 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,31 + 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 - imsg=30 - do i=0,NSPM-1,2 - write(14,1030) i,i/12000.0,p0(i) -1030 format(i5,f10.6,f10.3) - enddo - - end program msk32d +999 end program msk32d