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