WSJT-X/libm65/specjtms.f90

subroutine specjtms(k)

! Starting code for a JTMS3 decoder.

  parameter (NSMAX=30*48000)
  parameter (NFFT=8192,NH=NFFT/2)
  integer*2 id
  real x(NFFT),w(NFFT)
  real p(24)
  real chansym(258),softsym(341)
  integer nsum(24)
  complex cx(NFFT),cx2(NFFT),cx0(NFFT)
  complex covx(NH)
  real s1a(NH),s2a(580)
  logical first,window
  common/mscom/id(1440000),s1(215,703),s2(215,703)
  data first/.true./
  save

  if(first) then
     pi=4.0*atan(1.0)
     twopi=2.0*pi
     do i=1,nfft
        w(i)=(sin(i*pi/nfft))**2
     enddo
     df=48000.0/nfft
     ja=nint(2600.0)/df
     jb=nint(3400.0)/df
     iz=3000.0/df
     covx=0.
     read(10,3001) chansym
3001 format(50f1.0)
     chansym=2.0*chansym - 1.0
     window=.false.
     first=.false.
  endif

  ib=k
  ia=k-4095
  i0=ib-8191
  sq=0.
  do i=ia,ib
     sq=sq + (0.001*id(i))**2
  enddo
  write(13,1010) t,sq,db(sq)
1010 format(3f12.3)
  if(k.lt.8192) return

  x(1:nfft)=0.001*id(i0:ib)

  fac=2.0/nfft
  cx=fac*x
  if(window) cx=w*cx
  call four2a(cx,nfft,1,-1,1)                    !Forward c2c FFT

  do i=1,iz                                      !Save spectrum for plot
     s1a(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
  if(window) then
     cx(1:nh)=cx(1:nh)+covx(1:nh)               !Add previous segment's 2nd half
     covx(1:nh)=cx(nh+1:nfft)                    !Save 2nd half
  endif

  t=k/48000.0
  cx2=cx*cx
  call four2a(cx2,nfft,1,-1,1)                   !Forward c2c FFT of cx2

  spk0=0.
  do j=ja,jb
     sq=1.e-6*(real(cx2(j))**2 + aimag(cx2(j))**2)
     s2a(j)=sq
     f=(j-1)*df
     if(sq.gt.spk0) then
        spk0=sq
        f0=0.5*(f-3000.0)
        phi0=0.5*atan2(aimag(cx2(j)),real(cx2(j)))
     endif
     write(15,1020) (j-1)*df,sq
1020 format(f10.3,f12.3)
  enddo

  slimit=2.0
!  slimit=87.5
!  if(spk0.gt.slimit) then
  if(abs(spk0-87.3).lt.0.1) then
     write(*,1030) t,f0,phi0,spk0
1030 format('t:',f6.2,'   f0:',f7.1,'   phi0:',f6.2,'   spk0:',f8.1)
     do i=1,iz
        write(16,1040) i*df,s1a(i),db(s1a(i))
1040    format(3f12.3)
     enddo
     do j=ja,jb
        f=(j-1)*df
        f0a=0.5*(f-3000.0)
        write(17,1050) f0a,s2a(j)
1050    format(2f12.3)
     enddo

     phi=phi0
     phi=3.9
     dphi=2.0*pi*(f0+1500.0 -1.1)/48000.0
     p=0.
     nsum=0
     do i=1,nfft
        phi=phi+dphi
        if(phi.gt.twopi) phi=phi-twopi        
        cx0(i)=cx(i)*cmplx(cos(phi),-sin(phi))
        pha=atan2(aimag(cx0(i)),real(cx0(i)))
        write(18,1060) i,cx0(i),pha
1060    format(i6,5f12.3)
        j=mod(i-1,24) + 1
!        p(j)=p(j)+abs(cx0(i))
        p(j)=p(j) + real(cx0(i))**2 + aimag(cx0(i))**2
        nsum(j)=nsum(j)+1
     enddo

     do i=1,24
        p(i)=p(i)/nsum(i)
        write(20,1070) i,p(i)
1070    format(i6,f12.3)
     enddo

     do i=16,nfft,24
        amp=abs(cx0(i))
        pha=atan2(aimag(cx0(i)),real(cx0(i)))
        j=(i+23)/24
        write(21,1060) j,cx0(i),pha,pha+twopi,amp
        softsym(j)=real(cx0(i))
     enddo

!     do iter=1,5
     chansym=cshift(chansym,-86)
     do lag=0,83
        sum=dot_product(chansym,softsym(lag+1:lag+258))
        if(abs(sum).gt.smax) then
           smax=abs(sum)
           lagpk=lag
        endif
        write(22,1080) lag,sum
1080    format(i3,f12.3)
     enddo
!     chansym=cshift(chansym,43)
!     enddo

     do i=1,258
        prod=-chansym(i)*softsym(lagpk+i)
        write(23,1090) i,prod,chansym(i),softsym(lagpk+i)
1090    format(i5,3f10.3)
     enddo

     do i=1,258,6
        write(24,1100) (i+5)/6,int(chansym(i)),softsym(lagpk+i)
1100    format(2i5,f8.1)
     enddo
  endif

end subroutine specjtms