diff --git a/lib/qra/qra64/qra64sim.f90 b/lib/qra/qra64/qra64sim.f90 index 278c1797d..2360ef7da 100644 --- a/lib/qra/qra64/qra64sim.f90 +++ b/lib/qra/qra64/qra64sim.f90 @@ -14,14 +14,18 @@ program qra64sim complex cdat(NMAX) !Generated complex waveform complex cspread(0:NFFT-1) !Complex amplitude for Rayleigh fading complex z + complex c00(0:720000) !Analytic signal for dat() real*8 f0,dt,twopi,phi,dphi,baud,fsample,freq - character msg*22,fname*11,csubmode*1,arg*12 + character msg*22,fname*11,csubmode*1,arg*12,cd*1 character msgsent*22 + logical lsync + data lsync/.false./ nargs=iargc() - if(nargs.ne. 7) then - print *, 'Usage: qra64sim "msg" A-E Nsigs fDop DT Nfiles SNR' - print *, 'Example qra64sim "K1ABC W9XYZ EN37" A 10 0.2 0.0 1 0' + if(nargs.ne.8) then + print*,'Usage: qra64sim "msg" A-E Nsigs fDop DT Nfiles Sync SNR' + print*,'Example qra64sim "K1ABC W9XYZ EN37" A 10 0.2 0.0 1 T -26' + print*,'Sync = T to include sync test.' go to 999 endif call getarg(1,msg) @@ -36,8 +40,10 @@ program qra64sim call getarg(6,arg) read(arg,*) nfiles call getarg(7,arg) + if(arg(1:1).eq.'T' .or. arg(1:1).eq.'1') lsync=.true. + call getarg(8,arg) read(arg,*) snrdb - + if(mode64.ge.8) nsigs=1 rms=100. fsample=12000.d0 !Sample rate (Hz) @@ -54,6 +60,7 @@ program qra64sim write(*,1000) 1000 format('File Sig Freq A-E S/N DT Dop Message'/60('-')) + nsync=0 do ifile=1,nfiles !Loop over requested number of files write(fname,1002) ifile !Output filename 1002 format('000000_',i4.4) @@ -165,6 +172,30 @@ program qra64sim 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) + + if(lsync) then + cd=' ' + if(ifile.eq.nfiles) cd='d' + nf1=200 + nf2=3000 + nfqso=nint(f0) + ntol=100 + minsync=0 + emedelay=0.0 + call ana64(dat,npts,c00) + call sync64(c00,nf1,nf2,nfqso,ntol,minsync,mode64,emedelay,xdt2,f02, & + jpk0,sync,sync2,width) + terr=1.01/(8.0*baud) + ferr=1.01*mode64*baud + if(abs(xdt2-xdt).lt.terr .and. abs(f02-f0).lt.ferr) nsync=nsync+1 + open(40,file='sync64.out',status='unknown',position='append') + write(40,1030) ifile,64,csubmode,snrdb,fspread,xdt2-xdt,f02-f0, & + width,sync,sync2,nsync,cd +1030 format(i4,i3,1x,a1,2f7.1,f7.2,4f8.1,i5,1x,a1) + close(40) + endif enddo + if(lsync) write(*,1040) snrdb,nfiles,nsync +1040 format('SNR:',f6.1,' nfiles:',i5,' nsynced:',i5) 999 end program qra64sim diff --git a/lib/qra/qra65/qra65sim.f90 b/lib/qra/qra65/qra65sim.f90 index e91324d27..62ec63f57 100644 --- a/lib/qra/qra65/qra65sim.f90 +++ b/lib/qra/qra65/qra65sim.f90 @@ -14,13 +14,16 @@ program qra65sim complex cspread(0:NMAX-1) !Complex amplitude for Rayleigh fading complex z real*8 f0,dt,twopi,phi,dphi,baud,fsample,freq - character msg*22,fname*17,csubmode*1,arg*12 + character msg*22,fname*17,csubmode*1,arg*12,cd*1 character msgsent*22 - + logical lsync + data lsync/.false./ + nargs=iargc() - if(nargs.ne.8) then - print *, 'Usage: qra65sim "msg" A-E freq fDop DT TRp Nfiles SNR' - print *, 'Example: qra65sim "K1ABC W9XYZ EN37" A 1500 0.2 0.0 15 1 -10' + if(nargs.ne.9) then + print *, 'Usage: qra65sim "msg" A-E freq fDop DT TRp Nfiles Sync SNR' + print *, 'Example: qra65sim "K1ABC W9XYZ EN37" A 1500 0.0 0.0 60 1 T -26' + print*,'Sync = T to include sync test.' go to 999 endif call getarg(1,msg) @@ -37,8 +40,15 @@ program qra65sim call getarg(7,arg) read(arg,*) nfiles call getarg(8,arg) + if(arg(1:1).eq.'T' .or. arg(1:1).eq.'1') lsync=.true. + call getarg(9,arg) read(arg,*) snrdb + if(nfiles.lt.0) then + nfiles=-nfiles + lsync=.true. + endif + if(ntrperiod.eq.15) then nsps=1800 else if(ntrperiod.eq.30) then @@ -75,6 +85,7 @@ program qra65sim write(*,1000) 1000 format('File TR Freq Mode S/N DT Dop Message'/60('-')) + nsync=0 do ifile=1,nfiles !Loop over requested number of files if(ntrperiod.lt.60) then write(fname,1002) ifile !Output filename @@ -169,10 +180,23 @@ program qra65sim write(10) h,iwave(1:npts) !Save the .wav file close(10) -! do i=1,NMAX -! write(15,3020) i/12000.0,iwave(i) -!3020 format(f10.6,i8) -! enddo + if(lsync) then + cd=' ' + if(ifile.eq.nfiles) cd='d' + nfqso=nint(f0) + ntol=100 + call sync_qra65(iwave,npts,mode65,nsps,nfqso,ntol,xdt2,f02,snr2) + terr=1.01/(8.0*baud) + ferr=1.01*mode65*baud + if(abs(xdt2-xdt).lt.terr .and. abs(f02-f0).lt.ferr) nsync=nsync+1 + open(40,file='sync65.out',status='unknown',position='append') + write(40,1030) ifile,65,csubmode,snrdb,fspread,xdt2-xdt,f02-f0, & + snr2,nsync,cd +1030 format(i4,i3,1x,a1,2f7.1,f7.2,2f8.1,i5,1x,a1) + close(40) + endif enddo + if(lsync) write(*,1040) snrdb,nfiles,nsync +1040 format('SNR:',f6.1,' nfiles:',i5,' nsynced:',i5) 999 end program qra65sim diff --git a/lib/qra65_decode.f90 b/lib/qra65_decode.f90 index 7778844bf..c5ab8521c 100644 --- a/lib/qra65_decode.f90 +++ b/lib/qra65_decode.f90 @@ -93,8 +93,9 @@ contains ! AP control could be done differently, but this works well: maxaptype=0 - if(ndepth.eq.2) maxaptype=3 - if(ndepth.eq.3) maxaptype=5 +! if(ndepth.eq.2) maxaptype=3 +! if(ndepth.eq.3) maxaptype=5 + if(ndepth.ge.2) maxaptype=5 !### if(nc1.ne.nc1z .or. nc2.ne.nc2z .or. ng2.ne.ng2z .or. & maxaptype.ne.maxaptypez) then diff --git a/lib/qra_loops.f90 b/lib/qra_loops.f90 index d59bbe05f..59b9458b1 100644 --- a/lib/qra_loops.f90 +++ b/lib/qra_loops.f90 @@ -69,7 +69,6 @@ subroutine qra_loops(c00,npts2,mode,mode64,nsubmode,nFadingModel,minsync, & endif enddo ! ibw (b90 loop) if(iand(ndepth,3).lt.3 .and. irc.ge.0) go to 100 - if(irc.eq.0) go to 100 enddo ! idt (DT loop) enddo ! idf (f0 loop) @@ -85,5 +84,9 @@ subroutine qra_loops(c00,npts2,mode,mode64,nsubmode,nFadingModel,minsync, & ibw=ibwkeep endif -200 return +200 continue + write(53,3053) idt,idf,ibw,irc,b90,xdt,f0,snr2 +3053 format(4i5,f7.1,f7.2,2f7.1) + + return end subroutine qra_loops diff --git a/lib/sync_qra65.f90 b/lib/sync_qra65.f90 index c93aed626..b18aa8870 100644 --- a/lib/sync_qra65.f90 +++ b/lib/sync_qra65.f90 @@ -1,6 +1,8 @@ subroutine sync_qra65(iwave,nmax,mode65,nsps,nfqso,ntol,xdt,f0,snr1) -! Look for the sync vector in a QRA65 signal. +! Detect and align with the QRA65 sync vector, returning time and frequency +! offsets and SNR estimate. + ! Input: iwave(0:nmax-1) Raw data ! mode65 Tone spacing 1 2 4 8 16 (A-E) ! nsps Samples per symbol at 12000 Sa/s @@ -10,7 +12,7 @@ subroutine sync_qra65(iwave,nmax,mode65,nsps,nfqso,ntol,xdt,f0,snr1) ! f0 Frequency of sync tone ! snr1 Relative SNR of sync signal - parameter (NSTEP=4) !Quarter-symbol steps + parameter (NSTEP=8) !Step size nsps/NSTEP integer*2 iwave(0:nmax-1) !Raw data integer isync(22) !Indices of sync symbols integer ijpk(2) !Indices i and j at peak of ccf @@ -28,7 +30,7 @@ subroutine sync_qra65(iwave,nmax,mode65,nsps,nfqso,ntol,xdt,f0,snr1) iz=5000.0/df !Uppermost frequency bin, at 5000 Hz txt=85.0*nsps/12000.0 jz=(txt+1.0)*12000.0/istep !Number of quarter-symbol steps - if(nsps.ge.7680) jz=(txt+2.0)*12000.0/istep !For TR 60 s and higher + if(nsps.ge.6912) jz=(txt+2.0)*12000.0/istep !For TR 60 s and higher allocate(s1(iz,jz)) allocate(c0(0:nfft-1)) @@ -41,11 +43,11 @@ subroutine sync_qra65(iwave,nmax,mode65,nsps,nfqso,ntol,xdt,f0,snr1) endif fac=1/32767.0 - do j=1,jz !Compute symbol spectra at quarter-symbol steps + do j=1,jz !Compute symbol spectra at step size ia=(j-1)*istep ib=ia+nsps-1 k=-1 - do i=ia,ib,2 + do i=ia,ib,2 !Load iwave data into complex array c0, for r2c FFT xx=iwave(i) yy=iwave(i+1) k=k+1 @@ -57,11 +59,12 @@ subroutine sync_qra65(iwave,nmax,mode65,nsps,nfqso,ntol,xdt,f0,snr1) s1(i,j)=real(c0(i))**2 + aimag(c0(i))**2 enddo ! For large Doppler spreads, should we smooth the spectra here? + call smo121(s1(1:iz,j),iz) enddo i0=nint(nfqso/df) !Target QSO frequency call pctile(s1(i0-64:i0+192,1:jz),129*jz,40,base) - s1=s1/base !Maybe should subtract 1.0 here? + s1=s1/base - 1.0 ! Apply fast AGC s1max=20.0 !Empirical choice @@ -70,24 +73,22 @@ subroutine sync_qra65(iwave,nmax,mode65,nsps,nfqso,ntol,xdt,f0,snr1) if(smax.gt.s1max) s1(i0-64:i0+192,j)=s1(i0-64:i0+192,j)*s1max/smax enddo - dt4=nsps/(NSTEP*12000.0) !1/4 of symbol duration - j0=0.5/dt4 - if(nsps.ge.7680) j0=1.0/dt4 !Nominal index for start of signal + dtstep=nsps/(NSTEP*12000.0) !Step size in seconds + j0=0.5/dtstep + if(nsps.ge.6192) j0=1.0/dtstep !Nominal index for start of signal ccf=0. ia=min(64,nint(ntol/df)) - lag1=-1.0/dt4 - lag2=4.0/dt4 + 0.9999 + lag1=-1.0/dtstep + lag2=4.0/dtstep + 0.9999 - do i=-ia,ia - do lag=lag1,lag2 - do k=1,85 - n=NSTEP*(k-1) + 1 - j=n+lag+j0 - if(j.ge.1 .and. j.le.jz) then - ccf(i,lag)=ccf(i,lag) + sync(k)*s1(i0+i,j) - endif - enddo + do lag=lag1,lag2 + do k=1,85 + n=NSTEP*(k-1) + 1 + j=n+lag+j0 + if(j.ge.1 .and. j.le.jz) then + ccf(-ia:ia,lag)=ccf(-ia:ia,lag) + sync(k)*s1(i0-ia:i0+ia,j) + endif enddo enddo @@ -95,7 +96,7 @@ subroutine sync_qra65(iwave,nmax,mode65,nsps,nfqso,ntol,xdt,f0,snr1) ipk=ijpk(1)-65 jpk=ijpk(2)-27 f0=nfqso + ipk*df - xdt=jpk*dt4 + xdt=jpk*dtstep snr1=maxval(ccf)/22.0 return