diff --git a/lib/jt9.f90 b/lib/jt9.f90 index 1f23f3c5d..2139c66aa 100644 --- a/lib/jt9.f90 +++ b/lib/jt9.f90 @@ -249,7 +249,7 @@ program jt9 call timer('symspec ',1) endif nhsym0=nhsym - if(nhsym.ge.181 .and. mode.ne.240 .and. mode.ne.241) exit + if(nhsym.ge.181 .and. mode.ne.240 .and. mode.ne.241 .and. mode.ne.66) exit endif enddo close(unit=wav%lun) diff --git a/lib/qra/qra66/qra66sim.f90 b/lib/qra/qra66/qra66sim.f90 index 0ed23af30..41aa55c50 100644 --- a/lib/qra/qra66/qra66sim.f90 +++ b/lib/qra/qra66/qra66sim.f90 @@ -74,7 +74,7 @@ program qra66sim h=default_header(12000,npts) write(*,1000) -1000 format('File Freq A-E S/N DT Dop Message'/60('-')) +1000 format('File TR Freq Mode S/N DT Dop Message'/60('-')) do ifile=1,nfiles !Loop over requested number of files write(fname,1002) ifile !Output filename @@ -91,11 +91,12 @@ program qra66sim 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,f0,csubmode,xsnr,xdt,fspread,msgsent -1020 format(i4,f10.3,2x,a1,2x,f5.1,f6.2,f6.1,1x,a22) + write(*,1020) ifile,ntrperiod,f0,csubmode,xsnr,xdt,fspread,msgsent +1020 format(i4,i6,f7.1,2x,a1,2x,f5.1,f6.2,f6.1,1x,a22) phi=0.d0 dphi=0.d0 - k=(xdt+0.5)*12000 !Start audio at t = xdt + 0.5 s + 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 @@ -137,16 +138,6 @@ program qra66sim 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. diff --git a/lib/qra66_decode.f90 b/lib/qra66_decode.f90 index c1d6dd804..3996ffdf4 100644 --- a/lib/qra66_decode.f90 +++ b/lib/qra66_decode.f90 @@ -34,7 +34,7 @@ contains use timer_module, only: timer use packjt use, intrinsic :: iso_c_binding - parameter (NMAX=60*12000) !### Needs to be 300*12000 ### + parameter (NMAX=300*12000) !### Needs to be 300*12000 ### class(qra66_decoder), intent(inout) :: this procedure(qra66_decode_callback) :: callback character(len=12) :: mycall, hiscall @@ -43,7 +43,7 @@ contains integer*2 iwave(NMAX) !Raw data integer dat4(12) logical lapdx,ltext - complex c0(0:NMAX-1) !Analytic signal, 6000 S/s + complex, allocatable :: c0(:) !Analytic signal, 6000 S/s real s3(-64:127,63) real s3a(-64:127,63) real a(5) @@ -52,6 +52,8 @@ contains nfft1=ntrperiod*12000 nfft2=ntrperiod*6000 + allocate (c0(0:nfft1-1)) + if(ntrperiod.eq.15) then nsps=1800 else if(ntrperiod.eq.30) then @@ -67,11 +69,12 @@ contains endif baud=12000.0/nsps df1=12000.0/nfft1 - print*,'aaa',ntrperiod,hisgrid,nsps,baud - do i=1,NMAX - write(61,3061) i/12000.0,iwave(i)/32767.0 -3061 format(2f12.6) - enddo + +! do i=1,NMAX +! write(61,3061) i/12000.0,iwave(i)/32767.0 +!3061 format(2f12.6) +! enddo + this%callback => callback if(nutc.eq.-999) print*,lapdx,nfa,nfb,nfqso !Silence warning @@ -107,19 +110,20 @@ contains ! Downsample to give complex data at 6000 S/s fac=2.0/nfft1 - c0=fac*iwave + c0=fac*iwave(1:nfft1) call four2a(c0,nfft1,1,-1,1) !Forward c2c FFT c0(nfft2/2+1:nfft2)=0. !Zero the top half c0(0)=0.5*c0(0) call four2a(c0,nfft2,1,1,1) !Inverse c2c FFT call timer('sync66 ',0) - call sync66(iwave,60*12000,nsps,nfqso,ntol,xdt,f0,snr1) !### 300*12000 ### + call sync66(iwave,ntrperiod*12000,nsps,nfqso,ntol,xdt,f0,snr1) call timer('sync66 ',1) - jpk=(xdt+0.5)*6000 - 384 !### Empirical ### + jpk=(xdt+0.5)*6000 - 384 !### Empirical ### + if(ntrperiod.ge.60) jpk=(xdt+1.0)*6000 - 384 !### TBD ### if(jpk.lt.0) jpk=0 a=0. a(1)=-(f0 + 2.0*baud) !Data tones start 2 bins higher - call twkfreq(c0,c0,15*6000,6000.0,a) + call twkfreq(c0,c0,ntrperiod*6000,6000.0,a) xdt=jpk/6000.0 - 0.5 call spec66(c0(jpk:),nsps/2,s3) diff --git a/lib/sync66.f90 b/lib/sync66.f90 index 6fbd89b27..0c966ebf3 100644 --- a/lib/sync66.f90 +++ b/lib/sync66.f90 @@ -5,7 +5,7 @@ subroutine sync66(iwave,nmax,nsps,nfqso,ntol,xdt,f0,snr1) integer b11(11) !Barker 11 code integer ijpk(2) !Indices i and j at peak of sync_sig real, allocatable :: s1(:,:) !Symbol spectra - real, allocatable :: x(:) !Work array; 2FSK sync modulation + real, allocatable :: x(:) !Work array; demoduated 2FSK sync signal real sync(4*85) !sync vector real sync_sig(-64:64,-15:15) complex, allocatable :: c0(:) !Complex spectrum of symbol @@ -13,11 +13,13 @@ subroutine sync66(iwave,nmax,nsps,nfqso,ntol,xdt,f0,snr1) data sync(1)/99.0/ save sync - nfft=2*NSPS + nfft=2*nsps df=12000.0/nfft istep=nsps/NSTEP - iz=5000.0/df - jz=352 + 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 allocate(s1(iz,jz)) allocate(x(jz)) allocate(c0(0:nsps)) @@ -27,12 +29,13 @@ subroutine sync66(iwave,nmax,nsps,nfqso,ntol,xdt,f0,snr1) do k=1,22 kk=k if(kk.gt.11) kk=k-11 - sync(16*k-15)=2.0*b11(kk) - 1.0 + i=4*NSTEP*(k-1) + 1 + sync(i)=2.0*b11(kk) - 1.0 enddo endif fac=1/32767.0 - do j=1,JZ !Compute symbol spectra + do j=1,jz !Compute symbol spectra at quarter-symbol steps ia=(j-1)*istep ib=ia+nsps-1 k=-1 @@ -44,34 +47,43 @@ subroutine sync66(iwave,nmax,nsps,nfqso,ntol,xdt,f0,snr1) enddo c0(k+1:nfft/2)=0. call four2a(c0,nfft,1,-1,0) !r2c FFT - do i=1,IZ + do i=1,iz s1(i,j)=real(c0(i))**2 + aimag(c0(i))**2 enddo enddo i0=nint(nfqso/df) - call pctile(s1(i0-64:i0+192,1:JZ),129*JZ,40,base) + call pctile(s1(i0-64:i0+192,1:jz),129*jz,40,base) s1=s1/base s1max=20.0 ! Apply AGC - do j=1,JZ + do j=1,jz x(j)=maxval(s1(i0-64:i0+192,j)) if(x(j).gt.s1max) s1(i0-64:i0+192,j)=s1(i0-64:i0+192,j)*s1max/x(j) enddo - dt4=nsps/(NSTEP*12000.0) - j0=0.5/dt4 - sync_sig=0. ia=min(64,nint(ntol/df)) + dt4=nsps/(NSTEP*12000.0) !duration of 1/4 symbol + lag2=nint(0.5/dt4) + if(nsps.ge.7680) lag2=nint(1.0/dt4) + lag1=-lag2 + + jadd=11 + if(nsps.ge.3600) jadd=7 + if(nsps.ge.7680) jadd=6 + if(nsps.ge.16000) jadd=3 + if(nsps.ge.41472) jadd=1 + do i=-ia,ia - x=s1(i0+2+i,:)-s1(i0+i,:) + x=s1(i0+2+i,:)-s1(i0+i,:) !Do the 2FSK demodulation +! do lag=lag1,lag2 do lag=-15,15 -! Make this simpler: just add the 22 nonzero values? - do n=1,4*85 - j=n+lag+11 - if(j.ge.1 .and. j.le.JZ) sync_sig(i,lag)=sync_sig(i,lag) + sync(n)*x(j) + do k=1,22 + n=4*NSTEP*(k-1) + 1 + j=n+lag+jadd + if(j.ge.1 .and. j.le.jz) sync_sig(i,lag)=sync_sig(i,lag) + sync(n)*x(j) enddo enddo enddo @@ -83,6 +95,10 @@ subroutine sync66(iwave,nmax,nsps,nfqso,ntol,xdt,f0,snr1) ! Use peakup() here? f0=nfqso + ii*df jdt=jj + +! j0=0.5/dt4 +! if(nsps.ge.7680) j0=1.0/dt4 + tsym=nsps/12000.0 xdt=jdt*tsym/4.0 @@ -100,4 +116,3 @@ subroutine sync66(iwave,nmax,nsps,nfqso,ntol,xdt,f0,snr1) return end subroutine sync66 -