program q65sim ! Generate simulated Q65 data for testing the decoder. use wavhdr use packjt parameter (NMAX=300*12000) !Total samples in .wav file type(hdr) h !Header for .wav file integer*2 iwave(NMAX) !Generated waveform integer itone(85) !Channel symbols (values 0-65) integer ntone(85,10) !Channel symbols for up to 10 messages integer y(63) !Codeword integer istart !averaging compatible start seconds integer imins !minutes for 15s period timestamp integer isecs !seconds for 15s period timestamp real*4 xnoise(NMAX) !Generated random noise real*4 dat(NMAX) !Generated real data complex cdat(NMAX) !Generated complex waveform complex cspread(0:NMAX-1) !Complex amplitude for Rayleigh fading complex z real*8 f00,f0,dt,twopi,phi,dphi,baud,fsample,freq character fname*17,csubmode*1,arg*12,c2*2 character*37 msg,msgsent,imsg(10) nargs=iargc() if(nargs.ne.11) then print*,'Usage: q65sim "msg" A-E freq fDop DT f1 Stp TRp Nsig Nfile SNR' print*,'Example: q65sim "K1ABC W9XYZ EN37" A 1500 0.0 0.0 0.0 1 60 1 1 -26' print*,'Example: q65sim "ST" A 1500 0.0 0.0 0.0 1 60 1 -26' print*,' fDop = Doppler spread' print*,' f1 = Drift or Doppler rate (Hz/min)' print*,' Stp = Step size (Hz)' print*,' Stp = 0 implies no Doppler tracking' print*,' Nsig = number of generated signals, 1 - 10' print*,' Creates filenames which increment to permit averaging in first period' print*,' If msg = ST program produces a single tone at freq' go to 999 endif call getarg(1,msg) call getarg(2,csubmode) mode65=2**(ichar(csubmode)-ichar('A')) call getarg(3,arg) read(arg,*) f00 call getarg(4,arg) read(arg,*) fspread call getarg(5,arg) read(arg,*) xdt call getarg(6,arg) read(arg,*) f1 call getarg(7,arg) read(arg,*) nstp call getarg(8,arg) read(arg,*) ntrperiod call getarg(9,arg) read(arg,*) nsig call getarg(10,arg) read(arg,*) nfiles call getarg(11,arg) read(arg,*) snrdb if(ntrperiod.eq.15) then nsps=1800 else if(ntrperiod.eq.30) then nsps=3600 else if(ntrperiod.eq.60) then nsps=7200 else if(ntrperiod.eq.120) then nsps=16000 else if(ntrperiod.eq.300) then nsps=41472 else print*,'Invalid TR period' go to 999 endif rms=100. fsample=12000.d0 !Sample rate (Hz) npts=fsample*ntrperiod !Total samples in .wav file nfft=npts nh=nfft/2 dt=1.d0/fsample !Sample interval (s) twopi=8.d0*atan(1.d0) nsym=85 !Number of channel symbols mode65=2**(ichar(csubmode) - ichar('A')) imsg(1)=msg if(nsig.ge.2) then i0=index(msg,' ') i0=i0 + index(msg(i0+1:),' ')-2 do i=1,nsig c2=char(ichar('A')+i-1)//char(ichar('A')+i-1) imsg(i)=msg(1:i0-1)//c2//msg(i0+2:) enddo endif ichk=0 do i=1,nsig msg=imsg(i) call genq65(msg,ichk,msgsent,itone,i3,n3) ntone(:,i)=itone enddo if(nsig.eq.1) then j=0 do i=1,85 if(itone(i).gt.0) then j=j+1 y(j)=itone(i)-1 endif enddo write(*,1001) y(1:13),y(1:13) 1001 format('Generated message'/'6-bit: ',13i3/'binary: ',13b6.6) write(*,1002) y 1002 format(/'Codeword:'/(20i3)) write(*,1003) itone 1003 format(/'Channel symbols:'/(20i3)) endif baud=12000.d0/nsps !Keying rate (6.67 baud fot 15-s sequences) h=default_header(12000,npts) write(*,1004) 1004 format('File TR Freq Mode S/N Dop DT f1 Stp Message'/70('-')) do ifile=1,nfiles !Loop over requested number of files istart = (ifile*ntrperiod*2) - (ntrperiod*2) if(ntrperiod.lt.30) then !wdg was 60 imins=istart/60 isecs=istart-(60*imins) write(fname,1005) imins,isecs !Construction of output filename for 15s periods with averaging 1005 format('000000_',i4.4, i2.2,'.wav') else write(fname,1106) istart/60 !Output filename to be compatible with averaging 30-300s periods 1106 format('000000_',i4.4,'.wav') endif open(10,file=trim(fname),access='stream',status='unknown') xnoise=0. if(snrdb.lt.90) then do i=1,npts xnoise(i)=gran() !Generate gaussian noise enddo endif cdat=0. bandwidth_ratio=2500.0/6000.0 sig=sqrt(2*bandwidth_ratio)*10.0**(0.05*snrdb) if(snrdb.gt.90.0) sig=1.0 write(*,1020) ifile,ntrperiod,f00,csubmode,snrdb,fspread,xdt,f1,nstp,trim(msgsent) 1020 format(i4,i6,f7.1,2x,a1,2x,f5.1,1x,f6.2,2f6.1,i4,2x,a) n=65.0*baud*mode65/100.0 + 0.9999 nfstep=100*n nf1=1500 - nfstep*(nsig-1)/2 do n=1,nsig f0=f00 if(nsig.ge.2) then f0=nf1 + (n-1)*nfstep itone=ntone(:,n) endif phi=0.d0 dphi=0.d0 k=(xdt+0.5)*12000 !Start audio at t=xdt+0.5 s (TR=15 and 30 s) if(ntrperiod.ge.60) k=(xdt+1.0)*12000 !TR 60+ at t = xdt + 1.0 s isym0=-99 do i=1,npts !Add this signal into cdat() isym=i/nsps + 1 if(isym.gt.nsym) exit if(isym.ne.isym0) then freq_drift=f1*i*dt/60.0 if(nstp.ne.0) freq_drift=freq_drift - nstp*nint(freq_drift/nstp) if (msg(1:2).eq.'ST') then freq = f0 + freq_drift else freq = f0 + freq_drift + itone(isym)*baud*mode65 endif 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 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 cdat=cspread*cdat !Apply Rayleigh fading endif ! psig=0. ! pnoise=0. ! do i=1,npts ! if(i.gt.12000 .and. i.lt.624000) then ! psig=psig + aimag(cdat(i))**2 ! pnoise=pnoise + xnoise(i)*xnoise(i) ! endif ! enddo ! pnoise=pnoise*bandwidth_ratio ! snr_2500=db(psig/pnoise) !Calibration confirmation! ! print*,'SNR_2500:',snr_2500 dat=aimag(cdat) + xnoise !Add generated AWGN noise fac=32767.0 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 q65sim