WSJT-X/libm65/specjtms.f90

subroutine specjtms(k,px,pxsmo,spk0,f0)

! Starting code for a JTMS3 decoder.

  parameter (NSMAX=30*48000)
  parameter (MAXFFT=8192)
  integer*2 id
  real x(MAXFFT)
  complex cx(MAXFFT),cx2(MAXFFT)
  logical first
  common/mscom/id(1440000),s1(215),s2(215)
  data first/.true./
  save
 
  if(first) then
     pi=4.0*atan(1.0)
     twopi=2.0*pi
     kstep=2048
     first=.false.
     sqsmo=0.
  endif

  t=k/48000.0
  nfft=4096
  df=48000.0/nfft
  nh=nfft/2
  ib=k
  ia=k-kstep+1
  i0=k-nfft+1
  sq=0.
  do i=ia,ib
     d=id(i)
     sq=sq + d*d
  enddo
  sq=sq/2048.0
  sqsmo=0.95*sqsmo + 0.05*sq
  rms=sqrt(sq)
  px=db(sq) - 23.0
  pxsmo=db(sqsmo) - 23.0
!  write(13,1010) t,rms,sq,px,pxsmo
!1010 format(5f12.3)
  if(k.lt.nfft) return


  x(1:nfft)=id(i0:ib)
  fac=0.002/nfft
  cx=fac*x
  call four2a(cx,nfft,1,-1,1)                    !Forward c2c FFT

  iz=nint(2500.0/df)

  do i=1,iz                                      !Save spectrum for plot
     s1(i)=real(cx(i+1))**2 + aimag(cx(i+1))**2
  enddo

  cx(1)=0.5*cx(1)
  cx(nh+2:nfft)=0.
  call four2a(cx,nfft,1,1,1)                     !Inverse c2c FFT

  cx2(1:nfft)=cx(1:nfft)*cx(1:nfft)
  cx2(nfft+1:)=0.0

  nfft=8192
  df=48000.0/nfft
  
  call four2a(cx2,nfft,1,-1,1)                   !Forward c2c FFT of cx2

  spk0=0.
  j0=nint(2.0*1428.57/df)
  ja=j0-107
  jb=j0+107
  do j=ja,jb
     sq=1.e-4*(real(cx2(j))**2 + aimag(cx2(j))**2)
     s2(j-ja+1)=sq
     if(sq.gt.spk0) then
        spk0=sq
        f=(j-1)*df
        f0=0.5*(f-3000.0)
        phi0=0.5*atan2(aimag(cx2(j)),real(cx2(j)))
     endif
  enddo
  spk0=0.5*db(spk0) - 7.0
!  write(14,3001) k/2048,spk0,f0,phi0
!3001 format(i8,3f12.3)
  call flush(14)

  return
end subroutine specjtms