1program fst4sim
2
3! Generate simulated signals for experimental slow FT4  mode
4
5   use wavhdr
6   use packjt77
7   include 'fst4_params.f90'               !Set various constants
8   type(hdr) h                                !Header for .wav file
9   logical*1 wspr_hint
10   character arg*12,fname*17
11   character msg37*37,msgsent37*37,c77*77
12   complex, allocatable :: c0(:)
13   complex, allocatable :: c(:)
14   real, allocatable :: wave(:)
15   integer hmod
16   integer itone(NN)
17   integer*1 msgbits(101)
18   integer*2, allocatable :: iwave(:)        !Generated full-length waveform
19
20! Get command-line argument(s)
21   nargs=iargc()
22   if(nargs.ne.9) then
23      print*,'Need 9 arguments, got ',nargs
24      print*,'Usage:    fst4sim "message"        TRsec f0   DT  fdop del nfiles snr W'
25      print*,'Examples: fst4sim "K1JT K9AN EN50"  60  1500 0.0   0.1 1.0   10   -15 F'
26      print*,'W (T or F) argument is hint to encoder to use WSPR message when there is abiguity'
27      go to 999
28   endif
29   call getarg(1,msg37)                   !Message to be transmitted
30   call getarg(2,arg)
31   read(arg,*) nsec                       !TR sequence length, seconds
32   call getarg(3,arg)
33   read(arg,*) f00                        !Frequency (only used for single-signal)
34   call getarg(4,arg)
35   read(arg,*) xdt                        !Time offset from nominal (s)
36   call getarg(5,arg)
37   read(arg,*) fspread                    !Watterson frequency spread (Hz)
38   call getarg(6,arg)
39   read(arg,*) delay                      !Watterson delay (ms)
40   call getarg(7,arg)
41   read(arg,*) nfiles                     !Number of files
42   call getarg(8,arg)
43   read(arg,*) snrdb                      !SNR_2500
44   call getarg(9,arg)
45   read(arg,*) wspr_hint                  !0:break ties as 77-bit 1:break ties as 50-bit
46
47   nfiles=abs(nfiles)
48   twopi=8.0*atan(1.0)
49   fs=12000.0                             !Sample rate (Hz)
50   dt=1.0/fs                              !Sample interval (s)
51   nsps=0
52   if(nsec.eq.15) nsps=720
53   if(nsec.eq.30) nsps=1680
54   if(nsec.eq.60) nsps=3888
55   if(nsec.eq.120) nsps=8200
56   if(nsec.eq.300) nsps=21504
57   if(nsec.eq.900) nsps=66560
58   if(nsec.eq.1800) nsps=134400
59   if(nsps.eq.0) then
60      print*,'Invalid TR sequence length.'
61      go to 999
62   endif
63   baud=12000.0/nsps                      !Keying rate (baud)
64   nmax=nsec*12000
65   nz=nsps*NN
66   txt=nz*dt                              !Transmission length (s)
67   tt=nsps*dt                             !Duration of symbols (s)
68   nwave=max(nmax,(NN+2)*nsps)
69   allocate( c0(0:nwave-1) )
70   allocate( c(0:nwave-1) )
71   allocate( wave(nwave) )
72   allocate( iwave(nmax) )
73
74   bandwidth_ratio=2500.0/(fs/2.0)
75   sig=sqrt(2*bandwidth_ratio) * 10.0**(0.05*snrdb)
76   if(snrdb.gt.90.0) sig=1.0
77
78   if(wspr_hint) then
79      i3=0
80      n3=6
81   else
82      i3=-1
83      n3=-1
84   endif
85   call pack77(msg37,i3,n3,c77)
86   if(i3.eq.0.and.n3.eq.6) iwspr=1
87   call genfst4(msg37,0,msgsent37,msgbits,itone,iwspr)
88   write(*,*)
89   write(*,'(a9,a37,a3,L2,a7,i2)') 'Message: ',msgsent37,'W:',wspr_hint,' iwspr:',iwspr
90   write(*,1000) f00,xdt,txt,snrdb
911000 format('f0:',f9.3,'   DT:',f6.2,'   TxT:',f6.1,'   SNR:',f6.1)
92   write(*,*)
93   if(i3.eq.1) then
94      write(*,*) '         mycall                         hiscall                    hisgrid'
95      write(*,'(28i1,1x,i1,1x,28i1,1x,i1,1x,i1,1x,15i1,1x,3i1)') msgbits(1:77)
96   else
97      write(*,'(a14)') 'Message bits: '
98      write(*,'(77i1,1x,24i1)') msgbits
99   endif
100   write(*,*)
101   write(*,'(a17)') 'Channel symbols: '
102   write(*,'(10i1)') itone
103   write(*,*)
104
105!   call sgran()
106
107   fsample=12000.0
108   hmod=1
109   icmplx=1
110   f0=f00+1.5*hmod*baud
111   call gen_fst4wave(itone,NN,nsps,nwave,fsample,hmod,f0,icmplx,c0,wave)
112   k=nint((xdt+1.0)/dt)
113   if(nsec.eq.15) k=nint((xdt+0.5)/dt)
114   c0=cshift(c0,-k)
115   if(k.gt.0) c0(0:k-1)=0.0
116   if(k.lt.0) c0(nmax+k:nmax-1)=0.0
117
118   do ifile=1,nfiles
119      c=c0
120      if(fspread.gt.0.0 .or. delay.ne.0.0) call watterson(c,nwave,NZ,fs,delay,fspread)
121      if(fspread.lt.0.0) call lorentzian_fading(c,nwave,fs,-fspread)
122      c=sig*c
123      wave=aimag(c)
124      if(snrdb.lt.90) then
125         do i=1,nmax                   !Add gaussian noise at specified SNR
126            xnoise=gran()
127            wave(i)=wave(i) + xnoise
128         enddo
129      endif
130      gain=100.0
131      if(snrdb.lt.90.0) then
132         wave=gain*wave
133      else
134         datpk=maxval(abs(wave))
135         fac=32766.9/datpk
136         wave=fac*wave
137      endif
138      if(any(abs(wave).gt.32767.0)) print*,"Warning - data will be clipped."
139      iwave=nint(wave(:size(iwave)))
140      h=default_header(12000,nmax)
141      if(nmax/12000.le.30) then
142         write(fname,1102) ifile
1431102     format('000000_',i6.6,'.wav')
144      else
145         write(fname,1104) ifile
1461104     format('000000_',i4.4,'.wav')
147      endif
148      open(10,file=trim(fname),status='unknown',access='stream')
149      write(10) h,iwave                !Save to *.wav file
150      close(10)
151      write(*,1110) ifile,xdt,f00,snrdb,fname
1521110  format(i4,f7.2,f8.2,f7.1,2x,a17)
153   enddo
154
155999 end program fst4sim
156