From c5e2593979cf0b51e85e6670eb7c33748d74d94c Mon Sep 17 00:00:00 2001 From: Steven Franke Date: Tue, 14 Apr 2020 10:34:00 -0500 Subject: [PATCH] Add routines necessary to generate simulated 'wspr4' wav files. --- CMakeLists.txt | 6 ++ lib/fsk4hf/gen_wspr4wave.f90 | 68 +++++++++++++++++++++ lib/fsk4hf/genwspr4.f90 | 97 +++++++++++++++++++++++++++++ lib/fsk4hf/wspr4_params.f90 | 16 +++++ lib/fsk4hf/wspr4sim.f90 | 115 +++++++++++++++++++++++++++++++++++ 5 files changed, 302 insertions(+) create mode 100644 lib/fsk4hf/gen_wspr4wave.f90 create mode 100644 lib/fsk4hf/genwspr4.f90 create mode 100644 lib/fsk4hf/wspr4_params.f90 create mode 100644 lib/fsk4hf/wspr4sim.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 8458c6937..a420ad23b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -610,6 +610,9 @@ set (wsjt_FSRCS lib/fsk4hf/encode174_74.f90 lib/fsk4hf/bpdecode174_74.f90 lib/fsk4hf/osd174_74.f90 + lib/fsk4hf/wspr4sim.f90 + lib/fsk4hf/genwspr4.f90 + lib/fsk4hf/gen_wspr4wave.f90 ) # temporary workaround for a gfortran v7.3 ICE on Fedora 27 64-bit @@ -1357,6 +1360,9 @@ target_link_libraries (ldpcsim174_91 wsjt_fort wsjt_cxx) add_executable (ldpcsim174_74 lib/fsk4hf/ldpcsim174_74.f90 wsjtx.rc) target_link_libraries (ldpcsim174_74 wsjt_fort wsjt_cxx) +add_executable (wspr4sim lib/fsk4hf/wspr4sim.f90 wsjtx.rc) +target_link_libraries (wspr4sim wsjt_fort wsjt_cxx) + endif(WSJT_BUILD_UTILS) # build the main application diff --git a/lib/fsk4hf/gen_wspr4wave.f90 b/lib/fsk4hf/gen_wspr4wave.f90 new file mode 100644 index 000000000..b19f57bf3 --- /dev/null +++ b/lib/fsk4hf/gen_wspr4wave.f90 @@ -0,0 +1,68 @@ +subroutine gen_wspr4wave(itone,nsym,nsps,fsample,f0,cwave,wave,icmplx,nwave) + + real wave(nwave) + complex cwave(nwave) + real, allocatable, save :: pulse(:) + real, allocatable :: dphi(:) + integer itone(nsym) + logical first + data first/.true./ + save pulse,first,twopi,dt,hmod + + if(first) then + allocate( pulse(3*nsps*fsample) ) + twopi=8.0*atan(1.0) + dt=1.0/fsample + hmod=1.0 +! Compute the smoothed frequency-deviation pulse + do i=1,3*nsps + tt=(i-1.5*nsps)/real(nsps) + pulse(i)=gfsk_pulse(1.0,tt) + enddo + first=.false. + endif + +! Compute the smoothed frequency waveform. +! Length = (nsym+2)*nsps samples, zero-padded + allocate( dphi(0:(nsym+2)*nsps-1) ) + dphi_peak=twopi*hmod/real(nsps) + dphi=0.0 + do j=1,nsym + ib=(j-1)*nsps + ie=ib+3*nsps-1 + dphi(ib:ie) = dphi(ib:ie) + dphi_peak*pulse(1:3*nsps)*itone(j) + enddo + +! Calculate and insert the audio waveform + phi=0.0 + dphi = dphi + twopi*f0*dt !Shift frequency up by f0 + wave=0. + if(icmplx.eq.1) cwave=0. + k=0 + do j=0,(nsym+2)*nsps-1 + k=k+1 + if(icmplx.eq.0) then + wave(k)=sin(phi) + else + cwave(k)=cmplx(cos(phi),sin(phi)) + endif + phi=mod(phi+dphi(j),twopi) + enddo + +! Compute the ramp-up and ramp-down symbols + if(icmplx.eq.0) then + wave(1:nsps)=wave(1:nsps) * & + (1.0-cos(twopi*(/(i,i=0,nsps-1)/)/(2.0*nsps)))/2.0 + k1=(nsym+1)*nsps+1 + wave(k1:k1+nsps-1)=wave(k1:k1+nsps-1) * & + (1.0+cos(twopi*(/(i,i=0,nsps-1)/)/(2.0*nsps)))/2.0 + else + cwave(1:nsps)=cwave(1:nsps) * & + (1.0-cos(twopi*(/(i,i=0,nsps-1)/)/(2.0*nsps)))/2.0 + k1=(nsym+1)*nsps+1 + cwave(k1:k1+nsps-1)=cwave(k1:k1+nsps-1) * & + (1.0+cos(twopi*(/(i,i=0,nsps-1)/)/(2.0*nsps)))/2.0 + endif + + return +end subroutine gen_wspr4wave diff --git a/lib/fsk4hf/genwspr4.f90 b/lib/fsk4hf/genwspr4.f90 new file mode 100644 index 000000000..d3747f41a --- /dev/null +++ b/lib/fsk4hf/genwspr4.f90 @@ -0,0 +1,97 @@ +subroutine genwspr4(msg0,ichk,msgsent,msgbits,i4tone) + +! Encode an FT4 message +! Input: +! - msg0 requested message to be transmitted +! - ichk if ichk=1, return only msgsent +! - msgsent message as it will be decoded +! - i4tone array of audio tone values, {0,1,2,3} + +! Frame structure: +! s16 + 87symbols + 2 ramp up/down = 105 total channel symbols +! r1 + s4 + d29 + s4 + d29 + s4 + d29 + s4 + r1 + +! Message duration: TxT = 105*13312/12000 = 116.48 s + +! use iso_c_binding, only: c_loc,c_size_t + + use packjt77 + include 'wspr4_params.f90' + character*37 msg0 + character*37 message !Message to be generated + character*37 msgsent !Message as it will be received + character*77 c77 + character*24 c24 + integer*4 i4tone(NN),itmp(ND) + integer*1 codeword(2*ND) + integer*1 msgbits(74),rvec(77) + integer icos4a(4),icos4b(4),icos4c(4),icos4d(4) + integer*4 ncrc24 + logical unpk77_success + data icos4a/0,1,3,2/ + data icos4b/1,0,2,3/ + data icos4c/2,3,1,0/ + data icos4d/3,2,0,1/ + data rvec/0,1,0,0,1,0,1,0,0,1,0,1,1,1,1,0,1,0,0,0,1,0,0,1,1,0,1,1,0, & + 1,0,0,1,0,1,1,0,0,0,0,1,0,0,0,1,0,1,0,0,1,1,1,1,0,0,1,0,1, & + 0,1,0,1,0,1,1,0,1,1,1,1,1,0,0,0,1,0,1/ + message=msg0 + + do i=1, 37 + if(ichar(message(i:i)).eq.0) then + message(i:37)=' ' + exit + endif + enddo + do i=1,37 !Strip leading blanks + if(message(1:1).ne.' ') exit + message=message(i+1:) + enddo + + i3=-1 + n3=-1 + call pack77(message,i3,n3,c77) + call unpack77(c77,0,msgsent,unpk77_success) !Unpack to get msgsent + msgbits=0 + read(c77,'(50i1)') msgbits(1:50) + call get_crc24(msgbits,ncrc24) + write(c24,'(b24.24)') ncrc24 + read(c24,'(24i1)') msgbits(51:74) + + if(ichk.eq.1) go to 999 + read(c77,'(74i1)',err=1) msgbits + if(unpk77_success) go to 2 +1 msgbits=0 + itone=0 + msgsent='*** bad message *** ' + go to 999 + +entry get_wspr4_tones_from_74bits(msgbits,i4tone) + +2 msgbits=mod(msgbits+rvec(1:74),2) + call encode174_74(msgbits,codeword) + +! Grayscale mapping: +! bits tone +! 00 0 +! 01 1 +! 11 2 +! 10 3 + + do i=1,ND + is=codeword(2*i)+2*codeword(2*i-1) + if(is.le.1) itmp(i)=is + if(is.eq.2) itmp(i)=3 + if(is.eq.3) itmp(i)=2 + enddo + + i4tone(1:4)=icos4a + i4tone(5:33)=itmp(1:29) + i4tone(34:37)=icos4b + i4tone(38:66)=itmp(30:58) + i4tone(67:70)=icos4c + i4tone(71:99)=itmp(59:87) + i4tone(100:103)=icos4d + +999 return +end subroutine genwspr4 diff --git a/lib/fsk4hf/wspr4_params.f90 b/lib/fsk4hf/wspr4_params.f90 new file mode 100644 index 000000000..d4e575382 --- /dev/null +++ b/lib/fsk4hf/wspr4_params.f90 @@ -0,0 +1,16 @@ +! WSPR4 +! LDPC(174,74)/CRC24 code, four 4x4 Costas arrays for sync, ramp-up and ramp-down symbols + +parameter (KK=50) !Information bits (50 + CRC24) +parameter (ND=87) !Data symbols +parameter (NS=16) !Sync symbols +parameter (NN=NS+ND) !Sync and data symbols (103) +parameter (NN2=NS+ND+2) !Total channel symbols (105) +parameter (NSPS=13312) !Samples per symbol at 12000 S/s +parameter (NZ=NSPS*NN) !Sync and Data samples (1,397,136) +parameter (NZ2=NSPS*NN2) !Total samples in shaped waveform (1,397,760) +parameter (NMAX=408*3456) !Samples in iwave (1,410,048) +parameter (NFFT1=4*NSPS, NH1=NFFT1/2) !Length of FFTs for symbol spectra +parameter (NSTEP=NSPS) !Coarse time-sync step size +parameter (NHSYM=(NMAX-NFFT1)/NSTEP) !Number of symbol spectra (1/4-sym steps) +parameter (NDOWN=32) !Downsample factor diff --git a/lib/fsk4hf/wspr4sim.f90 b/lib/fsk4hf/wspr4sim.f90 new file mode 100644 index 000000000..ecb359b2f --- /dev/null +++ b/lib/fsk4hf/wspr4sim.f90 @@ -0,0 +1,115 @@ +program wspr4sim + +! Generate simulated signals for experimental "wspr4" mode + + use wavhdr + use packjt77 + include 'wspr4_params.f90' !Set various constants + type(hdr) h !Header for .wav file + character arg*12,fname*17 + character msg37*37,msgsent37*37 + character c77*77 + complex c0(0:NMAX-1) + complex c(0:NMAX-1) + real wave(NMAX) + integer itone(NN) + integer*1 msgbits(74) + integer*2 iwave(NMAX) !Generated full-length waveform + +! Get command-line argument(s) + nargs=iargc() + if(nargs.ne.7) then + print*,'Usage: wspr4sim "message" f0 DT fdop del nfiles snr' + print*,'Examples: wspr4sim "K9AN EN50 37" 1500 0.0 0.1 1.0 10 -15' + go to 999 + endif + call getarg(1,msg37) !Message to be transmitted + call getarg(2,arg) + read(arg,*) f0 !Frequency (only used for single-signal) + call getarg(3,arg) + read(arg,*) xdt !Time offset from nominal (s) + call getarg(4,arg) + read(arg,*) fspread !Watterson frequency spread (Hz) + call getarg(5,arg) + read(arg,*) delay !Watterson delay (ms) + call getarg(6,arg) + read(arg,*) nfiles !Number of files + call getarg(7,arg) + read(arg,*) snrdb !SNR_2500 + + nfiles=abs(nfiles) + twopi=8.0*atan(1.0) + fs=12000.0 !Sample rate (Hz) + dt=1.0/fs !Sample interval (s) + hmod=1.0 !Modulation index (0.5 is MSK, 1.0 is FSK) + tt=NSPS*dt !Duration of symbols (s) + baud=1.0/tt !Keying rate (baud) + txt=NZ2*dt !Transmission length (s) + + bandwidth_ratio=2500.0/(fs/2.0) + sig=sqrt(2*bandwidth_ratio) * 10.0**(0.05*snrdb) + if(snrdb.gt.90.0) sig=1.0 + + call genwspr4(msg37,0,msgsent37,msgbits,itone) + write(*,*) + write(*,'(a9,a37,3x,a7,i1,a1,i1)') 'Message: ',msgsent37,'i3.n3: ',i3,'.',n3 + write(*,1000) f0,xdt,txt,snrdb +1000 format('f0:',f9.3,' DT:',f6.2,' TxT:',f6.1,' SNR:',f6.1) + write(*,*) + if(i3.eq.1) then + write(*,*) ' mycall hiscall hisgrid' + write(*,'(28i1,1x,i1,1x,28i1,1x,i1,1x,i1,1x,15i1,1x,3i1)') msgbits(1:77) + else + write(*,'(a14)') 'Message bits: ' + write(*,'(77i1)') msgbits + endif + write(*,*) + write(*,'(a17)') 'Channel symbols: ' + write(*,'(76i1)') itone + write(*,*) + + call sgran() + + fsample=12000.0 + icmplx=1 + call gen_wspr4wave(itone,NN,NSPS,fsample,f0,c0,wave,icmplx,NMAX) + + k=nint((xdt+0.5)/dt)-NSPS + c0=cshift(c0,-k) + if(k.gt.0) c0(0:k-1)=0.0 + if(k.lt.0) c0(NMAX+k:NMAX-1)=0.0 + + do ifile=1,nfiles + c=c0 + if(fspread.ne.0.0 .or. delay.ne.0.0) call watterson(c,NMAX,NZ,fs,delay,fspread) + c=sig*c + wave=real(c) + peak=maxval(abs(wave)) + nslots=1 + if(snrdb.lt.90) then + do i=1,NMAX !Add gaussian noise at specified SNR + xnoise=gran() + wave(i)=wave(i) + xnoise + enddo + endif + gain=100.0 + if(snrdb.lt.90.0) then + wave=gain*wave + else + datpk=maxval(abs(wave)) + fac=32766.9/datpk + wave=fac*wave + endif + if(any(abs(wave).gt.32767.0)) print*,"Warning - data will be clipped." + iwave=nint(wave) + h=default_header(12000,NMAX) + write(fname,1102) ifile +1102 format('000000_',i6.6,'.wav') + open(10,file=fname,status='unknown',access='stream') + write(10) h,iwave !Save to *.wav file + close(10) + write(*,1110) ifile,xdt,f0,snrdb,fname +1110 format(i4,f7.2,f8.2,f7.1,2x,a17) + enddo + +999 end program wspr4sim