diff --git a/CMakeLists.txt b/CMakeLists.txt index 79c8302e5..52b69b84a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -463,7 +463,7 @@ set (wsjt_FSRCS lib/genmsk_128_90.f90 lib/genmsk40.f90 lib/fsk4hf/genft2.f90 - lib/fsk4hf/genft4.f90 + lib/ft4/genft4.f90 lib/genqra64.f90 lib/ft8/genft8refsig.f90 lib/genwspr.f90 @@ -509,8 +509,10 @@ set (wsjt_FSRCS lib/msk144signalquality.f90 lib/msk144sim.f90 lib/mskrtd.f90 - lib/fsk4hf/ft4sim.f90 - lib/fsk4hf/ft4d.f90 + lib/fsk4hf/ft2sim.f90 + lib/fsk4hf/ft2d.f90 + lib/ft4/ft4sim.f90 + lib/ft4/ft4d.f90 lib/ft2/cdatetime.f90 lib/ft2/ft2_decode.f90 lib/77bit/my_hash.f90 @@ -553,10 +555,13 @@ set (wsjt_FSRCS lib/sync4.f90 lib/sync64.f90 lib/sync65.f90 + lib/ft4/getcandidates4.f90 lib/fsk4hf/getcandidates2.f90 lib/ft2/getcandidates2a.f90 + lib/ft4/syncft4.f90 lib/ft8/sync8.f90 lib/ft8/sync8d.f90 + lib/ft4/sync4d.f90 lib/sync9.f90 lib/sync9f.f90 lib/sync9w.f90 @@ -1265,10 +1270,16 @@ target_link_libraries (ft8sim wsjt_fort wsjt_cxx) add_executable (msk144sim lib/msk144sim.f90 wsjtx.rc) target_link_libraries (msk144sim wsjt_fort wsjt_cxx) -add_executable (ft4sim lib/fsk4hf/ft4sim.f90 wsjtx.rc) +add_executable (ft2sim lib/fsk4hf/ft2sim.f90 wsjtx.rc) +target_link_libraries (ft2sim wsjt_fort wsjt_cxx) + +add_executable (ft2d lib/fsk4hf/ft2d.f90 wsjtx.rc) +target_link_libraries (ft2d wsjt_fort wsjt_cxx) + +add_executable (ft4sim lib/ft4/ft4sim.f90 wsjtx.rc) target_link_libraries (ft4sim wsjt_fort wsjt_cxx) -add_executable (ft4d lib/fsk4hf/ft4d.f90 wsjtx.rc) +add_executable (ft4d lib/ft4/ft4d.f90 wsjtx.rc) target_link_libraries (ft4d wsjt_fort wsjt_cxx) endif(WSJT_BUILD_UTILS) diff --git a/lib/fsk4hf/ft4d.f90 b/lib/fsk4hf/ft4d.f90 index d41853075..d1bf262c1 100644 --- a/lib/fsk4hf/ft4d.f90 +++ b/lib/fsk4hf/ft4d.f90 @@ -1,317 +1,329 @@ program ft4d - use crc - use packjt77 - include 'ft4_params.f90' - character arg*8,message*37,c77*77,infile*80,fname*16,datetime*11 - character*37 decodes(100) - character*120 data_dir - character*90 dmsg - complex cd2(0:NMAX/16-1) !Complex waveform - complex cb(0:NMAX/16-1) - complex cd(0:76*20-1) !Complex waveform - complex c3(0:19),c2(0:19),c1(0:19),c0(0:19) - complex ccor(0:3,76) - complex csum,cterm,cc0,cc1,cc2,cc3,csync1,csync2 - complex csync(12) - real*8 fMHz + use crc + use packjt77 + include 'ft4_params.f90' + character arg*8,message*37,c77*77,infile*80,fname*16,datetime*11 + character*37 decodes(100) + character*120 data_dir + character*90 dmsg + complex cd2(0:NMAX/16-1) !Complex waveform + complex cb(0:NMAX/16-1) + complex cd(0:76*20-1) !Complex waveform + complex csum,cterm + complex ctwk(80),ctwk2(80) + complex csymb(20) + complex cs(0:3,NN) + real s4(0:3,NN) - real a(5) - real rxdata(128),llr(128) !Soft symbols - real llr2(128) - real sbits(152),sbits1(152),sbits3(152) - real ps(0:8191),psbest(0:8191) - real candidates(100,2) - real savg(NH1),sbase(NH1) - integer ihdr(11) - integer*2 iwave(NMAX) !Generated full-length waveform - integer*1 message77(77),apmask(128),cw(128) - integer*1 hbits(152),hbits1(152),hbits3(152) - integer*1 s12(12) - logical unpk77_success - data s12/0,0,0,1,1,1,1,1,1,0,0,0/ + real*8 fMHz + real ps(0:8191),psbest(0:8191) + real bmeta(152),bmetb(152),bmetc(152) + real s(NH1,NHSYM) + real a(5) + real llr(128),llr2(128),llra(128),llrb(128),llrc(128) + real s2(0:255) + real candidate(3,100) + real savg(NH1),sbase(NH1) + integer ihdr(11) + integer icos4(0:3) + integer*2 iwave(NMAX) !Generated full-length waveform + integer*1 message77(77),apmask(128),cw(128) + integer*1 hbits(152),hbits1(152),hbits3(152) + integer*1 s12(12) + integer graymap(0:3) + integer ip(1) + logical unpk77_success + logical one(0:511,0:7) ! 256 4-symbol sequences, 8 bits + data s12/1,1,1,2,2,2,2,2,2,1,1,1/ + data icos4/0,1,3,2/ + data graymap/0,1,3,2/ + save one - fs=12000.0/NDOWN !Sample rate - dt=1/fs !Sample interval after downsample (s) - tt=NSPS*dt !Duration of "itone" symbols (s) - baud=1.0/tt !Keying rate for "itone" symbols (baud) - txt=NZ*dt !Transmission length (s) - twopi=8.0*atan(1.0) - h=1.000 !h=0.8 seems to be optimum for AWGN sensitivity (not for fading) + fs=12000.0/NDOWN !Sample rate + dt=1/fs !Sample interval after downsample (s) + tt=NSPS*dt !Duration of "itone" symbols (s) + baud=1.0/tt !Keying rate for "itone" symbols (baud) + txt=NZ*dt !Transmission length (s) + twopi=8.0*atan(1.0) + h=1.0 !h=0.8 seems to be optimum for AWGN sensitivity (not for fading) - dphi=twopi/2*baud*h*dt*16 ! dt*16 is samp interval after downsample - dphi0=-3*dphi - dphi1=-dphi - dphi2=+dphi - dphi3=+3*dphi - phi0=0.0 - phi1=0.0 - phi2=0.0 - phi3=0.0 - do i=0,19 - c3(i)=cmplx(cos(phi3),sin(phi3)) - c2(i)=cmplx(cos(phi2),sin(phi2)) - c1(i)=cmplx(cos(phi1),sin(phi1)) - c0(i)=cmplx(cos(phi0),sin(phi0)) - phi3=mod(phi3+dphi3,twopi) - phi2=mod(phi2+dphi2,twopi) - phi1=mod(phi1+dphi1,twopi) - phi0=mod(phi0+dphi0,twopi) - enddo - the=twopi*h/2.0 - cc3=cmplx(cos(3*the),+sin(3*the)) - cc2=cmplx(cos(the),+sin(the)) - cc1=cmplx(cos(the),-sin(the)) - cc0=cmplx(cos(3*the),-sin(3*the)) + one=.false. + do i=0,255 + do j=0,7 + if(iand(i,2**j).ne.0) one(i,j)=.true. + enddo + enddo + + nargs=iargc() + if(nargs.lt.1) then + print*,'Usage: ft4d [-a ] [-f fMHz] file1 [file2 ...]' + go to 999 + endif + iarg=1 + data_dir="." + call getarg(iarg,arg) + if(arg(1:2).eq.'-a') then + call getarg(iarg+1,data_dir) + iarg=iarg+2 + endif + call getarg(iarg,arg) + if(arg(1:2).eq.'-f') then + call getarg(iarg+1,arg) + read(arg,*) fMHz + iarg=iarg+2 + endif + ncoh=1 + + do ifile=iarg,nargs + call getarg(ifile,infile) + j2=index(infile,'.wav') + open(10,file=infile,status='old',access='stream') + read(10,end=999) ihdr,iwave + read(infile(j2-4:j2-1),*) nutc + datetime=infile(j2-11:j2-1) + close(10) + candidate=0.0 + ncand=0 + + nfqso=1500 + nfa=500 + nfb=2700 + syncmin=1.0 + maxcand=100 +! call syncft4(iwave,nfa,nfb,syncmin,nfqso,maxcand,s,candidate,ncand,sbase) + + call getcandidates4(iwave,375.0,3000.0,0.2,2200.0,100,savg,candidate,ncand,sbase) + ndecodes=0 + do icand=1,ncand + f0=candidate(1,icand)-1.5*37.5 + xsnr=1.0 + if( f0.le.375.0 .or. f0.ge.(5000.0-375.0) ) cycle + call ft4_downsample(iwave,f0,cd2) ! downsample from 320 Sa/Symbol to 20 Sa/Symbol + sum2=sum(cd2*conjg(cd2))/(20.0*76) + if(sum2.gt.0.0) cd2=cd2/sqrt(sum2) - nargs=iargc() - if(nargs.lt.1) then - print*,'Usage: ft4d [-a ] [-f fMHz] file1 [file2 ...]' - go to 999 - endif - iarg=1 - data_dir="." - call getarg(iarg,arg) - if(arg(1:2).eq.'-a') then - call getarg(iarg+1,data_dir) - iarg=iarg+2 - endif - call getarg(iarg,arg) - if(arg(1:2).eq.'-f') then - call getarg(iarg+1,arg) - read(arg,*) fMHz - iarg=iarg+2 - endif - ncoh=1 - - do ifile=iarg,nargs - call getarg(ifile,infile) - j2=index(infile,'.wav') - open(10,file=infile,status='old',access='stream') - read(10,end=999) ihdr,iwave - read(infile(j2-4:j2-1),*) nutc - datetime=infile(j2-11:j2-1) - close(10) - candidates=0.0 - ncand=0 - call getcandidates2(iwave,375.0,3000.0,0.2,2200.0,100,savg,candidates,ncand,sbase) - ndecodes=0 - do icand=1,ncand - f0=candidates(icand,1) - xsnr=1.0 - if( f0.le.375.0 .or. f0.ge.(5000.0-375.0) ) cycle - call ft4_downsample(iwave,f0,cd2) ! downsample from 320s/Symbol to 20s/Symbol - s2=sum(cd2*conjg(cd2))/(20.0*76) - if(s2.gt.0.0) cd2=cd2/sqrt(s2) ! 750 samples/second here - ibest=-1 - sybest=-99. - dfbest=-1. - do if=-30,+30 - df=if - a=0. - a(1)=-df - call twkfreq1(cd2,NMAX/16,fs,a,cb) - do istart=0,380 - csync1=0. - cterm=1 - do ib=1,12 - i1=(ib-1)*20+istart - if(s12(ib).eq.0) then - csync1=csync1+sum(cb(i1:i1+19)*conjg(c0(0:19)))*cterm - cterm=cterm*conjg(cc0) - else - csync1=csync1+sum(cb(i1:i1+19)*conjg(c3(0:19)))*cterm - cterm=cterm*conjg(cc3) - endif - enddo - if(abs(csync1).gt.sybest) then - ibest=istart - sybest=abs(csync1) - dfbest=df - endif - enddo - enddo - - a=0. -!dfbest=1500.0-f0 - a(1)=-dfbest - - call twkfreq1(cd2,NMAX/16,fs,a,cb) + ibest=-1 + smax=-99. + dfbest=-1. + do idf=-90,+90,5 + df=idf + a=0. + a(1)=df + ctwk=1. + call twkfreq1(ctwk,80,fs,a,ctwk2) + do istart=0,315 + call sync4d(cd2,istart,ctwk2,1,sync) + if(sync.gt.smax) then + smax=sync + ibest=istart + dfbest=df + endif + enddo + enddo + f0=f0+dfbest +!f0=1443.75 + call ft4_downsample(iwave,f0,cb) ! downsample from 320s/Symbol to 20s/Symbol + sum2=sum(abs(cb)**2)/(20.0*76) + if(sum2.gt.0.0) cb=cb/sqrt(sum2) !ibest=208 - ib=ibest + cd=cb(ibest:ibest+76*20-1) + do k=1,NN + i1=(k-1)*20 + csymb=cd(i1:i1+19) + call four2a(csymb,20,1,-1,1) + cs(0:3,k)=csymb(1:4)/1e2 + s4(0:3,k)=abs(csymb(1:4)) + enddo - cd=cb(ib:ib+76*20-1) - do nseq=1,1 - if( nseq.eq.1 ) then ! noncoherent single-symbol detection - sbits1=0.0 - do isym=1,76 - ib=(isym-1)*20 - ccor(3,isym)=sum(cd(ib:ib+19)*conjg(c3(0:19))) - ccor(2,isym)=sum(cd(ib:ib+19)*conjg(c2(0:19))) - ccor(1,isym)=sum(cd(ib:ib+19)*conjg(c1(0:19))) - ccor(0,isym)=sum(cd(ib:ib+19)*conjg(c0(0:19))) - sbits1(2*isym-1)= max(abs(ccor(2,isym)),abs(ccor(3,isym)))- & - max(abs(ccor(0,isym)),abs(ccor(1,isym))) - sbits1(2*isym) = max(abs(ccor(1,isym)),abs(ccor(2,isym)))- & - max(abs(ccor(0,isym)),abs(ccor(3,isym))) - hbits1(2*isym-1:2*isym)=0 - if(sbits1(2*isym-1).gt.0) hbits1(2*isym-1)=1 - if(sbits1(2*isym ).gt.0) hbits1(2*isym )=1 - enddo - sbits=sbits1 - hbits=hbits1 - sbits3=sbits1 - hbits3=hbits1 - elseif( nseq.ge.2 ) then - nbit=2*nseq-1 - numseq=2**(nbit) - ps=0 - do ibit=nbit/2+1,144-nbit/2 - ps=0.0 - pmax=0.0 - do iseq=0,numseq-1 - csum=0.0 - cterm=1.0 - k=1 - do i=nbit-1,0,-1 - ibb=iand(iseq/(2**i),1) - csum=csum+ccor(ibb,ibit-(nbit/2+1)+k)*cterm - if(ibb.eq.0) cterm=cterm*cc0 - if(ibb.eq.1) cterm=cterm*cc1 - k=k+1 - enddo - ps(iseq)=abs(csum) - if( ps(iseq) .gt. pmax ) then - pmax=ps(iseq) - ibflag=1 - endif - enddo - if( ibflag .eq. 1 ) then - psbest=ps - ibflag=0 - endif - call getbitmetric(2**(nbit/2),psbest,numseq,sbits3(ibit)) - hbits3(ibit)=0 - if(sbits3(ibit).gt.0) hbits3(ibit)=1 - enddo - sbits=sbits3 - hbits=hbits3 - endif - nsync_qual=count(hbits(1:24).eq.(/0,0,0,0,0,0,1,0,1,0,1,0,1,0,1,0,1,0,0,0,0,0,0,0/)) -! if(nsync_qual.lt.10) exit - rxdata=sbits(25:152) - rxav=sum(rxdata(1:128))/128.0 - rx2av=sum(rxdata(1:128)*rxdata(1:128))/128.0 - rxsig=sqrt(rx2av-rxav*rxav) - rxdata=rxdata/rxsig - sigma=0.80 - llr(1:128)=2*rxdata/(sigma*sigma) - apmask=0 - max_iterations=40 - do ibias=0,0 - llr2=llr - if(ibias.eq.1) llr2=llr+0.4 - if(ibias.eq.2) llr2=llr-0.4 - call bpdecode128_90(llr2,apmask,max_iterations,message77,cw,nharderror,niterations) - if(nharderror.ge.0) exit - enddo - if(sum(message77).eq.0) cycle - if( nharderror.ge.0 ) then - write(c77,'(77i1)') message77(1:77) - call unpack77(c77,1,message,unpk77_success) - idupe=0 - do i=1,ndecodes - if(decodes(i).eq.message) idupe=1 - enddo - if(idupe.eq.1) goto 888 - ndecodes=ndecodes+1 - decodes(ndecodes)=message - nsnr=nint(xsnr) - freq=f0+dfbest -1210 format(a11,2i4,f6.2,f12.7,2x,a22,i3) - write(*,1212) datetime(8:11),nsnr,ibest/750.0,freq,message,'*',nseq,nharderror,nsync_qual -1212 format(a4,i4,2x,f5.3,f11.1,2x,a22,a1,i5,i5,i5) - goto 888 - endif - enddo ! nseq -888 continue - enddo !candidate list - enddo !files +! sync quality check + is1=0 + is2=0 + is3=0 + do k=1,4 + ip=maxloc(s4(:,k)) + if(icos4(k-1).eq.(ip(1)-1)) is1=is1+1 + ip=maxloc(s4(:,k+36)) + if(icos4(k-1).eq.(ip(1)-1)) is2=is2+1 + ip=maxloc(s4(:,k+72)) + if(icos4(k-1).eq.(ip(1)-1)) is3=is3+1 + enddo +! hard sync sum - max is 12 + nsync=is1+is2+is3 - write(*,1120) + do nseq=1,3 + if(nseq.eq.1) nsym=1 + if(nseq.eq.2) nsym=2 + if(nseq.eq.3) nsym=4 + nt=2**(2*nsym) + do ks=1,76,nsym + amax=-1.0 + do i=0,nt-1 + i1=i/64 + i2=iand(i,63)/16 + i3=iand(i,15)/4 + i4=iand(i,3) + if(nsym.eq.1) then + s2(i)=abs(cs(graymap(i4),ks)) + elseif(nsym.eq.2) then + s2(i)=abs(cs(graymap(i3),ks)+cs(graymap(i4),ks+1)) + elseif(nsym.eq.4) then + s2(i)=abs(cs(graymap(i1),ks ) + & + cs(graymap(i2),ks+1) + & + cs(graymap(i3),ks+2) + & + cs(graymap(i4),ks+3) & + ) + else + print*,"Error - nsym must be 1, 2, or 4." + endif + enddo + ipt=1+(ks-1)*2 + if(nsym.eq.1) ibmax=1 + if(nsym.eq.2) ibmax=3 + if(nsym.eq.4) ibmax=7 + do ib=0,ibmax + bm=maxval(s2(0:nt-1),one(0:nt-1,ibmax-ib)) - & + maxval(s2(0:nt-1),.not.one(0:nt-1,ibmax-ib)) + if(ipt+ib .gt.152) cycle + if(nsym.eq.1) then + bmeta(ipt+ib)=bm + elseif(nsym.eq.2) then + bmetb(ipt+ib)=bm + elseif(nsym.eq.4) then + bmetc(ipt+ib)=bm + endif + enddo + enddo + enddo + + call normalizebmet(bmeta,152) + call normalizebmet(bmetb,152) + call normalizebmet(bmetc,152) + + hbits=0 + where(bmeta.ge.0) hbits=1 + ns1=count(hbits( 1: 8).eq.(/0,0,0,1,1,0,1,1/)) + ns2=count(hbits( 73: 80).eq.(/0,0,0,1,1,0,1,1/)) + ns3=count(hbits(145:152).eq.(/0,0,0,1,1,0,1,1/)) + nsync_qual=ns1+ns2+ns3 + + sigma=0.7 + llra(1:64)=bmeta(9:72) + llra(65:128)=bmeta(81:144) + llra=2*llra/sigma**2 + llrb(1:64)=bmetb(9:72) + llrb(65:128)=bmetb(81:144) + llrb=2*llrb/sigma**2 + llrc(1:64)=bmetc(9:72) + llrc(65:128)=bmetc(81:144) + llrc=2*llrc/sigma**2 + + do isd=1,3 + if(isd.eq.1) llr=llra + if(isd.eq.2) llr=llrb + if(isd.eq.3) llr=llrc + apmask=0 + max_iterations=40 + do ibias=0,0 + llr2=llr + if(ibias.eq.1) llr2=llr+0.4 + if(ibias.eq.2) llr2=llr-0.4 + call bpdecode128_90(llr2,apmask,max_iterations,message77,cw,nharderror,niterations) + if(nharderror.ge.0) exit + enddo + if(sum(message77).eq.0) cycle + if( nharderror.ge.0 ) then + write(c77,'(77i1)') message77(1:77) + call unpack77(c77,1,message,unpk77_success) + idupe=0 + do i=1,ndecodes + if(decodes(i).eq.message) idupe=1 + enddo + if(idupe.eq.1) cycle + ndecodes=ndecodes+1 + decodes(ndecodes)=message + nsnr=nint(xsnr) + write(*,1212) datetime(8:11),nsnr,ibest/750.0,f0,message,'*',nharderror,nsync_qual,isd,niterations +1212 format(a4,i4,2x,f5.3,f11.1,2x,a22,a1,i5,i5,i5,i5) + endif + enddo ! sequence estimation + enddo !candidate list + enddo !files + + write(*,1120) 1120 format("") 999 end program ft4d subroutine getbitmetric(ib,ps,ns,xmet) - real ps(0:ns-1) - xm1=0 - xm0=0 - do i=0,ns-1 - if( iand(i/ib,1) .eq. 1 .and. ps(i) .gt. xm1 ) xm1=ps(i) - if( iand(i/ib,1) .eq. 0 .and. ps(i) .gt. xm0 ) xm0=ps(i) - enddo - xmet=xm1-xm0 - return + real ps(0:ns-1) + xm1=0 + xm0=0 + do i=0,ns-1 + if( iand(i/ib,1) .eq. 1 .and. ps(i) .gt. xm1 ) xm1=ps(i) + if( iand(i/ib,1) .eq. 0 .and. ps(i) .gt. xm0 ) xm0=ps(i) + enddo + xmet=xm1-xm0 + return end subroutine getbitmetric -subroutine downsample2(ci,f0,co) - parameter(NI=144*160,NH=NI/2,NO=NI/16) ! downsample from 200 samples per symbol to 10 - complex ci(0:NI-1),ct(0:NI-1) - complex co(0:NO-1) - fs=12000.0 - df=fs/NI - ct=ci - call four2a(ct,NI,1,-1,1) !c2c FFT to freq domain - i0=nint(f0/df) - ct=cshift(ct,i0) - co=0.0 - co(0)=ct(0) - b=8.0 - do i=1,NO/2 - arg=(i*df/b)**2 - filt=exp(-arg) - co(i)=ct(i)*filt - co(NO-i)=ct(NI-i)*filt - enddo - co=co/NO - call four2a(co,NO,1,1,1) !c2c FFT back to time domain - return -end subroutine downsample2 +subroutine downsample4(ci,f0,co) + parameter(NI=144*160,NH=NI/2,NO=NI/16) ! downsample from 200 samples per symbol to 10 + complex ci(0:NI-1),ct(0:NI-1) + complex co(0:NO-1) + fs=12000.0 + df=fs/NI + ct=ci + call four2a(ct,NI,1,-1,1) !c2c FFT to freq domain + i0=nint(f0/df) + ct=cshift(ct,i0) + co=0.0 + co(0)=ct(0) + b=8.0 + do i=1,NO/2 + arg=(i*df/b)**2 + filt=exp(-arg) + co(i)=ct(i)*filt + co(NO-i)=ct(NI-i)*filt + enddo + co=co/NO + call four2a(co,NO,1,1,1) !c2c FFT back to time domain + return +end subroutine downsample4 subroutine ft4_downsample(iwave,f0,c) ! Input: i*2 data in iwave() at sample rate 12000 Hz ! Output: Complex data in c(), sampled at 1200 Hz - include 'ft4_params.f90' - parameter (NFFT2=NMAX/16) - integer*2 iwave(NMAX) - complex c(0:NMAX/16-1) - complex c1(0:NFFT2-1) - complex cx(0:NMAX/2) - real x(NMAX) - equivalence (x,cx) + include 'ft4_params.f90' + parameter (NFFT2=NMAX/16) + integer*2 iwave(NMAX) + complex c(0:NMAX/16-1) + complex c1(0:NFFT2-1) + complex cx(0:NMAX/2) + real x(NMAX) + equivalence (x,cx) - BW=6.0*75 - df=12000.0/NMAX - x=iwave - call four2a(x,NMAX,1,-1,0) !r2c FFT to freq domain - ibw=nint(BW/df) - i0=nint(f0/df) - c1=0. - c1(0)=cx(i0) - do i=1,NFFT2/2 - arg=(i-1)*df/bw - win=exp(-arg*arg) - c1(i)=cx(i0+i)*win - c1(NFFT2-i)=cx(i0-i)*win - enddo - c1=c1/NFFT2 - call four2a(c1,NFFT2,1,1,1) !c2c FFT back to time domain - c=c1(0:NMAX/16-1) - return + BW=6.0*75 + df=12000.0/NMAX + x=iwave + call four2a(x,NMAX,1,-1,0) !r2c FFT to freq domain + ibw=nint(BW/df) + i0=nint(f0/df) + c1=0. + c1(0)=cx(i0) + do i=1,NFFT2/2 + arg=(i-1)*df/bw + win=exp(-arg*arg) + c1(i)=cx(i0+i)*win + c1(NFFT2-i)=cx(i0-i)*win + enddo + c1=c1/NFFT2 + call four2a(c1,NFFT2,1,1,1) !c2c FFT back to time domain + c=c1(0:NMAX/16-1) + return end subroutine ft4_downsample diff --git a/lib/fsk4hf/ft4_params.f90 b/lib/ft4/ft4_params.f90 similarity index 64% rename from lib/fsk4hf/ft4_params.f90 rename to lib/ft4/ft4_params.f90 index 3f3e9769d..07f7d66e9 100644 --- a/lib/fsk4hf/ft4_params.f90 +++ b/lib/ft4/ft4_params.f90 @@ -6,9 +6,9 @@ parameter (ND=64) !Data symbols parameter (NS=12) !Sync symbols (12) parameter (NN=NS+ND) !Total channel symbols (76) parameter (NSPS=320) !Samples per symbol at 12000 S/s -parameter (NZ=NSPS*NN) !Samples in full 1.92 s waveform (23040) -parameter (NMAX=2.5*12000) !Samples in iwave (36,000) -parameter (NFFT1=640, NH1=NFFT1/2) !Length of FFTs for symbol spectra -parameter (NSTEP=NSPS/4) !Rough time-sync step size +parameter (NZ=NSPS*NN) !Samples in full 2.03 s message frame (24320) +parameter (NMAX=2.5*12000) !Samples in iwave (30,000) +parameter (NFFT1=1280, NH1=NFFT1/2) !Length of FFTs for symbol spectra +parameter (NSTEP=NSPS/4) !Coarse time-sync step size parameter (NHSYM=NMAX/NSTEP-3) !Number of symbol spectra (1/4-sym steps) parameter (NDOWN=16) !Downsample factor diff --git a/lib/ft4/ft4d.f90 b/lib/ft4/ft4d.f90 new file mode 100644 index 000000000..6701e6b2e --- /dev/null +++ b/lib/ft4/ft4d.f90 @@ -0,0 +1,329 @@ +program ft4d + + use crc + use packjt77 + include 'ft4_params.f90' + character arg*8,message*37,c77*77,infile*80,fname*16,datetime*11 + character*37 decodes(100) + character*120 data_dir + character*90 dmsg + complex cd2(0:NMAX/16-1) !Complex waveform + complex cb(0:NMAX/16-1) + complex cd(0:76*20-1) !Complex waveform + complex csum,cterm + complex ctwk(80),ctwk2(80) + complex csymb(20) + complex cs(0:3,NN) + real s4(0:3,NN) + + real*8 fMHz + real ps(0:8191),psbest(0:8191) + real bmeta(152),bmetb(152),bmetc(152) + real s(NH1,NHSYM) + real a(5) + real llr(128),llr2(128),llra(128),llrb(128),llrc(128) + real s2(0:255) + real candidate(3,100) + real savg(NH1),sbase(NH1) + integer ihdr(11) + integer icos4(0:3) + integer*2 iwave(NMAX) !Generated full-length waveform + integer*1 message77(77),apmask(128),cw(128) + integer*1 hbits(152),hbits1(152),hbits3(152) + integer*1 s12(12) + integer graymap(0:3) + integer ip(1) + logical unpk77_success + logical one(0:511,0:7) ! 256 4-symbol sequences, 8 bits + data s12/1,1,1,2,2,2,2,2,2,1,1,1/ + data icos4/0,1,3,2/ + data graymap/0,1,3,2/ + save one + + fs=12000.0/NDOWN !Sample rate + dt=1/fs !Sample interval after downsample (s) + tt=NSPS*dt !Duration of "itone" symbols (s) + baud=1.0/tt !Keying rate for "itone" symbols (baud) + txt=NZ*dt !Transmission length (s) + twopi=8.0*atan(1.0) + h=1.0 !h=0.8 seems to be optimum for AWGN sensitivity (not for fading) + + one=.false. + do i=0,255 + do j=0,7 + if(iand(i,2**j).ne.0) one(i,j)=.true. + enddo + enddo + + nargs=iargc() + if(nargs.lt.1) then + print*,'Usage: ft4d [-a ] [-f fMHz] file1 [file2 ...]' + go to 999 + endif + iarg=1 + data_dir="." + call getarg(iarg,arg) + if(arg(1:2).eq.'-a') then + call getarg(iarg+1,data_dir) + iarg=iarg+2 + endif + call getarg(iarg,arg) + if(arg(1:2).eq.'-f') then + call getarg(iarg+1,arg) + read(arg,*) fMHz + iarg=iarg+2 + endif + ncoh=1 + + do ifile=iarg,nargs + call getarg(ifile,infile) + j2=index(infile,'.wav') + open(10,file=infile,status='old',access='stream') + read(10,end=999) ihdr,iwave + read(infile(j2-4:j2-1),*) nutc + datetime=infile(j2-11:j2-1) + close(10) + candidate=0.0 + ncand=0 + + nfqso=1500 + nfa=500 + nfb=2700 + syncmin=1.0 + maxcand=100 +! call syncft4(iwave,nfa,nfb,syncmin,nfqso,maxcand,s,candidate,ncand,sbase) + + call getcandidates4(iwave,375.0,3000.0,0.2,2200.0,100,savg,candidate,ncand,sbase) + ndecodes=0 + do icand=1,ncand + f0=candidate(1,icand)-1.5*37.5 + xsnr=1.0 + if( f0.le.375.0 .or. f0.ge.(5000.0-375.0) ) cycle + call ft4_downsample(iwave,f0,cd2) ! downsample from 320 Sa/Symbol to 20 Sa/Symbol + sum2=sum(cd2*conjg(cd2))/(20.0*76) + if(sum2.gt.0.0) cd2=cd2/sqrt(sum2) + +! 750 samples/second here + ibest=-1 + smax=-99. + dfbest=-1. + do idf=-90,+90,5 + df=idf + a=0. + a(1)=df + ctwk=1. + call twkfreq1(ctwk,80,fs,a,ctwk2) + do istart=0,315 + call sync4d(cd2,istart,ctwk2,1,sync) + if(sync.gt.smax) then + smax=sync + ibest=istart + dfbest=df + endif + enddo + enddo + + f0=f0+dfbest +f0=1443.75 + call ft4_downsample(iwave,f0,cb) ! downsample from 320s/Symbol to 20s/Symbol + sum2=sum(abs(cb)**2)/(20.0*76) + if(sum2.gt.0.0) cb=cb/sqrt(sum2) +ibest=208 + cd=cb(ibest:ibest+76*20-1) + do k=1,NN + i1=(k-1)*20 + csymb=cd(i1:i1+19) + call four2a(csymb,20,1,-1,1) + cs(0:3,k)=csymb(1:4)/1e2 + s4(0:3,k)=abs(csymb(1:4)) + enddo + +! sync quality check + is1=0 + is2=0 + is3=0 + do k=1,4 + ip=maxloc(s4(:,k)) + if(icos4(k-1).eq.(ip(1)-1)) is1=is1+1 + ip=maxloc(s4(:,k+36)) + if(icos4(k-1).eq.(ip(1)-1)) is2=is2+1 + ip=maxloc(s4(:,k+72)) + if(icos4(k-1).eq.(ip(1)-1)) is3=is3+1 + enddo +! hard sync sum - max is 12 + nsync=is1+is2+is3 + + do nseq=1,3 + if(nseq.eq.1) nsym=1 + if(nseq.eq.2) nsym=2 + if(nseq.eq.3) nsym=4 + nt=2**(2*nsym) + do ks=1,76,nsym + amax=-1.0 + do i=0,nt-1 + i1=i/64 + i2=iand(i,63)/16 + i3=iand(i,15)/4 + i4=iand(i,3) + if(nsym.eq.1) then + s2(i)=abs(cs(graymap(i4),ks)) + elseif(nsym.eq.2) then + s2(i)=abs(cs(graymap(i3),ks)+cs(graymap(i4),ks+1)) + elseif(nsym.eq.4) then + s2(i)=abs(cs(graymap(i1),ks ) + & + cs(graymap(i2),ks+1) + & + cs(graymap(i3),ks+2) + & + cs(graymap(i4),ks+3) & + ) + else + print*,"Error - nsym must be 1, 2, or 4." + endif + enddo + ipt=1+(ks-1)*2 + if(nsym.eq.1) ibmax=1 + if(nsym.eq.2) ibmax=3 + if(nsym.eq.4) ibmax=7 + do ib=0,ibmax + bm=maxval(s2(0:nt-1),one(0:nt-1,ibmax-ib)) - & + maxval(s2(0:nt-1),.not.one(0:nt-1,ibmax-ib)) + if(ipt+ib .gt.152) cycle + if(nsym.eq.1) then + bmeta(ipt+ib)=bm + elseif(nsym.eq.2) then + bmetb(ipt+ib)=bm + elseif(nsym.eq.4) then + bmetc(ipt+ib)=bm + endif + enddo + enddo + enddo + + call normalizebmet(bmeta,152) + call normalizebmet(bmetb,152) + call normalizebmet(bmetc,152) + + hbits=0 + where(bmeta.ge.0) hbits=1 + ns1=count(hbits( 1: 8).eq.(/0,0,0,1,1,0,1,1/)) + ns2=count(hbits( 73: 80).eq.(/0,0,0,1,1,0,1,1/)) + ns3=count(hbits(145:152).eq.(/0,0,0,1,1,0,1,1/)) + nsync_qual=ns1+ns2+ns3 + + sigma=0.7 + llra(1:64)=bmeta(9:72) + llra(65:128)=bmeta(81:144) + llra=2*llra/sigma**2 + llrb(1:64)=bmetb(9:72) + llrb(65:128)=bmetb(81:144) + llrb=2*llrb/sigma**2 + llrc(1:64)=bmetc(9:72) + llrc(65:128)=bmetc(81:144) + llrc=2*llrc/sigma**2 + + do isd=1,3 + if(isd.eq.1) llr=llra + if(isd.eq.2) llr=llrb + if(isd.eq.3) llr=llrc + apmask=0 + max_iterations=40 + do ibias=0,0 + llr2=llr + if(ibias.eq.1) llr2=llr+0.4 + if(ibias.eq.2) llr2=llr-0.4 + call bpdecode128_90(llr2,apmask,max_iterations,message77,cw,nharderror,niterations) + if(nharderror.ge.0) exit + enddo + if(sum(message77).eq.0) cycle + if( nharderror.ge.0 ) then + write(c77,'(77i1)') message77(1:77) + call unpack77(c77,1,message,unpk77_success) + idupe=0 + do i=1,ndecodes + if(decodes(i).eq.message) idupe=1 + enddo + if(idupe.eq.1) cycle + ndecodes=ndecodes+1 + decodes(ndecodes)=message + nsnr=nint(xsnr) + write(*,1212) datetime(8:11),nsnr,ibest/750.0,f0,message,'*',nharderror,nsync_qual,isd,niterations +1212 format(a4,i4,2x,f5.3,f11.1,2x,a22,a1,i5,i5,i5,i5) + endif + enddo ! sequence estimation + enddo !candidate list + enddo !files + + write(*,1120) +1120 format("") + +999 end program ft4d + +subroutine getbitmetric(ib,ps,ns,xmet) + real ps(0:ns-1) + xm1=0 + xm0=0 + do i=0,ns-1 + if( iand(i/ib,1) .eq. 1 .and. ps(i) .gt. xm1 ) xm1=ps(i) + if( iand(i/ib,1) .eq. 0 .and. ps(i) .gt. xm0 ) xm0=ps(i) + enddo + xmet=xm1-xm0 + return +end subroutine getbitmetric + +subroutine downsample4(ci,f0,co) + parameter(NI=144*160,NH=NI/2,NO=NI/16) ! downsample from 200 samples per symbol to 10 + complex ci(0:NI-1),ct(0:NI-1) + complex co(0:NO-1) + fs=12000.0 + df=fs/NI + ct=ci + call four2a(ct,NI,1,-1,1) !c2c FFT to freq domain + i0=nint(f0/df) + ct=cshift(ct,i0) + co=0.0 + co(0)=ct(0) + b=8.0 + do i=1,NO/2 + arg=(i*df/b)**2 + filt=exp(-arg) + co(i)=ct(i)*filt + co(NO-i)=ct(NI-i)*filt + enddo + co=co/NO + call four2a(co,NO,1,1,1) !c2c FFT back to time domain + return +end subroutine downsample4 + +subroutine ft4_downsample(iwave,f0,c) + +! Input: i*2 data in iwave() at sample rate 12000 Hz +! Output: Complex data in c(), sampled at 1200 Hz + + include 'ft4_params.f90' + parameter (NFFT2=NMAX/16) + integer*2 iwave(NMAX) + complex c(0:NMAX/16-1) + complex c1(0:NFFT2-1) + complex cx(0:NMAX/2) + real x(NMAX) + equivalence (x,cx) + + BW=6.0*75 + df=12000.0/NMAX + x=iwave + call four2a(x,NMAX,1,-1,0) !r2c FFT to freq domain + ibw=nint(BW/df) + i0=nint(f0/df) + c1=0. + c1(0)=cx(i0) + do i=1,NFFT2/2 + arg=(i-1)*df/bw + win=exp(-arg*arg) + c1(i)=cx(i0+i)*win + c1(NFFT2-i)=cx(i0-i)*win + enddo + c1=c1/NFFT2 + call four2a(c1,NFFT2,1,1,1) !c2c FFT back to time domain + c=c1(0:NMAX/16-1) + return +end subroutine ft4_downsample + diff --git a/lib/fsk4hf/ft4sim.f90 b/lib/ft4/ft4sim.f90 similarity index 96% rename from lib/fsk4hf/ft4sim.f90 rename to lib/ft4/ft4sim.f90 index c5fb3e572..92b09bdd0 100644 --- a/lib/fsk4hf/ft4sim.f90 +++ b/lib/ft4/ft4sim.f90 @@ -18,6 +18,8 @@ program ft4sim integer itone(NN) integer*1 msgbits(77) integer*2 iwave(NMAX) !Generated full-length waveform + integer icos4(4) + data icos4/0,1,3,2/ ! Get command-line argument(s) nargs=iargc() @@ -48,7 +50,7 @@ program ft4sim twopi=8.0*atan(1.0) fs=12000.0 !Sample rate (Hz) dt=1.0/fs !Sample interval (s) - hmod=1.000 !Modulation index (0.5 is MSK, 1.0 is FSK) + hmod=1.0 !Modulation index (0.5 is MSK, 1.0 is FSK) tt=NSPS*dt !Duration of symbols (s) baud=1.0/tt !Keying rate (baud) txt=NZ*dt !Transmission length (s) @@ -114,10 +116,6 @@ program ft4sim c0=cshift(c0,-k) ia=k -do i=0,NMAX-1 -write(21,*) i,real(c0(i)),imag(c0(i)),dphi(i) -enddo - do ifile=1,nfiles c=c0 if(fspread.ne.0.0 .or. delay.ne.0.0) call watterson(c,NMAX,NWAVE,fs,delay,fspread) diff --git a/lib/fsk4hf/genft4.f90 b/lib/ft4/genft4.f90 similarity index 82% rename from lib/fsk4hf/genft4.f90 rename to lib/ft4/genft4.f90 index 3aed650e5..450efcd4d 100644 --- a/lib/fsk4hf/genft4.f90 +++ b/lib/ft4/genft4.f90 @@ -15,12 +15,15 @@ subroutine genft4(msg0,ichk,msgsent,i4tone) character*37 message !Message to be generated character*37 msgsent !Message as it will be received character*77 c77 - integer*4 i4tone(76) + integer*4 i4tone(76),itmp(64) integer*1 codeword(128) integer*1 msgbits(77) integer*1 s12(12) + integer icos4(4) real*8 xi(864),xq(864),pi,twopi - data s12/0,0,0,3,3,3,3,3,3,0,0,0/ +! data s12/1,1,1,2,2,2,2,2,2,1,1,1/ + data icos4/0,1,3,2/ + logical unpk77_success twopi=8.*atan(1.0) @@ -56,13 +59,17 @@ subroutine genft4(msg0,ichk,msgsent,i4tone) ! 10 3 !Create 144-bit channel vector: - i4tone(1:12)=s12 do i=1,64 is=codeword(2*i)+2*codeword(2*i-1) - if(is.le.1) i4tone(12+i)=is - if(is.eq.2) i4tone(12+i)=3 - if(is.eq.3) i4tone(12+i)=2 + if(is.le.1) itmp(i)=is + if(is.eq.2) itmp(i)=3 + if(is.eq.3) itmp(i)=2 enddo + i4tone(1:4)=icos4 + i4tone(5:36)=itmp(1:32) + i4tone(37:40)=icos4 + i4tone(41:72)=itmp(33:64) + i4tone(73:76)=icos4 999 return end subroutine genft4 diff --git a/lib/fsk4hf/getcandidates4.f90 b/lib/ft4/getcandidates4.f90 similarity index 83% rename from lib/fsk4hf/getcandidates4.f90 rename to lib/ft4/getcandidates4.f90 index 428b8faf5..9d75187fe 100644 --- a/lib/fsk4hf/getcandidates4.f90 +++ b/lib/ft4/getcandidates4.f90 @@ -1,9 +1,7 @@ subroutine getcandidates4(id,fa,fb,syncmin,nfqso,maxcand,savg,candidate, & ncand,sbase) -! For now, hardwired to find the largest peak in the average spectrum - - include 'ft2_params.f90' + include 'ft4_params.f90' real s(NH1,NHSYM) real savg(NH1),savsm(NH1) real sbase(NH1) @@ -11,9 +9,7 @@ subroutine getcandidates4(id,fa,fb,syncmin,nfqso,maxcand,savg,candidate, & complex cx(0:NH1) real candidate(3,maxcand) integer*2 id(NMAX) - integer*1 s8(8) integer indx(NH1) - data s8/0,1,1,1,0,0,1,0/ equivalence (x,cx) ! Compute symbol spectra, stepping by NSTEP steps. @@ -33,10 +29,9 @@ subroutine getcandidates4(id,fa,fb,syncmin,nfqso,maxcand,savg,candidate, & savg=savg + s(1:NH1,j) !Average spectrum enddo savsm=0. - do i=2,NH1-1 - savsm(i)=sum(savg(i-1:i+1))/3. + do i=6,NH1-5 + savsm(i)=sum(savg(i-5:i+5))/11. enddo - nfa=fa/df nfb=fb/df np=nfb-nfa+1 @@ -58,6 +53,8 @@ subroutine getcandidates4(id,fa,fb,syncmin,nfqso,maxcand,savg,candidate, & if(xmax.gt.1.2) then ncand=ncand+1 candidate(1,ncand)=f0 + candidate(2,ncand)=-99.9 + candidate(3,ncand)=xmax endif return end subroutine getcandidates4 diff --git a/lib/ft4/sync4d.f90 b/lib/ft4/sync4d.f90 new file mode 100644 index 000000000..efbe9b94f --- /dev/null +++ b/lib/ft4/sync4d.f90 @@ -0,0 +1,54 @@ +subroutine sync4d(cd0,i0,ctwk,itwk,sync) + +! Compute sync power for a complex, downsampled FT4 signal. +! 20 samples per symbol + + include 'ft4_params.f90' + parameter(NP=NMAX/16) + complex cd0(0:NP-1) + complex csync(80) + complex csync2(80) + complex ctwk(80) + complex z1,z2,z3 + logical first + integer icos4(0:3) + data icos4/0,1,3,2/ + data first/.true./ + save first,twopi,fs2,dt2,taus,baud,csync + + p(z1)=real(z1)**2 + aimag(z1)**2 !Statement function for power + +! Set some constants and compute the csync array. + if( first ) then + twopi=8.0*atan(1.0) + fs2=12000.0/NDOWN !Sample rate after downsampling + dt2=1/fs2 !Corresponding sample interval + taus=20*dt2 !Symbol duration + baud=1.0/taus !Keying rate + k=1 + phi=0.0 + do i=0,3 +! dphi=(twopi/2.0)*(2*icos4(i)-3)*baud*dt2 + dphi=twopi*icos4(i)*baud*dt2 + do j=1,20 + csync(k)=cmplx(cos(phi),sin(phi)) !Waveform for 7x7 Costas array + phi=mod(phi+dphi,twopi) + k=k+1 + enddo + enddo + first=.false. + endif + + sync=0 + i1=i0 !three Costas arrays + i2=i0+36*20-1 + i3=i0+72*20-1 + csync2=csync + if(itwk.eq.1) csync2=ctwk*csync2 !Tweak the frequency + if(i1.ge.0 .and. i1+79.le.NP-1) z1=sum(cd0(i1:i1+79)*conjg(csync2)) + if(i2.ge.0 .and. i2+79.le.NP-1) z2=sum(cd0(i2:i2+79)*conjg(csync2)) + if(i3.ge.0 .and. i3+79.le.NP-1) z3=sum(cd0(i3:i3+79)*conjg(csync2)) + sync = sync + p(z1) + p(z2) + p(z3) + + return +end subroutine sync4d diff --git a/lib/ft4/syncft4.f90 b/lib/ft4/syncft4.f90 new file mode 100644 index 000000000..db21b1e1e --- /dev/null +++ b/lib/ft4/syncft4.f90 @@ -0,0 +1,145 @@ +subroutine syncft4(iwave,nfa,nfb,syncmin,nfqso,maxcand,s,candidate, & + ncand,sbase) + + include 'ft4_params.f90' +! Search over +/- 2.5s relative to 0.5s TX start time. + parameter (JZ=20) + complex cx(0:NH1) + real s(NH1,NHSYM) + real savg(NH1) + real sbase(NH1) + real x(NFFT1) + real sync2d(NH1,-JZ:JZ) + real red(NH1) + real candidate0(3,maxcand) + real candidate(3,maxcand) + real dd(NMAX) + integer jpeak(NH1) + integer indx(NH1) + integer ii(1) + integer*2 iwave(NMAX) + integer icos4(0:3) + data icos4/0,1,3,2/ !Costas 4x4 tone pattern + equivalence (x,cx) + + dd=iwave/1e3 +! Compute symbol spectra, stepping by NSTEP steps. + savg=0. + tstep=NSTEP/12000.0 + df=12000.0/NFFT1 !3.125 Hz + fac=1.0/300.0 + do j=1,NHSYM + ia=(j-1)*NSTEP + 1 + ib=ia+NSPS-1 + x(1:NSPS)=fac*dd(ia:ib) + x(NSPS+1:)=0. + call four2a(x,NFFT1,1,-1,0) !r2c FFT + do i=1,NH1 + s(i,j)=real(cx(i))**2 + aimag(cx(i))**2 + enddo + savg=savg + s(1:NH1,j) !Average spectrum + enddo + + call baseline(savg,nfa,nfb,sbase) + + ia=max(1,nint(nfa/df)) + ib=nint(nfb/df) + nssy=NSPS/NSTEP ! # steps per symbol + nfos=NFFT1/NSPS ! # frequency bin oversampling factor + jstrt=0.25/tstep + candidate0=0. + k=0 + + do i=ia,ib + do j=-JZ,+JZ + ta=0. + tb=0. + tc=0. + t0a=0. + t0b=0. + t0c=0. + do n=0,3 + m=j+jstrt+nssy*n + if(m.ge.1.and.m.le.NHSYM) then + ta=ta + s(i+nfos*icos4(n),m) + t0a=t0a + sum(s(i:i+nfos*3:nfos,m)) + endif + tb=tb + s(i+nfos*icos4(n),m+nssy*36) + t0b=t0b + sum(s(i:i+nfos*3:nfos,m+nssy*36)) + if(m+nssy*72.le.NHSYM) then + tc=tc + s(i+nfos*icos4(n),m+nssy*72) + t0c=t0c + sum(s(i:i+nfos*3:nfos,m+nssy*72)) + endif + enddo + t=ta+tb+tc + t0=t0a+t0b+t0c + t0=(t0-t)/3.0 + sync_abc=t/t0 + t=tb+tc + t0=t0b+t0c + t0=(t0-t)/3.0 + sync_bc=t/t0 + sync2d(i,j)=max(sync_abc,sync_bc) + enddo + enddo + + red=0. + do i=ia,ib + ii=maxloc(sync2d(i,-JZ:JZ)) - 1 - JZ + j0=ii(1) + jpeak(i)=j0 + red(i)=sync2d(i,j0) + enddo + iz=ib-ia+1 + call indexx(red(ia:ib),iz,indx) + ibase=indx(nint(0.40*iz)) - 1 + ia + if(ibase.lt.1) ibase=1 + if(ibase.gt.nh1) ibase=nh1 + base=red(ibase) + red=red/base + do i=1,min(maxcand,iz) + n=ia + indx(iz+1-i) - 1 + if(red(n).lt.syncmin.or.isnan(red(n)).or.k.eq.maxcand) exit + k=k+1 +! candidate0(1,k)=n*df+37.5*1.5 + candidate0(1,k)=n*df + candidate0(2,k)=(jpeak(n)-1)*tstep + candidate0(3,k)=red(n) + enddo + ncand=k + +! Put nfqso at top of list, and save only the best of near-dupe freqs. + do i=1,ncand + if(abs(candidate0(1,i)-nfqso).lt.10.0) candidate0(1,i)=-candidate0(1,i) + if(i.ge.2) then + do j=1,i-1 + fdiff=abs(candidate0(1,i))-abs(candidate0(1,j)) + if(abs(fdiff).lt.4.0) then + if(candidate0(3,i).ge.candidate0(3,j)) candidate0(3,j)=0. + if(candidate0(3,i).lt.candidate0(3,j)) candidate0(3,i)=0. + endif + enddo + endif + enddo + + fac=20.0/maxval(s) + s=fac*s + +! Sort by sync +! call indexx(candidate0(3,1:ncand),ncand,indx) +! Sort by frequency + call indexx(candidate0(1,1:ncand),ncand,indx) + k=1 +! do i=ncand,1,-1 + do i=1,ncand + j=indx(i) +! if( candidate0(3,j) .ge. syncmin .and. candidate0(2,j).ge.-1.5 ) then + if( candidate0(3,j) .ge. syncmin ) then + candidate(2:3,k)=candidate0(2:3,j) + candidate(1,k)=abs(candidate0(1,j)) + k=k+1 + endif + enddo + ncand=k-1 + return +end subroutine syncft4