From d48fb58ffa389083b2bd40224e5c065ac2519dda Mon Sep 17 00:00:00 2001 From: Bill Somerville Date: Mon, 1 Jan 2018 22:02:31 +0000 Subject: [PATCH] Enhance jt65sim to allow 11025 Hz rate, selectable base frequency and o/p gain offset. Needs a review. git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@8388 ab8295b8-cf94-4d9e-aec4-7959e3be5d79 --- lib/jt65sim.f90 | 568 +++++++++++++++++++++++++----------------------- 1 file changed, 293 insertions(+), 275 deletions(-) diff --git a/lib/jt65sim.f90 b/lib/jt65sim.f90 index b02a64ca9..eca533b47 100644 --- a/lib/jt65sim.f90 +++ b/lib/jt65sim.f90 @@ -1,275 +1,293 @@ -program jt65sim - -! Generate simulated JT65 data for testing WSJT-X - - use wavhdr - use packjt - use options - parameter (NMAX=54*12000) ! = 648,000 - parameter (NFFT=10*65536,NH=NFFT/2) - type(hdr) h !Header for .wav file - integer*2 iwave(NMAX) !Generated waveform - integer*4 itone(126) !Channel symbols (values 0-65) - integer dgen(12) !Twelve 6-bit data symbols - integer sent(63) !RS(63,12) codeword - real*4 xnoise(NMAX) !Generated random noise - real*4 dat(NMAX) !Generated real data - complex cdat(NMAX) !Generated complex waveform - complex cspread(0:NFFT-1) !Complex amplitude for Rayleigh fading - complex z - real*8 f0,dt,twopi,phi,dphi,baud,fsample,freq,sps - character msg*22,fname*11,csubmode*1,c,optarg*500,numbuf*32 -! character call1*5,call2*5 - logical :: display_help=.false.,seed_prngs=.true. - type (option) :: long_options(9) = [ & - option ('help',.false.,'h','Display this help message',''), & - option ('sub-mode',.true.,'m','sub mode, default MODE=A','MODE'), & - option ('num-sigs',.true.,'n','number of signals per file, default SIGNALS=10','SIGNALS'), & - option ('doppler-spread',.true.,'d','Doppler spread, default SPREAD=0.0','SPREAD'), & - option ('time-offset',.true.,'t','Time delta, default SECONDS=0.0','SECONDS'), & - option ('num-files',.true.,'f','Number of files to generate, default FILES=1','FILES'), & - option ('no-prng-seed',.false.,'p','Do not seed PRNGs (use for reproducible tests)',''), & - option ('strength',.true.,'s','S/N in dB (2500Hz reference b/w), default SNR=0','SNR'), & - option ('message',.true.,'M','Message text','Message') ] - - integer nprc(126) !Sync pattern - data nprc/1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0, & - 0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1, & - 0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1, & - 0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1, & - 1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1, & - 0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1, & - 1,1,1,1,1,1/ - -! Default parameters: - csubmode='A' - mode65=1 - nsigs=10 - fspread=0. - xdt=0. - snrdb=0. - nfiles=1 - msg="K1ABC W9XYZ EN37" - - do - call getopt('hm:n:d:t:f:ps:M:',long_options,c,optarg,narglen,nstat,noffset,nremain,.true.) - if( nstat .ne. 0 ) then - exit - end if - select case (c) - case ('h') - display_help = .true. - case ('m') - read (optarg(:narglen), *) csubmode - if(csubmode.eq.'A') mode65=1 - if(csubmode.eq.'B') mode65=2 - if(csubmode.eq.'C') mode65=4 - case ('n') - read (optarg(:narglen), *,err=10) nsigs - case ('d') - read (optarg(:narglen), *,err=10) fspread - case ('t') - read (optarg(:narglen), *) numbuf - if (numbuf(1:1) == '\') then - read (numbuf(2:), *,err=10) xdt - else - read (numbuf, *,err=10) xdt - end if - case ('f') - read (optarg(:narglen), *,err=10) nfiles - case ('p') - seed_prngs=.false. - case ('s') - read (optarg(:narglen), *) numbuf - if (numbuf(1:1) == '\') then - read (numbuf(2:), *,err=10) snrdb - else - read (numbuf, *,err=10) snrdb - end if - case ('M') - read (optarg(:narglen), '(A)',err=10) msg - write(*,*) msg - end select - cycle -10 display_help=.true. - print *, 'Optional argument format error for option -', c - end do - - if(display_help .or. nstat.lt.0 .or. nremain.ge.1) then - print *, '' - print *, 'Usage: jt65sim [OPTIONS]' - print *, '' - print *, ' Generate one or more simulated JT65 signals in .WAV file(s)' - print *, '' - print *, 'Example: jt65sim -m B -n 10 -d 0.2 -s \\-24.5 -t 0.0 -f 4' - print *, '' - print *, 'OPTIONS: NB Use \ (\\ on *nix shells) to escape -ve arguments' - print *, '' - do i = 1, size (long_options) - call long_options(i) % print (6) - end do - go to 999 - endif - - if (seed_prngs) then - call init_random_seed() ! seed Fortran RANDOM_NUMBER generator - call sgran() ! see C rand generator (used in gran) - end if - - rms=100. - fsample=12000.d0 !Sample rate (Hz) - dt=1.d0/fsample !Sample interval (s) - twopi=8.d0*atan(1.d0) - npts=54*12000 !Total samples in .wav file - baud=11025.d0/4096.d0 !Keying rate - sps=12000.d0/baud !Samples per symbol, at fsample=12000 Hz - nsym=126 !Number of channel symbols - h=default_header(12000,npts) - dfsig=2000.0/nsigs !Freq spacing between sigs in file (Hz) - - 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//'.wav',access='stream',status='unknown') - - xnoise=0. - cdat=0. - if(snrdb.lt.90) then - do i=1,npts - xnoise(i)=gran() !Generate gaussian noise - enddo - endif - - do isig=1,nsigs !Generate requested number of sigs - if(mod(nsigs,2).eq.0) f0=1500.0 + dfsig*(isig-0.5-nsigs/2) - if(mod(nsigs,2).eq.1) f0=1500.0 + dfsig*(isig-(nsigs+1)/2) - xsnr=snrdb - if(snrdb.eq.0.0) xsnr=-19 - isig - if(csubmode.eq.'B' .and. snrdb.eq.0.0) xsnr=-21 - isig - if(csubmode.eq.'C' .and. snrdb.eq.0.0) xsnr=-21 - isig - -!### -! call1="K1ABC" -! ic3=65+mod(isig-1,26) -! ic2=65+mod((isig-1)/26,26) -! ic1=65 -! call2="W9"//char(ic1)//char(ic2)//char(ic3) -! write(msg,1010) call1,call2,nint(xsnr) -!1010 format(a5,1x,a5,1x,i3.2) -!### - call packmsg(msg,dgen,itype,.false.) !Pack message into 12 six-bit bytes - call rs_encode(dgen,sent) !Encode using RS(63,12) - call interleave63(sent,1) !Interleave channel symbols - call graycode65(sent,63,1) !Apply Gray code - - k=0 - do j=1,nsym !Insert sync and data into itone() - if(nprc(j).eq.0) then - k=k+1 - itone(j)=sent(k)+2 - else - itone(j)=0 - endif - enddo - - bandwidth_ratio=2500.0/6000.0 - sig=sqrt(2*bandwidth_ratio)*10.0**(0.05*xsnr) - if(xsnr.gt.90.0) sig=1.0 - write(*,1020) ifile,isig,f0,csubmode,xsnr,xdt,fspread,msg -1020 format(i4,i4,f10.3,2x,a1,2x,f5.1,f6.2,f5.1,1x,a22) - - phi=0.d0 - dphi=0.d0 - k=12000 + xdt*12000 !Start audio at t = xdt + 1.0 s - isym0=-99 - do i=1,npts !Add this signal into cdat() - isym=floor(i/sps)+1 - if(isym.gt.nsym) exit - if(isym.ne.isym0) then - freq=f0 + itone(isym)*baud*mode65 - dphi=twopi*freq*dt - isym0=isym - endif - phi=phi + dphi - if(phi.gt.twopi) phi=phi-twopi - xphi=phi - z=cmplx(cos(xphi),sin(xphi)) - k=k+1 - if(k.ge.1) cdat(k)=cdat(k) + sig*z - enddo - enddo - - if(fspread.ne.0) then !Apply specified Doppler spread - df=12000.0/nfft - twopi=8*atan(1.0) - cspread(0)=1.0 - cspread(NH)=0. - -! The following options were added 3/15/2016 to make the half-power tone -! widths equal to the requested Doppler spread. (Previously we effectively -! used b=1.0 and Gaussian shape, which made the tones 1.665 times wider.) -! b=2.0*sqrt(log(2.0)) !Gaussian (before 3/15/2016) -! b=2.0 !Lorenzian 3/15 - 3/27 - b=6.0 !Lorenzian 3/28 onward - - 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(exp(-x*x)) !Gaussian - a=sqrt(1.111/(1.0+x*x)-0.1) !Lorentzian - call random_number(r1) - phi1=twopi*r1 - z=a*cmplx(cos(phi1),sin(phi1)) - endif - cspread(i)=z - z=0. - if(x.lt.50.0) then - call random_number(r2) - phi2=twopi*r2 - z=a*cmplx(cos(phi2),sin(phi2)) - endif - cspread(NFFT-i)=z - enddo - - do i=0,NFFT-1 - f=i*df - if(i.gt.NH) f=(i-nfft)*df - s=real(cspread(i))**2 + aimag(cspread(i))**2 -! write(13,3000) i,f,s,cspread(i) -!3000 format(i5,f10.3,3f12.6) - enddo -! s=real(cspread(0))**2 + aimag(cspread(0))**2 -! write(13,3000) 1024,0.0,s,cspread(0) - - 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 - cdat=cspread(1:npts)*cdat !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 - - endif - - dat=aimag(cdat) + xnoise !Add the generated noise - fac=32767.0/nsigs - if(snrdb.ge.90.0) iwave(1:npts)=nint(fac*dat(1:npts)) - if(snrdb.lt.90.0) iwave(1:npts)=nint(rms*dat(1:npts)) - write(10) h,iwave(1:npts) !Save the .wav file - close(10) - enddo - -999 end program jt65sim +program jt65sim + +! Generate simulated JT65 data for testing WSJT-X + + use wavhdr + use packjt + use options + parameter (NMAX=54*12000) ! = 648,000 @12kHz + parameter (NFFT=10*65536,NH=NFFT/2) + type(hdr) h !Header for .wav file + integer*2 iwave(NMAX) !Generated waveform + integer*4 itone(126) !Channel symbols (values 0-65) + integer dgen(12) !Twelve 6-bit data symbols + integer sent(63) !RS(63,12) codeword + real*4 xnoise(NMAX) !Generated random noise + real*4 dat(NMAX) !Generated real data + complex cdat(NMAX) !Generated complex waveform + complex cspread(0:NFFT-1) !Complex amplitude for Rayleigh fading + complex z + real*8 f0,dt,twopi,phi,dphi,baud,fsample,freq,sps + character msg*22,fname*11,csubmode*1,c,optarg*500,numbuf*32 +! character call1*5,call2*5 + logical :: display_help=.false.,seed_prngs=.true. + type (option) :: long_options(12) = [ & + option ('help',.false.,'h','Display this help message',''), & + option ('sub-mode',.true.,'m','sub mode, default MODE=A','MODE'), & + option ('num-sigs',.true.,'n','number of signals per file, default SIGNALS=10','SIGNALS'), & + option ('f0',.true.,'F','base frequency offset, default F0=1500.0','F0'), & + option ('doppler-spread',.true.,'d','Doppler spread, default SPREAD=0.0','SPREAD'), & + option ('time-offset',.true.,'t','Time delta, default SECONDS=0.0','SECONDS'), & + option ('num-files',.true.,'f','Number of files to generate, default FILES=1','FILES'), & + option ('no-prng-seed',.false.,'p','Do not seed PRNGs (use for reproducible tests)',''), & + option ('strength',.true.,'s','S/N in dB (2500Hz reference b/w), default SNR=0','SNR'), & + option ('11025',.false.,'S','Generate at 11025Hz sample rate, default 12000Hz',''), & + option ('gain-offset',.true.,'G','Gain offset in dB, default GAIN=0dB','GAIN'), & + option ('message',.true.,'M','Message text','Message') ] + + integer nprc(126) !Sync pattern + data nprc/1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0, & + 0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1, & + 0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1, & + 0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1, & + 1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1, & + 0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1, & + 1,1,1,1,1,1/ + +! Default parameters: + csubmode='A' + mode65=1 + nsigs=10 + bf0=1500. + fspread=0. + xdt=0. + snrdb=0. + nfiles=1 + nsample_rate=12000 + gain_offset=0. + msg="K1ABC W9XYZ EN37" + + do + call getopt('hm:n:F:d:t:f:ps:SG:M:',long_options,c,optarg,narglen,nstat,noffset,nremain,.true.) + if( nstat .ne. 0 ) then + exit + end if + select case (c) + case ('h') + display_help = .true. + case ('m') + read (optarg(:narglen), *) csubmode + if(csubmode.eq.'A') mode65=1 + if(csubmode.eq.'B') mode65=2 + if(csubmode.eq.'C') mode65=4 + case ('n') + read (optarg(:narglen), *,err=10) nsigs + case ('F') + read (optarg(:narglen), *,err=10) bf0 + case ('d') + read (optarg(:narglen), *,err=10) fspread + case ('t') + read (optarg(:narglen), *) numbuf + if (numbuf(1:1) == '\') then !'\' + read (numbuf(2:), *,err=10) xdt + else + read (numbuf, *,err=10) xdt + end if + case ('f') + read (optarg(:narglen), *,err=10) nfiles + case ('p') + seed_prngs=.false. + case ('s') + read (optarg(:narglen), *) numbuf + if (numbuf(1:1) == '\') then !'\' + read (numbuf(2:), *,err=10) snrdb + else + read (numbuf, *,err=10) snrdb + end if + case ('S') + nsample_rate=11025 + case ('G') + read (optarg(:narglen), *) numbuf + if (numbuf(1:1) == '\') then !'\' + read (numbuf(2:), *, err=10) gain_offset + else + read (numbuf, *, err=10) gain_offset + end if + case ('M') + read (optarg(:narglen), '(A)',err=10) msg + write(*,*) msg + end select + cycle +10 display_help=.true. + print *, 'Optional argument format error for option -', c + end do + + if(display_help .or. nstat.lt.0 .or. nremain.ge.1) then + print *, '' + print *, 'Usage: jt65sim [OPTIONS]' + print *, '' + print *, ' Generate one or more simulated JT65 signals in .WAV file(s)' + print *, '' + print *, 'Example: jt65sim -m B -n 10 -d 0.2 -s \\-24.5 -t 0.0 -f 4' + print *, '' + print *, 'OPTIONS: NB Use \ (\\ on *nix shells) to escape -ve arguments' + print *, '' + do i = 1, size (long_options) + call long_options(i) % print (6) + end do + go to 999 + endif + + if (seed_prngs) then + call init_random_seed() ! seed Fortran RANDOM_NUMBER generator + call sgran() ! see C rand generator (used in gran) + end if + + rms=100. * 10. ** (gain_offset / 20.) + + fsample=nsample_rate !Sample rate (Hz) + dt=1.d0/fsample !Sample interval (s) + twopi=8.d0*atan(1.d0) + npts=54*nsample_rate !Total samples in .wav file + baud=11025.d0/4096.d0 !Keying rate + sps=real(nsample_rate)/baud !Samples per symbol, at fsample=NSAMPLE_RATE Hz + nsym=126 !Number of channel symbols + h=default_header(nsample_rate,npts) + dfsig=2000.0/nsigs !Freq spacing between sigs in file (Hz) + + 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//'.wav',access='stream',status='unknown') + + xnoise=0. + cdat=0. + if(snrdb.lt.90) then + do i=1,npts + xnoise(i)=gran() !Generate gaussian noise + enddo + endif + + do isig=1,nsigs !Generate requested number of sigs + if(mod(nsigs,2).eq.0) f0=bf0 + dfsig*(isig-0.5-nsigs/2) + if(mod(nsigs,2).eq.1) f0=bf0 + dfsig*(isig-(nsigs+1)/2) + xsnr=snrdb + if(snrdb.eq.0.0) xsnr=-19 - isig + if(csubmode.eq.'B' .and. snrdb.eq.0.0) xsnr=-21 - isig + if(csubmode.eq.'C' .and. snrdb.eq.0.0) xsnr=-21 - isig + +!### +! call1="K1ABC" +! ic3=65+mod(isig-1,26) +! ic2=65+mod((isig-1)/26,26) +! ic1=65 +! call2="W9"//char(ic1)//char(ic2)//char(ic3) +! write(msg,1010) call1,call2,nint(xsnr) +!1010 format(a5,1x,a5,1x,i3.2) +!### + call packmsg(msg,dgen,itype,.false.) !Pack message into 12 six-bit bytes + call rs_encode(dgen,sent) !Encode using RS(63,12) + call interleave63(sent,1) !Interleave channel symbols + call graycode65(sent,63,1) !Apply Gray code + + k=0 + do j=1,nsym !Insert sync and data into itone() + if(nprc(j).eq.0) then + k=k+1 + itone(j)=sent(k)+2 + else + itone(j)=0 + endif + enddo + + bandwidth_ratio=2500.0/6000.0 + sig=sqrt(2*bandwidth_ratio)*10.0**(0.05*xsnr) + if(xsnr.gt.90.0) sig=1.0 + write(*,1020) ifile,isig,f0,csubmode,xsnr,xdt,fspread,msg +1020 format(i4,i4,f10.3,2x,a1,2x,f5.1,f6.2,f5.1,1x,a22) + + phi=0.d0 + dphi=0.d0 + k=nsample_rate + xdt*nsample_rate !Start audio at t = xdt + 1.0 s + isym0=-99 + do i=1,npts !Add this signal into cdat() + isym=floor(i/sps)+1 + if(isym.gt.nsym) exit + if(isym.ne.isym0) then + freq=f0 + itone(isym)*baud*mode65 + dphi=twopi*freq*dt + isym0=isym + endif + phi=phi + dphi + if(phi.gt.twopi) phi=phi-twopi + xphi=phi + z=cmplx(cos(xphi),sin(xphi)) + k=k+1 + if(k.ge.1) cdat(k)=cdat(k) + sig*z + enddo + enddo + + if(fspread.ne.0) then !Apply specified Doppler spread + df=real(nsample_rate)/nfft + twopi=8*atan(1.0) + cspread(0)=1.0 + cspread(NH)=0. + +! The following options were added 3/15/2016 to make the half-power tone +! widths equal to the requested Doppler spread. (Previously we effectively +! used b=1.0 and Gaussian shape, which made the tones 1.665 times wider.) +! b=2.0*sqrt(log(2.0)) !Gaussian (before 3/15/2016) +! b=2.0 !Lorenzian 3/15 - 3/27 + b=6.0 !Lorenzian 3/28 onward + + 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(exp(-x*x)) !Gaussian + a=sqrt(1.111/(1.0+x*x)-0.1) !Lorentzian + call random_number(r1) + phi1=twopi*r1 + z=a*cmplx(cos(phi1),sin(phi1)) + endif + cspread(i)=z + z=0. + if(x.lt.50.0) then + call random_number(r2) + phi2=twopi*r2 + z=a*cmplx(cos(phi2),sin(phi2)) + endif + cspread(NFFT-i)=z + enddo + + do i=0,NFFT-1 + f=i*df + if(i.gt.NH) f=(i-nfft)*df + s=real(cspread(i))**2 + aimag(cspread(i))**2 +! write(13,3000) i,f,s,cspread(i) +!3000 format(i5,f10.3,3f12.6) + enddo +! s=real(cspread(0))**2 + aimag(cspread(0))**2 +! write(13,3000) 1024,0.0,s,cspread(0) + + 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 + cdat(1:npts)=cspread(1:npts)*cdat(1:npts) !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 + + endif + + dat=aimag(cdat) + xnoise !Add the generated noise + fac=32767.0/nsigs + if(snrdb.ge.90.0) iwave(1:npts)=nint(fac*dat(1:npts)) + if(snrdb.lt.90.0) iwave(1:npts)=nint(rms*dat(1:npts)) + write(10) h,iwave(1:npts) !Save the .wav file + close(10) + enddo + +999 end program jt65sim