diff --git a/lib/extract.f90 b/lib/extract.f90 index d1a483d35..78da285f7 100644 --- a/lib/extract.f90 +++ b/lib/extract.f90 @@ -209,7 +209,7 @@ subroutine extract(s3,nadd,mode65,ntrials,naggressive,ndepth,nflip, & tmp(i)=correct(64-i) enddo correct(1:63)=tmp(1:63) - call interleave63(correct,63,1) + call interleave63(correct,1) call graycode65(correct,63,1) call unpackmsg(dat4,decoded) !Unpack the user message ncount=0 diff --git a/lib/ft8/filt8.f90 b/lib/ft8/filt8.f90 index aeacf2a74..abe797ba7 100644 --- a/lib/ft8/filt8.f90 +++ b/lib/ft8/filt8.f90 @@ -9,7 +9,7 @@ subroutine filt8(f0,nslots,width,wave) equivalence (x,cx) x=wave - call four2a(x,NFFT,1,-1,0) !r2c + call four2a(cx,NFFT,1,-1,0) !r2c df=12000.0/NFFT fa=f0 - 0.5*6.25 fb=f0 + 7.5*6.25 + (nslots-1)*60.0 diff --git a/lib/ft8/ft8b.f90 b/lib/ft8/ft8b.f90 index 70040d1af..ec575426c 100644 --- a/lib/ft8/ft8b.f90 +++ b/lib/ft8/ft8b.f90 @@ -418,7 +418,7 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,nzhsym,lapon, & endif nbadcrc=0 ! If we get this far: valid codeword, valid (i3,n3), nonquirky message. call get_ft8_tones_from_77bits(message77,itone) - if(lsubtract) call subtractft8(dd0,itone,f1,xdt) + if(lsubtract) call subtractft8(dd0,itone,f1,xdt,.false.) ! write(21,3001) nzhsym,npasses,nqsoprogress,ipass,iaptype,lsubtract, & ! f1,xdt,msg37(1:22); flush(21) !3001 format(5i3,L3,f7.1,f7.2,2x,a22) diff --git a/lib/ft8/subtractft8.f90 b/lib/ft8/subtractft8.f90 index ec0beea9a..00b526934 100644 --- a/lib/ft8/subtractft8.f90 +++ b/lib/ft8/subtractft8.f90 @@ -1,4 +1,4 @@ -subroutine subtractft8(dd,itone,f0,dt) +subroutine subtractft8(dd0,itone,f0,dt,ldt) ! Subtract an ft8 signal ! @@ -7,60 +7,98 @@ subroutine subtractft8(dd,itone,f0,dt) ! Complex amp : cfilt(t) = LPF[ dd(t)*CONJG(cref(t)) ] ! Subtract : dd(t) = dd(t) - 2*REAL{cref*cfilt} - use timer_module, only: timer - parameter (NMAX=15*12000,NFRAME=1920*79) - parameter (NFFT=NMAX,NFILT=1400) - real*4 dd(NMAX), window(-NFILT/2:NFILT/2), xjunk - complex cref,camp,cfilt,cw + parameter (NFFT=NMAX,NFILT=2800) + real dd(NMAX),dd0(NMAX) + real window(-NFILT/2:NFILT/2) + real x(NFFT+2) + complex cx(0:NFFT/2) + complex cref,camp,cfilt,cw,z integer itone(79) - logical first + logical first,ldt data first/.true./ - common/heap8/cref(NFRAME),camp(NMAX),cfilt(NMAX),cw(NMAX),xjunk(NFRAME) - save first + common/heap8/cref(NFRAME),camp(NMAX),cfilt(NMAX),cw(NMAX) + equivalence (x,cx) + save first,/heap8/ - nstart=dt*12000+1 - nsym=79 - nsps=1920 - fs=12000.0 - icmplx=1 - bt=2.0 - call gen_ft8wave(itone,nsym,nsps,bt,fs,f0,cref,xjunk,icmplx,NFRAME) - camp=0. - do i=1,nframe - id=nstart-1+i - if(id.ge.1.and.id.le.NMAX) camp(i)=dd(id)*conjg(cref(i)) - enddo - - if(first) then -! Create and normalize the filter + if(first) then ! Create and normalize the filter pi=4.0*atan(1.0) fac=1.0/float(nfft) - sum=0.0 + sumw=0.0 do j=-NFILT/2,NFILT/2 window(j)=cos(pi*j/NFILT)**2 - sum=sum+window(j) + sumw=sumw+window(j) enddo cw=0. - cw(1:NFILT+1)=window/sum + cw(1:NFILT+1)=window/sumw cw=cshift(cw,NFILT/2+1) call four2a(cw,nfft,1,-1,1) cw=cw*fac first=.false. endif - cfilt=0.0 - cfilt(1:nframe)=camp(1:nframe) - call four2a(cfilt,nfft,1,-1,1) - cfilt(1:nfft)=cfilt(1:nfft)*cw(1:nfft) - call four2a(cfilt,nfft,1,1,1) +! Generate complex reference waveform + call gen_ft8wave(itone,79,1920,2.0,12000.0,f0,cref,xjunk,1,NFRAME) -! Subtract the reconstructed signal - do i=1,nframe - j=nstart+i-1 - if(j.ge.1 .and. j.le.NMAX) dd(j)=dd(j)-2*REAL(cfilt(i)*cref(i)) - enddo + if(ldt) then !Are we refining DT ? + sqa=sqf(-300) + sqb=sqf(300) + endif + sq0=sqf(0) !Do the subtraction with idt=0 + if(ldt) then + call peakup(sqa,sq0,sqb,dx) + if(abs(dx).gt.1.0) return !No acceptable minimum: do not subtract + i1=nint(300.0*dx) !First approximation of best idt + sqa=sqf(i1-60) + sqb=sqf(i1+60) + sq0=sqf(i1) + call peakup(sqa,sq0,sqb,dx) + if(abs(dx).gt.1.0) return !No acceptable minimum: do not subtract + i2=nint(60.0*dx) + i1 !Best estimate of idt + sq0=sqf(i2) !Do the subtraction with idt=i2 + endif + dd0=dd !Return dd0 with this signal subtracted return -end subroutine subtractft8 +contains + + real function sqf(idt) !Internal function: all variables accessible + nstart=dt*12000+1 + idt + camp=0. + dd=dd0 + do i=1,nframe + j=nstart-1+i + if(j.ge.1.and.j.le.NMAX) camp(i)=dd(j)*conjg(cref(i)) + enddo + + cfilt(1:nframe)=camp(1:nframe) + cfilt(nframe+1:)=0.0 + call four2a(cfilt,nfft,1,-1,1) + cfilt(1:nfft)=cfilt(1:nfft)*cw(1:nfft) + call four2a(cfilt,nfft,1,1,1) + + x=0. + do i=1,nframe + j=nstart+i-1 + if(j.ge.1 .and. j.le.NMAX) then + z=cfilt(i)*cref(i) + dd(j)=dd(j)-2.0*real(z) !Subtract the reconstructed signal + x(i)=dd(j) + endif + enddo + sqq=0. + if(ldt) then + call four2a(cx,NFFT,1,-1,0) !Forward FFT, r2c + df=12000.0/NFFT + ia=(f0-1.5*6.25)/df + ib=(f0+8.5*6.25)/df + do i=ia,ib + sqq=sqq + real(cx(i))*real(cx(i)) + aimag(cx(i))*aimag(cx(i)) + enddo + endif + sqf=sqq + return + end function sqf + +end subroutine subtractft8 diff --git a/lib/ft8_decode.f90 b/lib/ft8_decode.f90 index 64fb46ee2..d43f5948f 100644 --- a/lib/ft8_decode.f90 +++ b/lib/ft8_decode.f90 @@ -77,7 +77,7 @@ contains endif if(nzhsym.eq.50 .and. ndec_early.ge.1) then do i=1,ndec_early - call subtractft8(dd,itone_save(1,i),f1_save(i),xdt_save(i)) + call subtractft8(dd,itone_save(1,i),f1_save(i),xdt_save(i),.true.) enddo endif ifa=nfa diff --git a/lib/lpf1.f90 b/lib/lpf1.f90 index f2bb2377c..a620e9a55 100644 --- a/lib/lpf1.f90 +++ b/lib/lpf1.f90 @@ -11,7 +11,7 @@ subroutine lpf1(dd,jz,dat,jz2) fac=1.0/float(NFFT1) x(1:jz)=fac*dd(1:jz) x(jz+1:NFFT1)=0.0 - call four2a(x,NFFT1,1,-1,0) !Forwarxd FFT, r2c + call four2a(cx,NFFT1,1,-1,0) !Forwarxd FFT, r2c cx(NFFT2/2:)=0.0 ! df=11025.0/NFFT1 diff --git a/lib/refspectrum.f90 b/lib/refspectrum.f90 index bc5bede33..fdd302ae5 100644 --- a/lib/refspectrum.f90 +++ b/lib/refspectrum.f90 @@ -40,7 +40,7 @@ subroutine refspectrum(id2,bclear,brefspec,buseref,fname) if(brefspec) then x(0:NH-1)=0.001*id2(1:NH) x(NH:NFFT-1)=0.0 - call four2a(x,NFFT,1,-1,0) !r2c FFT + call four2a(cx,NFFT,1,-1,0) !r2c FFT do i=1,NH s(i)=s(i) + real(cx(i))**2 + aimag(cx(i))**2 @@ -134,7 +134,7 @@ subroutine refspectrum(id2,bclear,brefspec,buseref,fname) ! Make the filter causal for overlap and add. cx(0)=0.0 cx(1:NH)=fil(1:NH)/NFFT - call four2a(x,NFFT,1,1,-1) + call four2a(cx,NFFT,1,1,-1) x=cshift(x,-400) x(800:NH)=0.0 call four2a(cx,NFFT,1,-1,0) @@ -146,7 +146,7 @@ subroutine refspectrum(id2,bclear,brefspec,buseref,fname) x(0:NH-1)=id2(1:NH) x(NH:NFFT-1)=0.0 x=x/NFFT - call four2a(x,NFFT,1,-1,0) + call four2a(cx,NFFT,1,-1,0) cx=cfil*cx call four2a(cx,NFFT,1,1,-1) x(0:NH-1)=x(0:NH-1)+xs diff --git a/lib/timf2.f90 b/lib/timf2.f90 index 07f509831..bec391e2e 100644 --- a/lib/timf2.f90 +++ b/lib/timf2.f90 @@ -59,7 +59,7 @@ subroutine timf2(x0,k,nfft,nwindow,nb,peaklimit,x1, & x(0:nfft-1)=x0 if(nwindow.eq.2) x(0:nfft-1)=w(0:nfft-1)*x(0:nfft-1) - call four2a(x,nfft,1,-1,0) !First forward FFT, r2c + call four2a(cx,nfft,1,-1,0) !First forward FFT, r2c cxt(0:nh)=cx(0:nh) ! Identify frequencies with strong signals. diff --git a/lib/wav11.f90 b/lib/wav11.f90 index 3bba0bf8a..219b7880a 100644 --- a/lib/wav11.f90 +++ b/lib/wav11.f90 @@ -14,7 +14,7 @@ subroutine wav11(d2,npts,dd) jz=min(NZ12,npts) x(1:jz)=d2(1:jz) x(jz+1:)=0.0 - call four2a(x,nfft1,1,-1,0) !Forward FFT, r2c + call four2a(cx,nfft1,1,-1,0) !Forward FFT, r2c df=12000.0/NFFT1 ia=5000.0/df cx(ia:)=0.0 diff --git a/lib/wav12.f90 b/lib/wav12.f90 index 31fe1d599..ff655436a 100644 --- a/lib/wav12.f90 +++ b/lib/wav12.f90 @@ -34,7 +34,7 @@ subroutine wav12(d2,d1,npts,nbitsam2) x(1:jz)=d2(1:jz) x(jz+1:)=0.0 - call four2a(x,nfft1,1,-1,0) !Forwarxd FFT, r2c + call four2a(cx,nfft1,1,-1,0) !Forwarxd FFT, r2c cx(nfft1/2:)=0.0 call four2a(cx,nfft2,1,1,-1) !Inverse FFT, c2r