w4kek/WSJT-X
1
0
Fork 0
mirror of https://github.com/saitohirga/WSJT-X.git synced 2024-11-18 01:52:05 -05:00
Code Issues Packages Projects Releases Wiki Activity
0b7d77589e
WSJT-X/lib/jt65sim.f90
Bill Somerville dd1362b69a Rationalize NA contest mode 
Generic message packing and unpacking routines now understand antipode
grid contest messages.  These messages  are now recognized as standard
messages in  message response processing and  dealt with appropriately
when contest mode is selected and applicable (currently FT8 and MSK144
only).

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@8062 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
2017-09-01 20:10:35 +00:00

272 lines
9.8 KiB
Fortran
Raw Blame History

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(8) = [ &
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') ]
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
do
call getopt('hm:n:d:t:f:ps:',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
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)
msg="K1ABC W9XYZ EN37"
!###
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
Reference in New Issue View Git Blame Copy Permalink