WSJT-X/lib/ft4/subtractft4.f90
Steve Franke a9623703b3 FT4: Change Fortran code to NSPS=576 and make necessary accompanying
changes on the C++ side. Basically works except that Tx audio has
incorrect DT and audio is truncated at the end. Also, command line
decoding using JT9 is not as sensitive as decoding from within WSJT-X.
2019-05-22 17:02:15 -05:00

67 lines
1.6 KiB
Fortran

subroutine subtractft4(dd,itone,f0,dt)
! Subtract an ft4 signal
!
! Measured signal : dd(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
parameter (NMAX=21*3456,NSPS=576,NFFT=NMAX,NFILT=1400)
parameter (NFRAME=(103+2)*NSPS)
real*4 dd(NMAX), window(-NFILT/2:NFILT/2), xjunk
complex cref,camp,cfilt,cw
integer itone(103)
logical first
data first/.true./
common/heap8/cref(NFRAME),camp(NMAX),cfilt(NMAX),cw(NMAX),xjunk(NFRAME)
save first
nstart=dt*12000+1-NSPS
nsym=103
fs=12000.0
icmplx=1
bt=1.0
nss=NSPS
call gen_ft4wave(itone,nsym,nss,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
pi=4.0*atan(1.0)
fac=1.0/float(nfft)
sum=0.0
do j=-NFILT/2,NFILT/2
window(j)=cos(pi*j/NFILT)**2
sum=sum+window(j)
enddo
cw=0.
cw(1:NFILT+1)=window/sum
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)
! 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
return
end subroutine subtractft4