subroutine short65(data,jz,NFreeze,MouseDF,DFTolerance, + mode65,nspecialbest,nstest,dfsh,iderrbest,idriftbest, + snrdb,ss1,ss2,nwsh) C Checks to see if this might be a shorthand message. C This is done before zapping, downsampling, or normal decoding. parameter (NP2=60*11025) !Size of data array parameter (NFFT=16384) !FFT length parameter (NH=NFFT/2) !Step size parameter (NQ=NFFT/4) !Saved spectral points parameter (MAXSTEPS=60*11025/NH) !Max # of steps real data(jz) integer DFTolerance real s2(NH,MAXSTEPS) !2d spectrum real ss(NQ,4) !Save spectra in four phase bins real psavg(NQ) real sigmax(4) !Peak of spectrum at each phase real ss1(-224:224) !Lower magenta curve real ss2(-224:224) !Upper magenta curve real ssavg(-10:10) integer ipk(4) !Peak bin at each phase save nspecialbest=0 !Default return value nstest=0 df=11025.0/NFFT C Do 16 k FFTs, stepped by 8k. (*** Maybe should step by 4k? ***) call zero(psavg,NQ) nsteps=(jz-NH)/(4*NH) nsteps=4*nsteps !Number of steps do j=1,nsteps k=(j-1)*NH + 1 call ps(data(k),NFFT,s2(1,j)) !Get power spectra call add(psavg,s2(1,j),psavg,NQ) enddo call flat1(psavg,s2,NQ,nsteps,NH,MAXSTEPS) nfac=40*mode65 dtstep=0.5/df fac=dtstep/(60.0*df) C Define range of frequencies to be searched fa= 670.46 fb=1870.46 ia=fa/df ib=fb/df + 4.1*nfac !Upper tone is above sync tone by 4*nfac*df Hz if(ib.gt.NQ) ib=NQ if(NFreeze.eq.1) then fa=max( 670.46,1270.46+MouseDF-DFTolerance) fb=min(1870.46,1270.46+MouseDF+DFTolerance) endif ia2=fa/df ib2=fb/df + 4.1*nfac !Upper tone is above sync tone by 4*nfac*df Hz if(ib2.gt.NQ) ib2=NQ C Find strongest line in each of the 4 phases, repeating for each drift rate. sbest=0. snrbest=0. idz=6.0/df !Is this the right drift range? do idrift=-idz,idz drift=idrift*df*60.0/49.04 call zero(ss,4*NQ) !Clear the accumulating array do j=1,nsteps n=mod(j-1,4)+1 k=nint((j-nsteps/2)*drift*fac) + ia call add(ss(ia,n),s2(k,j),ss(ia,n),ib-ia+1) enddo do n=1,4 sigmax(n)=0. do i=ia2,ib2 sig=ss(i,n) if(sig.ge.sigmax(n)) then ipk(n)=i sigmax(n)=sig if(sig.ge.sbest) then sbest=sig nbest=n fdotsh=drift endif endif enddo enddo n2best=nbest+2 if(n2best.gt.4) n2best=nbest-2 xdf=min(ipk(nbest),ipk(n2best))*df - 1270.46 if(NFreeze.eq.1 .and. abs(xdf-mousedf).gt.DFTolerance) goto 10 idiff=abs(ipk(nbest)-ipk(n2best)) xk=float(idiff)/nfac k=nint(xk) iderr=nint((xk-k)*nfac) nspecial=0 maxerr=nint(0.008*abs(idiff) + 0.51) if(abs(iderr).le.maxerr .and. k.ge.2 .and. k.le.4) nspecial=k if(nspecial.gt.0) then call getsnr(ss(ia2,nbest),ib2-ia2+1,snr1) call getsnr(ss(ia2,n2best),ib2-ia2+1,snr2) snr=0.5*(snr1+snr2) if(snr.gt.snrbest) then snrbest=snr nspecialbest=nspecial nstest=snr/2.0 - 2.0 !Threshold set here if(nstest.lt.0) nstest=0 if(nstest.gt.10) nstest=10 dfsh=nint(xdf) fdotbest=fdotsh iderrbest=iderr idiffbest=idiff idriftbest=idrift snrdb=db(snr) - db(2500.0/df) - db(sqrt(nsteps/4.0))+1.8 n1=nbest n2=n2best ipk1=ipk(n1) ipk2=ipk(n2) endif endif if(nstest.eq.0) nspecial=0 10 enddo if(nstest.eq.0) nspecialbest=0 if(nstest.gt.0) then df4=4.0*df if(ipk1.gt.ipk2) then ntmp=n1 n1=n2 n2=ntmp ntmp=ipk1 ipk1=ipk2 ipk2=ntmp endif call zero(ss1,449) call zero(ss2,449) do i=ia2,ib2,4 f=df*i k=nint((f-1270.46)/df4) if(k.ge.-224 .and. k.le.224) then ss1(k)=0.3 * (ss(i-2,n1) + ss(i-1,n1) + ss(i,n1) + + ss(i+1,n1) + ss(i+2,n1)) ss2(k)=0.3 * (ss(i-2,n2) + ss(i-1,n2) + ss(i,n2) + + ss(i+1,n2) + ss(i+2,n2)) endif enddo ! kpk1=nint(0.25*ipk(n1)-472.0) kpk1=nint(0.25*ipk1-472.0) kpk2=kpk1 + nspecial*mode65*10 ssmax=0. do i=-10,10 ssavg(i)=ss1(kpk1+i) + ss2(kpk2+i) if(ssavg(i).gt.ssmax) then ssmax=ssavg(i) itop=i endif enddo base=0.25*(ssavg(-10)+ssavg(-9)+ssavg(9)+ssavg(10)) shalf=0.5*(ssmax+base) do k=1,8 if(ssavg(itop-k).lt.shalf) go to 110 enddo k=8 110 x=(ssavg(itop-(k-1))-shalf)/(ssavg(itop-(k-1))-ssavg(itop-k)) do k=1,8 if(ssavg(itop+k).lt.shalf) go to 120 enddo k=8 120 x=x+(ssavg(itop+(k-1))-shalf)/(ssavg(itop+(k-1))-ssavg(itop+k)) nwsh=nint(x*df4) endif return end