WSJT-X/lib/analytic.f90

subroutine analytic(d,npts,nfft,c,dpc,bseq,bdeq)

! Convert real data to analytic signal

  parameter (NFFTMAX=1024*1024)

  real d(npts)              ! passband signal
  real h(NFFTMAX/2)         ! real BPF magnitude
  real dpc(3),dpclast(3)    ! dynamic phase coeffs
  real spc(3),spclast(3)    ! static phase coeffs
  real sac(5),saclast(5)    ! amp coeffs
  real fp     

  complex corrs(NFFTMAX/2)  ! allpass static phase correction 
  complex corrd(NFFTMAX/2)   ! allpass overall phase correction
  complex c(NFFTMAX)        ! analytic signal

  logical*1 bseq            ! boolean static equalizer flag
  logical*1 bdeq            ! boolean dynamic equalizer flag
  logical*1 bseqlast      

  data nfft0/0/
  data bseqlast/.false./
  data spclast/0.0,0.0,0.0/
  data saclast/1.0,0.0,0.0,0.0,0.0/
  data spc/-0.952,0.768,-0.565/  ! baseline phase coeffs for TS2000
  data sac/1.0,0.05532,0.11438,0.12918,0.09274/ ! amp coeffs for TS2000

  save corrs,corrd,nfft0,h,sac,spc,saclast,spclast,dpclast,pi,t,beta

  df=12000.0/nfft
  nh=nfft/2
  if( nfft.ne.nfft0 ) then
     pi=4.0*atan(1.0)
     t=1.0/2000.0
     beta=0.1
     do i=1,nh+1
        ff=(i-1)*df
        f=ff-1500.0
        h(i)=1.0
        if(abs(f).gt.(1-beta)/(2*t) .and. abs(f).le.(1+beta)/(2*t)) then
           h(i)=h(i)*0.5*(1+cos((pi*t/beta )*(abs(f)-(1-beta)/(2*t))))
        elseif( abs(f) .gt. (1+beta)/(2*t) ) then
           h(i)=0.0
        endif
     enddo
     nfft0=nfft
  endif

  if( any(spclast .ne. spc) .or. any(saclast .ne. sac) .or. any(dpclast .ne. dpc) ) then
     spclast=spc
     dpclast=dpc
     saclast=sac
     do i=1,nh+1
        ff=(i-1)*df
        f=ff-1500.0
        fp=f/1000.0
        ps=fp*fp*(spc(1)+fp*(spc(2)+fp*spc(3)))
        amp=sac(1)+fp*(sac(2)+fp*(sac(3)+fp*(sac(4)+fp*sac(5))))
        corrs(i)=amp*cmplx(cos(ps),sin(ps))
        pd=fp*fp*(dpc(1)+fp*(dpc(2)+fp*dpc(3)))
        corrd(i)=cmplx(cos(pd),sin(pd))
     enddo
  endif

  fac=2.0/nfft
  c(1:npts)=fac*d(1:npts)
  c(npts+1:nfft)=0.
  call four2a(c,nfft,1,-1,1)               !Forward c2c FFT

  if( (.not. bseq) .and. (.not. bdeq) ) then
    c(1:nh+1)=h(1:nh+1)*c(1:nh+1)
  else if( bseq .and. (.not. bdeq) ) then
    c(1:nh+1)=h(1:nh+1)*corrs(1:nh+1)*c(1:nh+1)
  else if( (.not. bseq) .and. bdeq ) then
    c(1:nh+1)=h(1:nh+1)*corrd(1:nh+1)*c(1:nh+1)
  else if( bseq .and. bdeq ) then
    c(1:nh+1)=h(1:nh+1)*corrs(1:nh+1)*corrd(1:nh+1)*c(1:nh+1)
  endif

  c(1)=0.5*c(1)                            !Half of DC term
  c(nh+2:nfft)=0.                          !Zero the negative frequencies
  call four2a(c,nfft,1,1,1)                !Inverse c2c FFT
  return
end subroutine analytic