Add detailed comments to get_candidates_fst240.f90.

This commit is contained in:
Joe Taylor 2020-06-29 09:37:29 -04:00
parent 494481fa8a
commit 4fbed923ab

View File

@ -244,6 +244,7 @@ contains
endif endif
enddo enddo
ncand=ic ncand=ic
xsnr=0.
do icand=1,ncand do icand=1,ncand
sync=candidates(icand,2) sync=candidates(icand,2)
@ -351,8 +352,8 @@ contains
iaptype=0 iaptype=0
qual=0. qual=0.
fsig=fc_synced - 1.5*hmod*baud fsig=fc_synced - 1.5*hmod*baud
write(21,'(i6,7i6,f7.1,f9.2,3f7.1,1x,a37)') & !write(21,'(i6,7i6,f7.1,f9.2,3f7.1,1x,a37)') &
nutc,icand,itry,iaptype,ijitter,ntype,nsync_qual,nharderrors,dmin,sync,xsnr,xdt,fsig,msg ! nutc,icand,itry,iaptype,ijitter,ntype,nsync_qual,nharderrors,dmin,sync,xsnr,xdt,fsig,msg
call this%callback(nutc,smax1,nsnr,xdt,fsig,msg, & call this%callback(nutc,smax1,nsnr,xdt,fsig,msg, &
iaptype,qual,ntrperiod) iaptype,qual,ntrperiod)
goto 2002 goto 2002
@ -479,77 +480,79 @@ write(21,'(i6,7i6,f7.1,f9.2,3f7.1,1x,a37)') &
subroutine get_candidates_fst240(c_bigfft,nfft1,nsps,hmod,fs,fa,fb, & subroutine get_candidates_fst240(c_bigfft,nfft1,nsps,hmod,fs,fa,fb, &
ncand,candidates,base) ncand,candidates,base)
complex c_bigfft(0:nfft1/2) complex c_bigfft(0:nfft1/2) !Full length FFT of raw data
integer hmod integer hmod !Modulation index (submode)
integer indx(100),im(1) integer im(1) !For maxloc
real candidates(100,4) real candidates(100,4) !Candidate list
real candidates0(100,4) real s(18000) !Low resolution power spectrum
real snr_cand(100) real s2(18000) !CCF of s() with 4 tones
real s(18000) real xdb(-3:3) !Model 4-tone CCF peaks
real s2(18000)
real xdb(-3:3)
data xdb/0.25,0.50,0.75,1.0,0.75,0.50,0.25/ data xdb/0.25,0.50,0.75,1.0,0.75,0.50,0.25/
data nfft1z/-1/ data nfft1z/-1/
save nfft1z save nfft1z
nh1=nfft1/2 nh1=nfft1/2
df1=fs/nfft1 df1=fs/nfft1
baud=fs/nsps baud=fs/nsps !Keying rate
df2=baud/2.0 df2=baud/2.0
nd=df2/df1 nd=df2/df1 !s() sums this many bins of big FFT
ndh=nd/2 ndh=nd/2
ia=nint(max(100.0,fa)/df2) ia=nint(max(100.0,fa)/df2) !Low frequency search limit
ib=nint(min(4800.0,fb)/df2) ib=nint(min(4800.0,fb)/df2) !High frequency limit
signal_bw=4*(12000.0/nsps)*hmod signal_bw=4*(12000.0/nsps)*hmod
analysis_bw=min(4800.0,fb)-max(100.0,fa) analysis_bw=min(4800.0,fb)-max(100.0,fa)
noise_bw=10.0*signal_bw noise_bw=10.0*signal_bw !Is this a good compromise?
if(analysis_bw.gt.noise_bw) then if(analysis_bw.gt.noise_bw) then
ina=ia ina=ia
inb=ib inb=ib
else else
fcenter=(fa+fb)/2.0 fcenter=(fa+fb)/2.0 !If noise_bw > analysis_bw,
fl = max(100.0,fcenter-noise_bw/2.)/df2 fl = max(100.0,fcenter-noise_bw/2.)/df2 !we'll search over noise_bw
fh = min(4800.0,fcenter+noise_bw/2.)/df2 fh = min(4800.0,fcenter+noise_bw/2.)/df2
ina=nint(fl) ina=nint(fl)
inb=nint(fh) inb=nint(fh)
endif endif
s=0.
s=0. !Compute low-resloution power spectrum
do i=ina,inb ! noise analysis window includes signal analysis window do i=ina,inb ! noise analysis window includes signal analysis window
j0=nint(i*df2/df1) j0=nint(i*df2/df1)
do j=j0-ndh,j0+ndh do j=j0-ndh,j0+ndh
s(i)=s(i) + real(c_bigfft(j))**2 + aimag(c_bigfft(j))**2 s(i)=s(i) + real(c_bigfft(j))**2 + aimag(c_bigfft(j))**2
enddo enddo
enddo enddo
ina=max(ina,1+3*hmod)
ina=max(ina,1+3*hmod) !Don't run off the ends
inb=min(inb,18000-3*hmod) inb=min(inb,18000-3*hmod)
s2=0. s2=0.
do i=ina,inb do i=ina,inb !Compute CCF of s() and 4 tones
s2(i)=s(i-hmod*3) + s(i-hmod) +s(i+hmod) +s(i+hmod*3) s2(i)=s(i-hmod*3) + s(i-hmod) +s(i+hmod) +s(i+hmod*3)
enddo enddo
call pctile(s2(ina+hmod*3:inb-hmod*3),inb-ina+1-hmod*6,30,base) call pctile(s2(ina+hmod*3:inb-hmod*3),inb-ina+1-hmod*6,30,base)
s2=s2/base s2=s2/base !Normalize wrt noise level
thresh=1.25 thresh=1.25 !First candidate threshold
ncand=0 ncand=0
candidates=0 candidates=0
if(ia.lt.3) ia=3 if(ia.lt.3) ia=3
if(ib.gt.18000-2) ib=18000-2 if(ib.gt.18000-2) ib=18000-2
! Find candidates, using the CLEAN algorithm to remove a model of each one
! from s2() after it has been found.
pval=99.99 pval=99.99
do while(ncand.lt.100 .and. pval.gt.thresh) do while(ncand.lt.100 .and. pval.gt.thresh)
im=maxloc(s2(ia:ib)) im=maxloc(s2(ia:ib))
iploc=ia+im(1)-1 iploc=ia+im(1)-1 !Index of CCF peak
pval=s2(iploc) pval=s2(iploc) !Peak value
if(s2(iploc).gt.thresh) then if(s2(iploc).gt.thresh) then !Is this a possible candidate?
do i=-3,+3 do i=-3,+3 !Remove 0.9 of a model CCF at
k=iploc+2*hmod*i k=iploc+2*hmod*i !this frequency from s2()
if(k.ge.ia .and. k.le.ib) then if(k.ge.ia .and. k.le.ib) then
s2(k)=max(0.,s2(k)-0.9*pval*xdb(i)) s2(k)=max(0.,s2(k)-0.9*pval*xdb(i))
endif endif
enddo enddo
ncand=ncand+1 ncand=ncand+1
candidates(ncand,1)=df2*iploc candidates(ncand,1)=df2*iploc !Candidate frequency
candidates(ncand,2)=pval candidates(ncand,2)=pval !Rough estimate of SNR
endif endif
enddo enddo