diff --git a/lib/superfox/sfox_sync.f90 b/lib/superfox/sfox_sync.f90 index c20b341ff..17baf1d25 100644 --- a/lib/superfox/sfox_sync.f90 +++ b/lib/superfox/sfox_sync.f90 @@ -5,13 +5,12 @@ subroutine sfox_sync(iwave,fsample,isync,f,t) integer*2 iwave(0:NMAX-1) integer isync(44) integer ipeak(2) + integer ipeak2(1) complex, allocatable :: c(:) !Work array - real x(171) real, allocatable :: s(:,:) !Symbol spectra, stepped by NSTEP real, allocatable :: savg(:) !Average spectrum real, allocatable :: ccf(:,:) -! character*1 line(-15:15),mark(0:6),c1 -! data mark/' ','.','-','+','X','$','#'/ + real, allocatable :: s2(:) !Fine spectrum of sync tone nfft=nsps nh=nfft/2 @@ -27,11 +26,6 @@ subroutine sfox_sync(iwave,fsample,isync,f,t) lag1=-lagmax lag2=lagmax - x=0. - do i=1,NS - x(isync(i))=1.0 - enddo - allocate(s(0:nh/2,jz)) allocate(savg(0:nh/2)) allocate(c(0:nfft-1)) @@ -71,12 +65,12 @@ subroutine sfox_sync(iwave,fsample,isync,f,t) ccfmax=0. do lag=lag1,lag2 ccft=0. - do kk=1,NS - k=isync(kk) + do m=1,NS + k=isync(m) n=NSTEP*(k-1) + 1 j=n+lag+j0 if(j.ge.1 .and. j.le.jz) ccft=ccft + s(i,j) - enddo ! kk + enddo ! m ccft=ccft - NS*savg(i) ccf(i,lag)=ccft if(ccft.gt.ccfmax) then @@ -96,34 +90,51 @@ subroutine sfox_sync(iwave,fsample,isync,f,t) ipk=ipeak(1)-1+ia jpk=ipeak(2)-1+lag1 - dxi=0. dxj=0. - if(ipk.gt.ia .and. ipk.lt.ib) then - call peakup(ccf(ipk-1,jpk),ccf(ipk,jpk),ccf(ipk+1,jpk),dxi) - endif if(jpk.gt.lag1 .and. jpk.lt.lag2) then call peakup(ccf(ipk,jpk-1),ccf(ipk,jpk),ccf(ipk,jpk+1),dxj) endif f=ibest*df + bw/2 + dxi*df - t=lagbest*dtstep + dxj*dtstep -! write(*,4100) ibest,lagbest,f,dxi*df,t,dxj*dtstep -!4100 format(2i6,2f10.1,2f10.3) + t=(lagbest+dxj)*dtstep + t=t-0.01 !### Why is this needed? ### - nsum=0 - sq=0. - do lag=lag1,lag2 - if(abs(lag-lagbest).gt.3) then - sq=sq + ccf(ibest,lag)**2 - nsum=nsum+1 - endif - write(51,3051) lag*dtstep,ccf(ibest,lag) -3051 format(2f12.4) + nfft2=4*NSPS + deallocate(c) + allocate(c(0:nfft2-1)) + allocate(s2(0:nfft2-1)) + + i0=(t+0.5)*fsample + s2=0. + df2=fsample/nfft2 + do m=1,NS + i1=i0+(isync(m)-1)*NSPS + i2=i1+NSPS-1 + k=-1 + do i=i1,i2,2 !Load iwave data into complex array c0, for r2c FFT + if(i.gt.0) then + xx=iwave(i) + yy=iwave(i+1) + else + xx=0. + yy=0. + endif + k=k+1 + c(k)=fac*cmplx(xx,yy) + enddo + c(k+1:)=0. + call four2a(c,nfft2,1,-1,0) !r2c FFT + do i=1,nfft2/4 + s2(i)=s2(i) + real(c(i))**2 + aimag(c(i))**2 + enddo enddo - - rms=sqrt(sq/nsum) - snrsync=ccf(ibest,lagbest)/rms -! print*,'snr:',snrsync + ipeak2=maxloc(s2) + ipk=ipeak2(1)-1 + dxi=0. + if(ipk.gt.1 .and. ipk.lt.nfft/4) then + call peakup(s2(ipk-1),s2(ipk),s2(ipk+1),dxi) + endif + f=(ipk+dxi)*df2 + bw/2.0 return end subroutine sfox_sync diff --git a/lib/superfox/sfoxtest.f90 b/lib/superfox/sfoxtest.f90 index c7b0508fa..e23c24b8f 100644 --- a/lib/superfox/sfoxtest.f90 +++ b/lib/superfox/sfoxtest.f90 @@ -31,7 +31,6 @@ program sfoxtest integer, allocatable :: rxdat2(:) integer, allocatable :: rxprob2(:) integer, allocatable :: correct(:) - logical hard_sync character fname*17,arg*12,itu*2 data isync/ 1, 2, 5, 11, 19, 24, 26, 28, 29, 35, & @@ -39,13 +38,12 @@ program sfoxtest 80, 82, 84, 85, 92, 98, 103, 107, 109, 111, & 116, 122, 130, 131, 134, 136, 137, 140, 146, 154, & 159, 161, 163, 165/ - data nsb/1,2,4,7,11,16,22,29,37,39/ nargs=iargc() if(nargs.ne.11) then - print*,'Usage: sfoxtest f0 DT ITU M N K NS v hs nfiles snr' - print*,'Example: sfoxtest 1500 0.15 MM 7 127 48 33 0 F 10 -10' + print*,'Usage: sfoxtest f0 DT ITU M N K NS v st nfiles snr' + print*,'Example: sfoxtest 1500 0.15 MM 7 127 48 33 0 3 10 -10' print*,' f0=0 means f0, DT will assume suitable random values' print*,' LQ: Low Latitude Quiet' print*,' MM: Mid Latitude Moderate' @@ -53,7 +51,7 @@ program sfoxtest print*,' ... and similarly for LM LD MQ MD HQ HM' print*,' NS: number of sync symbols' print*,' v=1 for .wav files, 2 for verbose output, 3 for both' - print*,' hs = T for hard-wired sync' + print*,' st: Sync type, 0 for hard-wired, otherwise 1-3' print*,' snr=0 means loop over SNRs' go to 999 endif @@ -73,7 +71,7 @@ program sfoxtest call getarg(8,arg) read(arg,*) nv call getarg(9,arg) - hard_sync=arg(1:1).eq.'T' + read(arg,*) nstype call getarg(10,arg) read(arg,*) nfiles call getarg(11,arg) @@ -86,8 +84,6 @@ program sfoxtest call sfox_init(mm0,nn0,kk0,itu,fspread,delay,fsample,ns0) tsync=NSYNC/fsample txt=(NN+NS)*NSPS/fsample - nstype=nv/10 - nv=mod(nv,10) write(*,1000) MM,NN,KK,NSPS,baud,bw,itu,tsync,txt,nstype 1000 format('M:',i2,' Base code: (',i3,',',i3,') NSPS:',i5, & @@ -202,13 +198,12 @@ program sfoxtest if(snr.lt.90.0) iwave(1:NMAX)=nint(rms*dat(1:NMAX)) crcvd=sig*crcvd+cnoise - if(hard_sync) then - f=f1 ! + 5.0*(ran1(idum)-0.5) - t=xdt ! + 0.01*(ran1(idum)-0.5) + if(nstype.eq.0) then + f=f1 !Hard-wired sync + t=xdt else -! Find signal freq and DT call timer('sync ',0) - call sfox_sync(iwave,fsample,isync,f,t) + call sfox_sync(iwave,fsample,isync,f,t) ! Find signal freq and DT call timer('sync ',1) endif ferr=f-f1 @@ -223,7 +218,7 @@ program sfoxtest endif a=0. - a(1)=1500.0-f - baud + a(1)=1500.0-f - baud !Shift frequencies down by one bin call timer('twkfreq ',0) call twkfreq(crcvd,crcvd,NMAX,fsample,a) call timer('twkfreq ',1)