From ff46c5a0c42dd843e3196ddb57b2282348ae75bd Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Mon, 24 Feb 2020 16:18:42 -0500 Subject: [PATCH] In subtractft8.f90: refine DT for early decodes before subtracting them from dd(). --- lib/ft8/ft8b.f90 | 2 +- lib/ft8/subtractft8.f90 | 121 ++++++++++++++++++++++++++++------------ lib/ft8_decode.f90 | 2 +- 3 files changed, 86 insertions(+), 39 deletions(-) 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..cf485b657 100644 --- a/lib/ft8/subtractft8.f90 +++ b/lib/ft8/subtractft8.f90 @@ -1,66 +1,113 @@ -subroutine subtractft8(dd,itone,f0,dt) +subroutine subtractft8(dd0,itone,f0,dt,ldt) -! Subtract an ft8 signal -! -! Measured signal : dd(t) = a(t)cos(2*pi*f0*t+theta(t)) +! Subtract an ft8 signal. If ldt==.true., refine DT first. + +! Raw data : dd(t) = a(t)cos(2*pi*f0*t+theta(t)) +! Reference signal : cref(t) = exp( j*(2*pi*f0*t+phi(t)) ) + + parameter (NMAX=15*12000,NFRAME=1920*79) + real dd(NMAX),dd0(NMAX) + complex cref(NFRAME) + logical ldt + +! Generate complex reference waveform + call gen_ft8wave(itone,79,1920,2.0,12000.0,f0,cref,xjunk,1,NFRAME) + + if(ldt) then !Are we refining DT ? + sqa=sqf(dd0,cref,f0,dt,ldt,-300,dd) !Yes + sqb=sqf(dd0,cref,f0,dt,ldt,300,dd) + endif + sq0=sqf(dd0,cref,f0,dt,ldt,0,dd) !Do the subtraction with idt=0 + if(ldt) then + call peakup(sqa,sq0,sqb,dx) + if(abs(dx).gt.1.0) goto 100 !No acceptable minimum: do not subtract + i1=nint(300.0*dx) !First approximation of better idt + sqa=sqf(dd0,cref,f0,dt,ldt,i1-60,dd) + sqb=sqf(dd0,cref,f0,dt,ldt,i1+60,dd) + sq0=sqf(dd0,cref,f0,dt,ldt,i1,dd) + call peakup(sqa,sq0,sqb,dx) + if(abs(dx).gt.1.0) then !No acceptable minimum + sq0=sqf(dd0,cref,f0,dt,ldt,0,dd) !Use idt=0 for subtraction + go to 100 + endif + i2=nint(60.0*dx) + i1 !Best estimate of idt + sq0=sqf(dd0,cref,f0,dt,ldt,i2,dd) !Do the subtraction with idt=i2 + endif +100 dd0=dd !Return dd0 with signal subtracted + + return +end subroutine subtractft8 + +real function sqf(dd0,cref,f0,dt,ldt,idt,dd) + +! Raw data : dd0(t) = a(t)cos(2*pi*f0*t+theta(t)) ! Reference signal : cref(t) = exp( j*(2*pi*f0*t+phi(t)) ) ! Complex amp : cfilt(t) = LPF[ dd(t)*CONJG(cref(t)) ] -! Subtract : dd(t) = dd(t) - 2*REAL{cref*cfilt} - - use timer_module, only: timer - +! Subtract : dd(t) = dd0(t) - 2*REAL{cref*cfilt} + 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 - integer itone(79) - logical first + 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),cref(NFRAME) + complex camp,cfilt,cw,z + logical first,ldt data first/.true./ - common/heap8/cref(NFRAME),camp(NMAX),cfilt(NMAX),cw(NMAX),xjunk(NFRAME) - save first - - 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 + common/heap8/camp(NMAX),cfilt(NMAX),cw(NMAX) + equivalence (x,cx) + save first,/heap8/ 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 + + 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) -! Subtract the reconstructed signal + x=0. 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)) + 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 + sq=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 + sq=sq + real(cx(i))*real(cx(i)) + aimag(cx(i))*aimag(cx(i)) + enddo + endif + sqf=sq return -end subroutine subtractft8 - +end function sqf 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