1subroutine subtractft8(dd0,itone,f0,dt,lrefinedt) 2 3! Subtract an ft8 signal 4! 5! Measured signal : dd(t) = a(t)cos(2*pi*f0*t+theta(t)) 6! Reference signal : cref(t) = exp( j*(2*pi*f0*t+phi(t)) ) 7! Complex amp : cfilt(t) = LPF[ dd(t)*CONJG(cref(t)) ] 8! Subtract : dd(t) = dd(t) - 2*REAL{cref*cfilt} 9 10 parameter (NMAX=15*12000,NFRAME=1920*79) 11 parameter (NFFT=NMAX,NFILT=4000) 12 real dd(NMAX),dd0(NMAX) 13 real window(-NFILT/2:NFILT/2) 14 real x(NFFT+2) 15 real endcorrection(NFILT/2+1) 16 complex cx(0:NFFT/2) 17 complex cref,camp,cfilt,cw,z 18 integer itone(79) 19 logical first,lrefinedt,ldt 20 data first/.true./ 21 common/heap8/cref(NFRAME),camp(NMAX),cfilt(NMAX),cw(NMAX) 22 equivalence (x,cx) 23 save first,/heap8/,endcorrection 24 25 if(first) then ! Create and normalize the filter 26 pi=4.0*atan(1.0) 27 fac=1.0/float(nfft) 28 sumw=0.0 29 do j=-NFILT/2,NFILT/2 30 window(j)=cos(pi*j/NFILT)**2 31 sumw=sumw+window(j) 32 enddo 33 cw=0. 34 cw(1:NFILT+1)=window/sumw 35 cw=cshift(cw,NFILT/2+1) 36 call four2a(cw,nfft,1,-1,1) 37 cw=cw*fac 38 first=.false. 39 do j=1,NFILT/2+1 40 endcorrection(j)=1.0/(1.0-sum(window(j-1:NFILT/2))/sumw) 41 enddo 42 endif 43 44! Generate complex reference waveform cref 45 call gen_ft8wave(itone,79,1920,2.0,12000.0,f0,cref,xjunk,1,NFRAME) 46 47 ldt=lrefinedt 48 if(ldt) then !Are we refining DT ? 49 sqa=sqf(-90) 50 sqb=sqf(+90) 51 sq0=sqf(0) 52 call peakup(sqa,sq0,sqb,dx) 53 if(abs(dx).gt.1.0) return !No acceptable minimum: do not subtract 54 i2=nint(90.0*dx) !Best estimate of idt 55 ldt=.false. 56 sq0=sqf(i2) !Do the subtraction with idt=i2 57 else 58 sq0=sqf(0) !Do the subtraction with idt=0 59 endif 60 dd0=dd !Return dd0 with this signal subtracted 61 return 62 63contains 64 65 real function sqf(idt) !Internal function: all variables accessible 66 nstart=dt*12000+1 + idt 67 camp=0. 68 dd=dd0 69 do i=1,nframe 70 j=nstart-1+i 71 if(j.ge.1.and.j.le.NMAX) camp(i)=dd(j)*conjg(cref(i)) 72 enddo 73 74 cfilt(1:nframe)=camp(1:nframe) 75 cfilt(nframe+1:)=0.0 76 call four2a(cfilt,nfft,1,-1,1) 77 cfilt(1:nfft)=cfilt(1:nfft)*cw(1:nfft) 78 call four2a(cfilt,nfft,1,1,1) 79 cfilt(1:NFILT/2+1)=cfilt(1:NFILT/2+1)*endcorrection 80 cfilt(nframe:nframe-NFILT/2:-1)=cfilt(nframe:nframe-NFILT/2:-1)*endcorrection 81 x=0. 82 do i=1,nframe 83 j=nstart+i-1 84 if(j.ge.1 .and. j.le.NMAX) then 85 z=cfilt(i)*cref(i) 86 dd(j)=dd(j)-2.0*real(z) !Subtract the reconstructed signal 87 x(i)=dd(j) 88 endif 89 enddo 90 sqq=0. 91 if(ldt) then 92 call four2a(cx,NFFT,1,-1,0) !Forward FFT, r2c 93 df=12000.0/NFFT 94 ia=(f0-1.5*6.25)/df 95 ib=(f0+8.5*6.25)/df 96 do i=ia,ib 97 sqq=sqq + real(cx(i))*real(cx(i)) + aimag(cx(i))*aimag(cx(i)) 98 enddo 99 endif 100 sqf=sqq 101 return 102 end function sqf 103 104end subroutine subtractft8 105