diff --git a/lib/JTMSKsim.f90 b/lib/JTMSKsim.f90 index c984deb1f..23e430d0f 100644 --- a/lib/JTMSKsim.f90 +++ b/lib/JTMSKsim.f90 @@ -1,134 +1,93 @@ program JTMSKsim use wavhdr - parameter (NSPM=1404) - parameter (NFFT=256*1024) parameter (NMAX=15*12000) - type(hdr) h - complex cb11(0:NSPM-1) - complex cmsg(0:NSPM-1) - complex c(0:NFFT-1) - complex c1(0:NSPM-1) - complex c2(0:NSPM-1) - integer i4tone(NSPM) - integer*2 iwave(NMAX) - real w(0:255) - character*8 arg + real pings(0:NMAX-1) + character arg*8,msg*22,msgsent*22,fname*40 + character*3 rpt(0:7) + complex cmsg(0:1404-1) !Waveform of message (once) + complex cwave(0:NMAX-1) !Simulated received waveform + real*8 dt,twopi,freq,phi,dphi0,dphi1,dphi + type(hdr) h !Header for .wav file + integer*2 iwave(0:NMAX-1) + integer itone(234) !Message bits + integer b11(11) !Barker-11 code + data b11/1,1,1,0,0,0,1,0,0,1,0/ + data rpt /'26 ','27 ','28 ','R26','R27','R28','RRR','73 '/ nargs=iargc() - if(nargs.ne.2) then - print*,'Usage: JTMSK freq dB' + if(nargs.ne.5) then + print*,'Usage: JTMSKsim message freq width snr nfiles' + print*,' ' + print*,'Examples: JTMSKsim "K1ABC W9XYZ EN37" 1500 0.05 0.0 1' + print*,' JTMSKsim " R26" 1500 0.01 0.0 10' go to 999 endif - call getarg(1,arg) - read(arg,*) fmid + call getarg(1,msg) call getarg(2,arg) + read(arg,*) freq + call getarg(3,arg) + read(arg,*) width + call getarg(4,arg) read(arg,*) snrdb - - open(10,file='JTMSKcode.out',status='old') - do j=1,234 - read(10,*) junk,i4tone(j) - enddo - close(10) - - npts=NMAX - h=default_header(12000,npts) - twopi=8.0*atan(1.0) - dt=1.0/12000.0 - df=12000.0/NFFT - k=-1 - phi=0. - phi0=0. + call getarg(5,arg) + read(arg,*) nfiles sig=10.0**(0.05*snrdb) + twopi=8.d0*atan(1.d0) + h=default_header(12000,NMAX) - call init_random_seed() ! seed Fortran RANDOM_NUMBER generator - call sgran() ! see C rand generator (used in gran) + ichk=0 + call genmsk(msg,ichk,msgsent,itone,itype) !Check message type + if(itype.lt.1 .or. itype.gt.7) then + print*,'Illegal message' + go to 999 + endif + dt=1.d0/12000.d0 !Sample interval + dphi0=twopi*(freq-500.d0)*dt !Phase increment, lower tone + dphi1=twopi*(freq+500.d0)*dt !Phase increment, upper tone -! Generate the Tx waveform - cb11=0. - do j=1,234 - dphi=twopi*(fmid-500.0)*dt - if(i4tone(j).eq.1) dphi=twopi*(fmid+500.0)*dt - dphi0=twopi*1000.0*dt - if(i4tone(j).eq.1) dphi0=twopi*2000.0*dt + nsym=234 + if(itype.eq.7) nsym=35 + nspm=6*nsym !Samples per message + k=-1 + phi=0.d0 + do j=1,nsym + dphi=dphi0 + if(itone(j).eq.1) dphi=dphi1 do i=1,6 k=k+1 - phi=phi+dphi + phi=phi + dphi if(phi.gt.twopi) phi=phi-twopi - cmsg(k)=cmplx(cos(phi),sin(phi)) - phi0=phi0+dphi0 - if(phi0.gt.twopi) phi0=phi0-twopi -! if((k.ge.1 .and. k.le.66) .or. (k.ge.283 .and. k.le.348) .or. & -! (k.ge.769 .and.k.le.834)) cb11(k)=cmplx(cos(phi0),sin(phi0)) - if(k.ge.1 .and. k.le.66) cb11(k)=cmplx(cos(phi0),sin(phi0)) -! write(11,3001) k*dt,cmsg(k),phi -!3001 format(4f12.6) + xphi=phi + cmsg(k)=cmplx(cos(xphi),sin(xphi)) enddo enddo -! Generate noise with B=2500 Hz, rms=1.0 - c=0. - ia=nint(250.0/df) - ib=nint(2759.0/df) - do i=ia,ib - x=gran() - y=gran() - c(i)=cmplx(x,y) - enddo - call four2a(c,NFFT,1,1,1) - sq=0. - do i=0,npts-1 - sq=sq + real(c(i))**2 + aimag(c(i))**2 - enddo - rms=sqrt(0.5*sq/npts) - c=c/rms + call makepings(pings,NMAX,width,sig) - i0=3*12000 - ncopy=NSPM -! ncopy=0.5*NSPM - c(i0:i0+ncopy-1)=c(i0:i0+ncopy-1) + sig*cmsg(0:ncopy-1) - do i=1,npts - iwave(i)=100.0*real(c(i)) - enddo + do ifile=1,nfiles !Loop over requested number of files + write(fname,1002) ifile !Output filename +1002 format('000000_',i4.4) + open(10,file=fname(1:11)//'.wav',access='stream',status='unknown') - open(12,file='150901_000000.wav',status='unknown',access='stream') - write(12) h,iwave(1:npts) - - smax=0. - fpk=0. - jpk=0 - nfft1=256 - w=0. - do i=0,65 - w(i)=sin(i*twopi/132.0) - enddo - - do ia=0,npts-nfft1 - c1(0:nfft1-1)=c(ia:ia+nfft1-1)*conjg(cb11(0:nfft1-1)) - c2(0:nfft1-1)=w(0:nfft1-1)*c1(0:nfft1-1) - call four2a(c2,nfft1,1,-1,1) - do i=0,nfft1-1 -! write(21,1100) i,c1(i) -!1100 format(i6,2f12.3) + fac=sqrt(6000.0/2500.0) + j=-1 + do i=0,NMAX-1 + j=j+1 + if(j.ge.6*nsym) j=j-6*nsym + xx=0.707*gran() + yy=0.707*gran() + cwave(i)=pings(i)*cmsg(j) + fac*cmplx(xx,yy) + iwave(i)=30.0*real(cwave(i)) +! write(88,3003) i,i/12000.d0,cwave(i) +!3003 format(i8,f12.6,2f10.3) enddo - df1=12000.0/nfft1 - do i=-20,20 - j=i - if(i.lt.0) j=j+nfft1 - f=j*df1 - if(i.lt.0) f=f-12000.0 - s=1.e-3*(real(c2(j))**2 + aimag(c2(j))**2) - if(abs(ia-i0).lt.1404) write(22,1110) f,c2(j),s,ia -1110 format(4f12.3,i6) - if(s.gt.smax) then - smax=s - jpk=ia - fpk=f - endif - enddo + write(10) h,iwave !Save the .wav file + close(10) + + call jtmsk_short(cwave,NMAX,msg) + enddo - print*,smax,jpk*dt,fpk 999 end program JTMSKsim - diff --git a/lib/Makefile.jt65 b/lib/Makefile.jt65 index e64f874c0..7409b5c19 100644 --- a/lib/Makefile.jt65 +++ b/lib/Makefile.jt65 @@ -36,11 +36,11 @@ OBJS1 = astrosub.o astro0.o astro.o sun.o coord.o tmoonsub.o \ four2a.o grid2deg.o wisdom.o \ symspec.o analytic.o db.o \ encode232.o interleave9.o\ - entail.o fano232.o gran.o sync9.o decjt9.o \ - fil3.o decoder.o timer.o exp_decode65.o fqso_first.o \ + entail.o fano232.o gran.o sync9.o \ + fil3.o decoder.o fqso_first.o \ twkfreq.o symspec2.o shell.o sync65.o peakup.o slope.o xcor.o\ - fillcom.o chkss2.o zplot9.o flat1.o flat2.o \ - jt65a.o symspec65.o flat65.o ccf65.o decode65a.o \ + chkss2.o zplot9.o flat1.o flat2.o \ + symspec65.o flat65.o ccf65.o decode65a.o \ filbig.o fil6521.o afc65b.o decode65b.o setup65.o \ extract.o fchisq65.o demod64a.o chkhist.o interleave63.o ccf2.o \ move.o indexx.o graycode65.o twkfreq65.o smo.o smo121.o \ diff --git a/lib/genmsk.f90 b/lib/genmsk.f90 index 7eccb443d..63a49eef3 100644 --- a/lib/genmsk.f90 +++ b/lib/genmsk.f90 @@ -3,16 +3,18 @@ subroutine genmsk(msg0,ichk,msgsent,i4tone,itype) ! Encode a JTMSK message ! Input: ! - msg0 requested message to be transmitted -! - ichk if nonzero, return only msgsent +! - ichk if ichk=1, return only msgsent +! if ichk.ge.10000, set imsg=ichk-10000 for short msg ! - msgsent message as it will be decoded ! - i4tone array of audio tone values, 0 or 1 ! - itype message type -! 1 = standard message +! 1 = standard message "Call_1 Call_2 Grid/Rpt" ! 2 = type 1 prefix ! 3 = type 1 suffix ! 4 = type 2 prefix ! 5 = type 2 suffix ! 6 = free text (up to 13 characters) +! 7 = short message " Rpt" use iso_c_binding, only: c_loc,c_size_t use packjt @@ -51,15 +53,15 @@ subroutine genmsk(msg0,ichk,msgsent,i4tone,itype) enddo if(message(1:1).eq.'<') then - call genmsk_short(message,msgsent,i4tone,itype) + call genmsk_short(message,msgsent,ichk,i4tone,itype) if(itype.lt.0) go to 999 - i4tone(38)=-37 + i4tone(36)=-35 go to 999 endif call packmsg(message,i4Msg6BitWords,itype) !Pack into 12 6-bit bytes call unpackmsg(i4Msg6BitWords,msgsent) !Unpack to get msgsent - if(ichk.ne.0) go to 999 + if(ichk.eq.1) go to 999 call entail(i4Msg6BitWords,i1Msg8BitBytes) !Add tail, make 8-bit bytes ihash=nhash(c_loc(i1Msg8BitBytes),int(9,c_size_t),146) ihash=2*iand(ihash,32767) !Generate the CRC diff --git a/lib/genmsk_short.f90 b/lib/genmsk_short.f90 index b4860f1cf..b0292abfb 100644 --- a/lib/genmsk_short.f90 +++ b/lib/genmsk_short.f90 @@ -1,14 +1,21 @@ -subroutine genmsk_short(msg,msgsent,itone,itype) +subroutine genmsk_short(msg,msgsent,ichk,itone,itype) use hashing - character*22 msg0,msg,msgsent + character*22 msg,msgsent character*3 crpt,rpt(0:7) - integer itone(37) + logical first + integer itone(35) + integer ig24(0:4096-1) !Codewords for Golay (24,12) code integer b11(11) - integer golay_15_3(0:7) data b11/1,1,1,0,0,0,1,0,0,1,0/ !Barker 11 code - data rpt /'26 ','27 ','28 ','R26','R27','R28','RRR','73 '/ - data golay_15_3/00000,07025,11704,14025,19164,20909,26468,31765/ + data rpt /'26 ','27 ','28 ','R26','R27','R28','RRR','73 '/ + data first/.true./ + save first,ig24 + + if(first) then + call golay24_table(ig24) !Define the Golay(24,12) codewords + first=.false. + endif itype=-1 msgsent='*** bad message ***' @@ -22,25 +29,27 @@ subroutine genmsk_short(msg,msgsent,itone,itype) enddo go to 900 -10 itone(1:11)=b11 - irpt=i - ncodeword=golay_15_3(irpt) !Save the 15-bit codeword for report - do i=26,12,-1 !Insert codeword into itone array - itone(i)=iand(ncodeword,1) - ncodeword=ncodeword/2 - enddo - call hash(msg(2:i1-1),i1-2,ihash) - ihash=iand(ihash,1023) !10-bit hash for the two callsigns - n=ihash - do i=36,27,-1 - itone(i)=iand(n,1) +10 irpt=i !Report index, 0-7 + if(ichk.lt.10000) then + call hash(msg(2:i1-1),i1-2,ihash) + ihash=iand(ihash,511) !9-bit hash for the two callsigns + ig=8*ihash + irpt !12-bit message information + else + ig=ichk-10000 + endif + ncodeword=ig24(ig) + itone(1:11)=b11 !Insert the Barker-11 code + n=2**24 + do i=12,35 !Insert codeword into itone array n=n/2 - enddo - n=count(itone(1:36).eq.0) - itone(37)=1 - if(iand(n,1).eq.1) itone(37)=0 + itone(i)=0 + if(iand(ncodeword,n).ne.0) itone(i)=1 + enddo msgsent=msg itype=7 + n=count(itone(1:35).eq.0) + if(mod(n,2).ne.0) stop 'Parity error in genmsk_short.' + 900 return end subroutine genmsk_short diff --git a/lib/jtmsk_short.f90 b/lib/jtmsk_short.f90 new file mode 100644 index 000000000..6b0c50d22 --- /dev/null +++ b/lib/jtmsk_short.f90 @@ -0,0 +1,145 @@ +subroutine jtmsk_short(cdat,npts,msg) + + parameter (NMAX=15*12000,NSAVE=100) + character*22 msg,decoded,msgsent + character*3 rpt(0:7) + complex cdat(0:npts-1) + complex cw(0:209,0:4096) !Waveforms of possible messages + complex cb11(0:65) !Complex waveform of Barker 11 + complex z1,z2a,z2b + real*8 dt,twopi,freq,phi,dphi0,dphi1,dphi + real r1(0:NMAX-1) + real r2(0:4096) + real r1save(NSAVE) + integer*8 count0,count1,clkfreq + integer itone(234) !Message bits + integer jgood(NSAVE) + integer indx(NSAVE) + logical first + data rpt /'26 ','27 ','28 ','R26','R27','R28','RRR','73 '/ + data first/.true./ + save first,cw + +! msg=' 73 ' + + if(first) then + dt=1.d0/12000.d0 + twopi=8.d0*atan(1.d0) + freq=1500.d0 + dphi0=twopi*(freq-500.d0)*dt !Phase increment, lower tone + dphi1=twopi*(freq+500.d0)*dt !Phase increment, upper tone + nsym=35 !Number of symbols + nspm=6*nsym !Samples per message + + do imsg=0,4096 !Generate all possible message waveforms + ichk=0 + if(imsg.lt.4096) ichk=10000+imsg + call genmsk(msg,ichk,msgsent,itone,itype) !Encode the message + k=-1 + phi=0.d0 + do j=1,nsym + dphi=dphi0 + if(itone(j).eq.1) dphi=dphi1 + do i=1,6 + k=k+1 + phi=phi + dphi + if(phi.gt.twopi) phi=phi-twopi + xphi=phi + cw(k,imsg)=cmplx(cos(xphi),sin(xphi)) + enddo + enddo + enddo + first=.false. + endif + + r1thresh=0.40 + r2thresh=0.50 + rmax=0.9 + ngood=0 + nbad=0 + maxdecodes=999 + + cb11=cw(0:65,0) + call system_clock(count0,clkfreq) + + r1max=0. + do j=0,npts-210 !Find the B11 sync vectors + z1=0. + ss=0. + do i=0,65 + ss=ss + real(cdat(j+i))**2 + aimag(cdat(j+i))**2 + z1=z1 + cdat(j+i)*conjg(cb11(i)) !Signal matching B11 + enddo + ss=sqrt(ss/66.0)*66.0 + r1(j)=abs(z1)/(0.908*ss) !Goodness-of-fit to B11 + if(r1(j).gt.r1max) then + r1max=r1(j) + jpk=j + endif + enddo + + k=0 + do j=1,npts-211 + if(r1(j).gt.r1thresh .and. r1(j).ge.r1(j-1) .and. r1(j).ge.r1(j+1) ) then + k=k+1 + jgood(k)=j + r1save(k)=r1(j) + if(k.ge.NSAVE) exit + endif + enddo + kmax=k + call indexx(r1save,kmax,indx) + + r2bad=0. + do kk=1,kmax + k=indx(kmax+1-kk) + j=jgood(k) + u1=0. + u2=0. + r2max=0. + do imsg=0,4096 + ssa=0. + ssb=0. + do i=0,209 + ssa=ssa + real(cdat(j+i))**2 + aimag(cdat(j+i))**2 + ssb=ssb + real(cdat(j+i-144))**2 + aimag(cdat(j+i-144))**2 + enddo + z2a=dot_product(cw(0:209,imsg),cdat(j:j+209)) + z2b=dot_product(cw(0:65,imsg),cdat(j:j+65)) + & + dot_product(cw(66:209,imsg),cdat(j-144:j-1)) + ssa=sqrt(ssa/210.0)*210.0 + ssb=sqrt(ssb/210.0)*210.0 + r2(imsg)=max(abs(z2a)/(0.908*ssa),abs(z2b)/(0.908*ssb)) + if(r2(imsg).gt.r2max) then + r2max=r2(imsg) + ibest=imsg + u2=u1 + u1=r2max + endif + enddo + r1or2=r1(j)/r2max + if(r2max.ge.r2thresh .and. u2/u1.lt.rmax .and. r1or2.ge.0.91) then + t=j/12000.0 + n=0 + irpt=iand(ibest,7) + decoded="<...> "//rpt(irpt) + if(r2max.eq.r2(4096)) then + n=1 + decoded=msg(1:14)//rpt(irpt) + endif + if(n.eq.0) nbad=nbad+1 + if(n.eq.1) ngood=ngood+1 + if(n.eq.0 .and. r2max.gt.r2bad) r2bad=r2max + write(52,3020) k,t,ibest,r1(j),r2max,u2/u1,r1or2,n,decoded +3020 format(i3,f9.4,i5,4f7.2,i2,1x,a22) + if(ngood+nbad.ge.maxdecodes) exit + endif + enddo + + call system_clock(count1,clkfreq) + t=float(count1-count0)/float(clkfreq) + print "('Decode: ',f6.3)",t + print "('Worst false decode:',f6.3)",r2bad + print "('Good:',i3,' Bad:',i3)",ngood,nbad + return +end subroutine jtmsk_short diff --git a/lib/makepings.f90 b/lib/makepings.f90 new file mode 100644 index 000000000..1e2b12d2a --- /dev/null +++ b/lib/makepings.f90 @@ -0,0 +1,27 @@ +subroutine makepings(pings,npts,width,sig) + + real pings(npts) + real*8 t + real t0(14) + + iping0=-999 + dt=1.0/12000.0 + do i=1,14 + t0(i)=i !Make pings at t=1, 2, ... 14 s. + enddo + w=width + amp=sig + + do i=1,npts + iping=min(max(1,i/12000),14) + t=(i*dt-t0(iping))/w + if(t.lt.0.d0 .and. t.lt.10.0) then + fac=0. + else + fac=2.718*t*dexp(-t) + endif + pings(i)=fac*amp + enddo + + return +end subroutine makepings