diff --git a/CMakeLists.txt b/CMakeLists.txt index ed1f5c74d..3b6cc15a7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1128,6 +1128,9 @@ target_link_libraries (sumsim wsjt_fort wsjt_cxx) add_executable (q65sim lib/qra/q65/q65sim.f90) target_link_libraries (q65sim wsjt_fort wsjt_cxx) +add_executable (q65code lib/qra/q65/q65code.f90) +target_link_libraries (q65code wsjt_fort wsjt_cxx) + add_executable (test_q65 lib/test_q65.f90) target_link_libraries (test_q65 wsjt_fort wsjt_cxx) @@ -1563,7 +1566,8 @@ install (TARGETS jt9 wsprd fmtave fcal fmeasure ) if(WSJT_BUILD_UTILS) -install (TARGETS ft8code jt65code jt9code jt4code msk144code fst4sim q65sim +install (TARGETS ft8code jt65code jt9code jt4code msk144code + q65code fst4sim q65sim RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} COMPONENT runtime BUNDLE DESTINATION ${CMAKE_INSTALL_BINDIR} COMPONENT runtime ) diff --git a/lib/qra/q65/q65code.f90 b/lib/qra/q65/q65code.f90 new file mode 100644 index 000000000..4d16ec914 --- /dev/null +++ b/lib/qra/q65/q65code.f90 @@ -0,0 +1,236 @@ +module gf64math +! add and subtract in GF(2^6) based on primitive polynomial x^6+x+1 + + implicit none + integer, private :: gf64log(0:63) + integer, private :: gf64antilog(0:62) + +! table of the logarithms of the elements of GF(M) (log(0) never used) + data gf64log/ & + -1, 0, 1, 6, 2, 12, 7, 26, 3, 32, & + 13, 35, 8, 48, 27, 18, 4, 24, 33, 16, & + 14, 52, 36, 54, 9, 45, 49, 38, 28, 41, & + 19, 56, 5, 62, 25, 11, 34, 31, 17, 47, & + 15, 23, 53, 51, 37, 44, 55, 40, 10, 61, & + 46, 30, 50, 22, 39, 43, 29, 60, 42, 21, & + 20, 59, 57, 58/ + +! table of GF(M) elements given their logarithm + data gf64antilog/ & + 1, 2, 4, 8, 16, 32, 3, 6, 12, 24, & + 48, 35, 5, 10, 20, 40, 19, 38, 15, 30, & + 60, 59, 53, 41, 17, 34, 7, 14, 28, 56, & + 51, 37, 9, 18, 36, 11, 22, 44, 27, 54, & + 47, 29, 58, 55, 45, 25, 50, 39, 13, 26, & + 52, 43, 21, 42, 23, 46, 31, 62, 63, 61, & + 57, 49, 33/ + +contains + + integer function gf64_add(i1,i2) + implicit none + integer::i1 + integer::i2 + gf64_add=iand(ieor(i1,i2),63) + end function gf64_add + + integer function gf64_mult(i1,i2) + implicit none + integer::i1 + integer::i2 + integer::j + + if(i1.eq.0 .or. i2.eq.0) then + gf64_mult=0 + elseif(i1.eq.1) then + gf64_mult=i2 + elseif(i2.eq.1) then + gf64_mult=i1 + else + j=mod(gf64log(i1)+gf64log(i2),63) + gf64_mult=gf64antilog(j) + endif + end function gf64_mult + +end module gf64math + +module q65_generator + + integer generator(15,50) + data generator/ & + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, & + 0,20, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, & + 0,20, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, & + 0,20, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, & + 0,20, 0, 1, 1, 0, 0, 0,10, 0, 0, 0, 0, 1, 0, & + 0,20, 0, 1, 1, 0, 0, 0,10, 0, 0, 0,44, 1, 0, & + 0,20, 0, 1, 1, 0, 0, 0,10, 1, 0, 0,44, 1, 0, & + 0,20, 0, 1, 1, 0, 0, 0,10, 1, 0, 0,44, 1,14, & + 0,20, 0, 1, 1, 0, 0, 0,10, 1,31, 0,44, 1,14, & + 0,20, 0, 1, 1,33, 0, 0,10, 1,31, 0,44, 1,14, & + 56,20, 0, 1, 1,33, 0, 0,10, 1,31, 0,44, 1,14, & + 56,20, 0, 1, 1,33, 0, 1,10, 1,31, 0,44, 1,14, & + 56, 1, 0, 1, 1,33, 0, 1,10, 1,31, 0,44, 1,14, & + 56, 1, 0, 1, 1,33, 0, 1,10, 1,31,36,44, 1,14, & + 56, 1, 0, 1, 1,33, 0, 1,43, 1,31,36,44, 1,14, & + 56, 1, 0, 1, 1,33, 0, 1,43,17,31,36,44, 1,14, & + 56, 1, 0, 1, 1,33, 0, 1,43,17,31,36,36, 1,14, & + 56, 1, 0, 1, 1,33,53, 1,43,17,31,36,36, 1,14, & + 56, 1, 0,35, 1,33,53, 1,43,17,31,36,36, 1,14, & + 56, 1, 0,35, 1,33,53, 1,43,17,30,36,36, 1,14, & + 56, 1, 0,35, 1,33,53,52,43,17,30,36,36, 1,14, & + 56, 1, 0,35, 1,32,53,52,43,17,30,36,36, 1,14, & + 56, 1,60,35, 1,32,53,52,43,17,30,36,36, 1,14, & + 56, 1,60,35, 1,32,53,52,43,17,30,36,36,49,14, & + 56, 1,60,35, 1,32,53,52,43,17,30,36,37,49,14, & + 56, 1,60,35,54,32,53,52,43,17,30,36,37,49,14, & + 56, 1,60,35,54,32,53,52, 1,17,30,36,37,49,14, & + 1, 1,60,35,54,32,53,52, 1,17,30,36,37,49,14, & + 1, 0,60,35,54,32,53,52, 1,17,30,36,37,49,14, & + 1, 0,60,35,54,32,53,52, 1,17,30,37,37,49,14, & + 1, 0,61,35,54,32,53,52, 1,17,30,37,37,49,14, & + 1, 0,61,35,54,32,53,52, 1,48,30,37,37,49,14, & + 1, 0,61,35,54,32,53,52, 1,48,30,37,37,49,15, & + 1, 0,61,35,54, 0,53,52, 1,48,30,37,37,49,15, & + 1, 0,61,35,54, 0,52,52, 1,48,30,37,37,49,15, & + 1, 0,61,35,54, 0,52,52, 1,48,30,37,37, 0,15, & + 1, 0,61,35,54, 0,52,34, 1,48,30,37,37, 0,15, & + 1, 0,61,35,54, 0,52,34, 1,48,30,37, 0, 0,15, & + 1, 0,61,35,54, 0,52,34, 1,48,30,20, 0, 0,15, & + 1, 0, 0,35,54, 0,52,34, 1,48,30,20, 0, 0,15, & + 1, 0, 0,35,54, 0,52,34, 1, 0,30,20, 0, 0,15, & + 0, 0, 0,35,54, 0,52,34, 1, 0,30,20, 0, 0,15, & + 0, 0, 0,35,54, 0,52,34, 1, 0,38,20, 0, 0,15, & + 0, 0, 0,35, 0, 0,52,34, 1, 0,38,20, 0, 0,15, & + 0, 0, 0,35, 0, 0,52, 0, 1, 0,38,20, 0, 0,15, & + 0, 0, 0,35, 0, 0,52, 0, 1, 0,38,20, 0, 0, 0, & + 0, 0, 0,35, 0, 0,52, 0, 0, 0,38,20, 0, 0, 0, & + 0, 0, 0,35, 0, 0,52, 0, 0, 0,38, 0, 0, 0, 0, & + 0, 0, 0, 0, 0, 0,52, 0, 0, 0,38, 0, 0, 0, 0, & + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,38, 0, 0, 0, 0/ + +end module q65_generator + +subroutine q65_encode(message,codeword) + use gf64math + use q65_generator + integer message(15) + integer codeword(65) + integer i,j + + codeword=0 + codeword(1:15)=message + do i=1,15 + do j=16,65 + codeword(j)=gf64_add(codeword(j),gf64_mult(message(i),generator(i,j-15))) + enddo + enddo + + return +end + +subroutine get_q65crc12(mc2,ncrc1,ncrc2) +! + character c12*12,c6*6 + integer*1 mc(90),mc2(90),tmp(6) + integer*1 r(13),p(13) + integer ncrc +! polynomial for 12-bit CRC 0xF01 + data p/1,1,0,0,0,0,0,0,0,1,1,1,1/ + +! flip bit order of each 6-bit symbol for consistency with Nico's calculation + do i=0,14 + tmp=mc2(i*6+1:i*6+6) + mc(i*6+1:i*6+6)=tmp(6:1:-1) + enddo + +! divide by polynomial + r=mc(1:13) + do i=0,77 + r(13)=mc(i+13) + r=mod(r+r(1)*p,2) + r=cshift(r,1) + enddo + + write(c6,'(6b1)') r(6:1:-1) + read(c6,'(b6.6)') ncrc1 + read(c6,'(6b1)') mc2(79:84) + write(c6,'(6b1)') r(12:7:-1) + read(c6,'(b6.6)') ncrc2 + read(c6,'(6b1)') mc2(85:90) + +end subroutine get_q65crc12 + +program q65code + + use packjt77 + + implicit none + character msg*37,msgsent*27,arg*12,c12*12,c6*6 + integer message(15) + integer codeword(65), shortcodeword(63) + integer isync(22) + integer itone(85) + integer i,j,k + integer nargs,ncrc1,ncrc2,ncrc + integer i3,n3 + integer*1 mbits(90) + integer msymbols(15) + character*77 c77 + data isync/1,9,12,13,15,22,23,26,27,33,35,38,46,50,55,60,62,66,69,74,76,85/ + nargs=iargc() + if(nargs .ne. 1) then + print*,'Usage: q65sf "msg"' + goto 999 + endif + call getarg(1,msg) + + i3=-1 + n3=-1 + call pack77(msg,i3,n3,c77) + mbits=0 + read(c77,'(77i1)') mbits(1:77) + +! Message is 77 bits long. Add a 0 bit to create a 78-bit message and pad with +! 12 zeros to create 90-bit mbit array for CRC calculation. + call get_q65crc12(mbits,ncrc1,ncrc2) + +! Now have message in bits 1:78 and CRC in bits 79:90. +! Group message bits into 15 6-bit symbols: + do i=0,14 + write(c6,'(6i1)') mbits( (i*6+1):(i*6+6) ) + read(c6,'(b6.6)') message(i+1) + enddo + +! Encode to create a 65-symbol codeword + call q65_encode(message,codeword) + + write(*,*) 'Generated message plus CRC (90 bits)' + write(*,'(a8,15i4)') '6 bit : ',message + write(*,'(a8,90i1)') 'binary: ',mbits + write(*,*) ' ' + write(*,*) 'Codeword with CRC symbols (65 symbols)' + write(*,'(20i3)') codeword + +!Shorten the codeword by omitting the CRC symbols (symbols 14 and 15) + shortcodeword(1:13)=codeword(1:13) + shortcodeword(14:63)=codeword(16:65) + +!Insert sync symbols to create array of channel symbols + j=1 + k=0 + do i=1,85 + if(i.eq.isync(j)) then + j=j+1 + itone(i)=0 + else + k=k+1 + itone(i)=shortcodeword(k)+1 + endif + enddo + + write(*,*) ' ' + write(*,*) 'Channel symbols (85 total)' + write(*,'(20i3)') itone + +999 end program q65code diff --git a/map65/libm65/CMakeLists.txt b/map65/libm65/CMakeLists.txt index 0409eb412..ab40d521a 100644 --- a/map65/libm65/CMakeLists.txt +++ b/map65/libm65/CMakeLists.txt @@ -37,6 +37,7 @@ set (libm65_FSRCS ftnquit.f90 q65b.f90 gen65.f90 + gen_q65_cwave.f90 gen_q65_wave.f90 geocentric.f90 getdphi.f90 @@ -94,7 +95,7 @@ set_source_files_properties (${libm65_ka9q_CSRCS} PROPERTIES COMPILE_FLAGS -Wno- set (libm65_CSRCS ${libm65_ka9q_CSRCS} ftrsd2.c - gran.c +# gran.c igray.c tmoonsub.c usleep.c diff --git a/map65/libm65/gen_q65_cwave.f90 b/map65/libm65/gen_q65_cwave.f90 new file mode 100644 index 000000000..f517917b5 --- /dev/null +++ b/map65/libm65/gen_q65_cwave.f90 @@ -0,0 +1,70 @@ +subroutine gen_q65_cwave(msg,ntxfreq,ntone_spacing,msgsent,cwave,nwave) + +! Encodes a Q65 message to yield complex cwave() at fsample = 96000 Hz + + use packjt + parameter (NMAX=60*96000) + character*22 msg + character*22 msgsent !Message as it will be received + character*120 cmnd + character*80 wsjtx_dir + character*16 cjunk + real*8 t,dt,phi,f,f0,dfgen,dphi,twopi,tsym + complex cwave(NMAX) + integer itone(85) + integer icos7(0:6) + data icos7/2,5,6,0,4,1,3/ !Defines a 7x7 Costas array + data twopi/6.283185307179586476d0/ + save + + open(9,file='wsjtx_dir.txt',status='old') + read(9,'(a)') wsjtx_dir + close(9) + + msgsent=msg +! 1 2 3 4 5 +! 12345678901234567890123456789012345678901234567890123456789012345 + cmnd='q65sim "K1ABC W9XYZ EN37 " A 1500 0 0 0 0 60 0 99 >itone.txt' + cmnd(9:30)=msg + write(cmnd(35:38),'(i4)') ntxfreq + cmnd=trim(wsjtx_dir)//cmnd + call execute_command_line(cmnd) + open(9,file='itone.txt',status='old') + do i=1,99 + read(9,1000,end=999) cjunk +1000 format(a16) + if(cjunk.eq.'Channel symbols:') exit + enddo + read(9,1002) itone +1002 format(20i3) + close(9) + +! Set up necessary constants + nsym=85 + tsym=7200.d0/12000.d0 + dt=1.d0/96000.d0 + f0=ntxfreq + dfgen=ntone_spacing*12000.d0/7200.d0 + phi=0.d0 + dphi=twopi*dt*f0 + i=0 + nwave=85*7200*96000.d0/12000.d0 + t=0.d0 + j0=0 + do i=1,nwave + t=t+dt + j=t/tsym + 1 + if(j.gt.85) exit + if(j.ne.j0) then + f=f0 + itone(j)*dfgen + dphi=twopi*dt*f + j0=j + endif + phi=phi+dphi + if(phi.gt.twopi) phi=phi-twopi + xphi=phi + cwave(i)=cmplx(cos(xphi),-sin(xphi)) + enddo + +999 return +end subroutine gen_q65_cwave diff --git a/map65/libm65/mapsim.f90 b/map65/libm65/mapsim.f90 index 5142584ab..5aeede8cc 100644 --- a/map65/libm65/mapsim.f90 +++ b/map65/libm65/mapsim.f90 @@ -2,7 +2,7 @@ program mapsim ! Generate simulated data for testing of MAP65 - parameter (NMAX=96000*60) + parameter (NMAX=60*96000) real*4 d4(4,NMAX) !Floating-point data integer*2 id4(4,NMAX) !i*2 data, dual polarization integer*2 id2(2,NMAX) !i*2 data, single polarization @@ -10,35 +10,41 @@ program mapsim complex z,zx,zy real*8 fcenter,fsample,samfac,f,dt,twopi,phi,dphi character msg0*22,message*22,msgsent*22,arg*8,fname*13,mode*2 + logical bq65 nargs=iargc() - if(nargs.ne.9) then - print*,'Usage: mapsim DT "message" mode f1 f2 nsigs pol SNR nfiles' - print*,'Example: 2.5 "CQ K1ABC FN42" B -22 33 20 45 -20 1' + if(nargs.ne.10) then + print*,'Usage: mapsim "message" mode DT fa fb nsigs pol fDop SNR nfiles' + print*,'Example: mapsim "CQ K1ABC FN42" B 2.5 -20 20 21 45 0.0 -20 1' print*,' ' - print*,'Enter message = "" to use entries in msgs.txt.' - print*,'Enter pol = -1 to generate a range of polarization angles.' - print*,'Enter SNR = 0 to generate a range of SNRs.' + print*,' mode = A B C for JT65; QA-QE for Q65-60A' + print*,' fa = lowest freq in kHz, relative to center' + print*,' fb = highest freq in kHz, relative to center' + print*,' message = "" to use entries in msgs.txt.' + print*,' pol = -1 to generate a range of polarization angles.' + print*,' SNR = 0 to generate a range of SNRs.' go to 999 endif - call getarg(1,arg) - read(arg,*) dt0 !Time delay - call getarg(2,msg0) + call getarg(1,msg0) message=msg0 !Transmitted message - call getarg(3,mode) !JT65 sub-mode (A B C B2 C2) + call getarg(2,mode) !JT65 sub-mode (A B C QA-QE) + call getarg(3,arg) + read(arg,*) dt0 !Time delay call getarg(4,arg) - read(arg,*) f1 !Lowest freq (kHz, relative to fcenter) + read(arg,*) fa !Lowest freq (kHz, relative to fcenter) call getarg(5,arg) - read(arg,*) f2 !Highest freq + read(arg,*) fb !Highest freq call getarg(6,arg) read(arg,*) nsigs !Number of signals in each file call getarg(7,arg) read(arg,*) npol !Polarization in degrees - call getarg(8,arg) - read(arg,*) snrdb !S/N pol=npol + call getarg(8,arg) + read(arg,*) fdop !Doppler spread call getarg(9,arg) + read(arg,*) snrdb !S/N + call getarg(10,arg) read(arg,*) nfiles !Number of files rmsdb=25. @@ -49,17 +55,20 @@ program mapsim twopi=8.d0*atan(1.d0) rad=360.0/twopi samfac=1.d0 - mode65=1 - if(mode(1:1).eq.'B') mode65=2 - if(mode(1:1).eq.'C') mode65=4 - nfast=1 - if(mode(2:2).eq.'2') nfast=2 - npts=NMAX/nfast + bq65=(mode(1:1).eq.'Q') + ntone_spacing=1 + ntxfreq=1270 + fac=1.0/32767.0 + if(mode(1:1).eq.'B' .or. mode(2:2).eq.'B') ntone_spacing=2 + if(mode(1:1).eq.'C' .or. mode(2:2).eq.'C') ntone_spacing=4 + if(mode(2:2).eq.'D') ntone_spacing=8 + if(mode(2:2).eq.'E') ntone_spacing=16 + npts=NMAX open(12,file='msgs.txt',status='old') write(*,1000) -1000 format('File N freq S/N pol Message'/ & - '---------------------------------------------------') +1000 format('File N Mode DT freq pol fDop SNR Message'/ & + '--------------------------------------------------------------------') do ifile=1,nfiles nmin=ifile-1 @@ -72,31 +81,43 @@ program mapsim call noisegen(d4,npts) !Generate Gaussuian noise if(msg0.ne.' ') then - call cgen65(message,mode65,nfast,samfac,nsendingsh,msgsent,cwave,nwave) + if(bq65) then + call gen_q65_cwave(message,ntxfreq,ntone_spacing,msgsent, & + cwave,nwave) + else + call cgen65(message,ntone_spacing,samfac,nsendingsh,msgsent, & + cwave,nwave) + endif endif - rewind 12 + + if(fdop.gt.0.0) call dopspread(cwave,nwave,fdop) + do isig=1,nsigs if(msg0.eq.' ') then read(12,1004) message 1004 format(a22) - call cgen65(message,mode65,nfast,samfac,nsendingsh,msgsent, & - cwave,nwave) + if(bq65) then + call gen_q65_cwave(msg,ntxfreq,mode65,msgsent,cwave,nwave) + else + call cgen65(message,ntone_spacing,samfac,nsendingsh,msgsent, & + cwave,nwave) + endif endif - + if(npol.lt.0) pol=(isig-1)*180.0/nsigs a=cos(pol/rad) b=sin(pol/rad) - f=1000.0*(f1+f2)/2.0 - if(nsigs.gt.1) f=1000.0*(f1 + (isig-1)*(f2-f1)/(nsigs-1.0)) + f=1000.0*(fa+fb)/2.0 + if(nsigs.gt.1) f=1000.0*(fa + (isig-1)*(fb-fa)/(nsigs-1.0)) dphi=twopi*f*dt + 0.5*twopi snrdbx=snrdb if(snrdb.ge.-1.0) snrdbx=-15.0 - 15.0*(isig-1.0)/nsigs sig=sqrt(2.2*2500.0/96000.0) * 10.0**(0.05*snrdbx) - write(*,1020) ifile,isig,0.001*f,snrdbx,nint(pol),msgsent -1020 format(i3,i4,f8.3,f7.1,i5,2x,a22) + write(*,1020) ifile,isig,mode,dt0,0.001*f,nint(pol),fDop,snrdbx,msgsent +1020 format(i3,i3,2x,a2,f6.2,f8.3,i5,2f7.1,2x,a22) phi=0. ! i0=fsample*(3.5d0+0.05d0*(isig-1)) @@ -133,3 +154,55 @@ program mapsim enddo 999 end program mapsim + +subroutine dopspread(cwave,nwave,fspread) + + parameter (NMAX=60*96000) + parameter (NFFT=NMAX,NH=NFFT/2) + complex cwave(NMAX) + complex cspread(0:NMAX-1) + + twopi=8.0*atan(1.0) + df=96000.0/nfft + cspread(0)=1.0 + cspread(NH)=0. + b=6.0 !Use truncated Lorenzian shape for fspread + do i=1,NH + f=i*df + x=b*f/fspread + z=0. + a=0. + if(x.lt.3.0) then !Cutoff beyond x=3 + a=sqrt(1.111/(1.0+x*x)-0.1) !Lorentzian amplitude + phi1=twopi*rran() !Random phase + z=a*cmplx(cos(phi1),sin(phi1)) + endif + cspread(i)=z + z=0. + if(x.lt.3.0) then !Same thing for negative freqs + phi2=twopi*rran() + z=a*cmplx(cos(phi2),sin(phi2)) + endif + cspread(nfft-i)=z + enddo + + call four2a(cspread,nfft,1,1,1) !Transform to time domain + + sum=0. + do i=0,nfft-1 + p=real(cspread(i))**2 + aimag(cspread(i))**2 + sum=sum+p + enddo + avep=sum/nfft + fac=sqrt(1.0/avep) + cspread=fac*cspread !Normalize to constant avg power + cwave=cspread*cwave !Apply Rayleigh fading + +! do i=0,nfft-1 +! p=real(cspread(i))**2 + aimag(cspread(i))**2 +! write(14,3010) i,p,cspread(i) +!3010 format(i8,3f12.6) +! enddo + + return +end subroutine dopspread