Starting to implement a decoder for the wide/slow JT9 submodes.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@6575 ab8295b8-cf94-4d9e-aec4-7959e3be5d79

lib/Makefile.jt9w

@@ -0,0 +1,46 @@
+
+# Set paths
+EXE_DIR = ..\\..\\wsjtx_install
+QT_DIR = C:/wsjt-env/Qt5/5.2.1/mingw48_32
+FFTW3_DIR = ..
+
+INCPATH = -I${QT_DIR}/include/QtCore -I${QT_DIR}/include 
+
+# Compilers
+CC = gcc
+CXX = g++
+FC = gfortran
+AR = ar cr
+RANLIB = ranlib
+MKDIR = mkdir -p
+CP = cp
+RM = rm -f
+
+FFLAGS = -O2 -fbounds-check -Wall -Wno-conversion
+CFLAGS = -O2 -I. 
+
+# Default rules
+%.o: %.c
+	${CC} ${CFLAGS} -c $<
+%.o: %.f
+	${FC} ${FFLAGS} -c $<
+%.o: %.F
+	${FC} ${FFLAGS} -c $<
+%.o: %.f90
+	${FC} ${FFLAGS} -c $<
+%.o: %.F90
+	${FC} ${FFLAGS} -c $<
+
+all:	 jt9w
+
+OBJS1 = jt9w.o smo.o sync9w.o pctile.o shell.o lorentzian.o fchisq0.o \
+      softsym9w.o four2a.o interleave9.o jt9fano.o fano232.o packjt.o \
+      deg2grid.o grid2deg.o fmtmsg.o db.o
+      
+jt9w: $(OBJS1) 
+	$(FC) -o jt9w $(OBJS1) -lfftw3f
+
+.PHONY : clean
+
+clean:
+	$(RM) *.o JTMSKcode JTMSKcode.exe

lib/jt9w.f90

@@ -0,0 +1,72 @@
+program jt9w
+
+  parameter (NSMAX=6827,NZMAX=60*12000)
+  real ss(184,NSMAX)
+  real ref(NSMAX)
+  real ccfred(NSMAX)
+  real ccfblue(-9:18)
+  real a(5)
+  integer*2 id2(NZMAX)
+  character*12 arg
+  character*22 decoded
+  integer*1 i1SoftSymbols(207)
+
+  call getarg(1,arg)
+  read(arg,*) iutc
+
+  open(20,file='refspec.dat',status='old')
+  do i=1,NSMAX
+     read(20,*) j,freq,ref(i)
+  enddo
+
+  df=12000.0/16384.0
+  nsps=6912
+  tstep=nsps*0.5/12000.0
+  npts=52*12000
+  limit=10000
+  ntol=100
+
+  do ifile=1,999
+     read(60,end=999) nutc,nfqso,ntol,ndepth,nmode,nsubmode,nzhsym,ss,id2
+     if(nutc.ne.iutc) cycle
+
+     ia=nint((nfqso-ntol)/df)
+     ib=nint((nfqso+ntol)/df)
+     lag1=-int(2.5/tstep + 0.9999)
+     lag2=int(5.0/tstep + 0.9999)
+     nhsym=184
+     do iter=1,2
+        nadd=3
+        if(iter.eq.2) nadd=2*nint(0.375*a(4)) + 1
+        call sync9w(ss,nhsym,lag1,lag2,ia,ib,ccfred,ccfblue,ipk,lagpk,nadd)
+        sum1=sum(ccfblue) - ccfblue(lagpk-1)-ccfblue(lagpk) -ccfblue(lagpk+1)
+        sq=dot_product(ccfblue,ccfblue) - ccfblue(lagpk-1)**2 - &
+             ccfblue(lagpk)**2 - ccfblue(lagpk+1)**2
+        base=sum1/25.0
+        rms=sqrt(sq/24.0)
+        snr=(ccfblue(lagpk)-base)/rms
+        nsnr=db(snr)-29.7
+        xdt0=lagpk*tstep
+        
+        call lorentzian(ccfred(ia),ib-ia+1,a)
+        f0=(ia+a(3))*df
+!        write(*,3001) nadd,a,base,rms,snr,xdt,f0,ipk,lagpk
+!3001    format(i3,9f7.2,f7.1,2i5)
+
+        ccfblue=(ccfblue-base)/rms
+!        rewind 16
+!        do lag=lag1,lag2
+!           write(16,3002) lag*tstep,ccfblue(lag)
+!3002       format(2f10.3)
+!        enddo
+
+     enddo
+
+     call softsym9w(id2,npts,xdt0,f0,a(4)*df,nsubmode,xdt1,i1softsymbols)
+     call jt9fano(i1softsymbols,limit,nlim,decoded)
+     write(*,1100) nutc,nsnr,xdt1-1.0,nint(f0),decoded,nlim
+1100 format(i4.4,i4,f5.1,i5,1x,'@',1x,a22,i10)
+
+  enddo
+
+999 end program jt9w

lib/lorentzian.f90

@@ -1,11 +1,12 @@
 subroutine lorentzian(y,npts,a)
 
 ! Input:  y(npts); assume x(i)=i, i=1,npts
-! Output: a(4)
+! Output: a(5)
 !         a(1) = baseline
 !         a(2) = amplitude
 !         a(3) = x0
 !         a(4) = width
+!         a(5) = chisqr
 
   real y(npts)
   real a(5)

lib/softsym9w.f90

@@ -0,0 +1,127 @@
+subroutine softsym9w(id2,npts,xdt0,f0,width,nsubmode,xdt1,i1softsymbols)
+
+  parameter (NFFT=6912,NH=NFFT/2,NQ=NH/2)
+  real s(NQ)
+  real s2(0:8,85)
+  real s3(0:7,69)
+  real x(NFFT)
+  complex cx(0:NH)
+  integer*2 id2(60*12000)
+  integer*1 i1SoftSymbolsScrambled(207)
+  integer*1 i1softsymbols(207)
+  include 'jt9sync.f90'
+  equivalence (x,cx)
+
+  df=12000.0/NFFT
+!  i0a=1
+!  i0b=5*12000
+  i0a=max(0.0,(xdt0-1.0)*12000.0)
+  i0b=(xdt0+1.0)*12000.0
+  k1=nint((f0-0.5*width)/df)
+  k2=nint((f0+0.5*width)/df)
+  smax=0.
+  i0pk=0
+  i1softsymbols=0
+
+  do i0=i0a,i0b,432
+     s=0.
+     ssum=0.
+     do j=1,16
+        ia=i0 + (ii(j)-1)*nfft
+        ib=ia+NFFT-1
+        x=1.e-6*id2(ia:ib)
+        call four2a(x,nfft,1,-1,0)        !r2c FFT
+        do k=1,NQ
+           s(k)=s(k) + real(cx(k))**2 + aimag(cx(k))**2
+        enddo
+     enddo
+     ssum=ssum + sum(s(k1:k2))
+     if(ssum.gt.smax) then
+        smax=ssum
+        i0pk=i0
+     else 
+        if(ssum.lt.0.7*smax) exit
+     endif
+  end do
+  xdt1=(i0pk-1)/12000.0
+
+  if(i0pk.le.0) go to 999
+
+  m=0
+  do j=1,85
+     ia=i0pk + (j-1)*nfft
+     ib=ia+NFFT-1
+     x=1.e-6*id2(ia:ib)
+     call four2a(x,nfft,1,-1,0)        !r2c FFT
+     do k=1,NQ
+        s(k)=real(cx(k))**2 + aimag(cx(k))**2
+     enddo
+
+     dtone=df*(2**nsubmode)
+     do i=0,8
+        f=f0 + i*dtone
+        k1=nint((f-0.5*width)/df)
+        k2=nint((f+0.5*width)/df)
+        s2(i,j)=sum(s(k1:k2))                  !Symbol spectra, including sync
+     enddo
+
+     if(isync(j).eq.0) then
+        m=m+1
+        s3(0:7,m)=s2(1:8,j)                   !Symbol spectra, data only
+     endif
+
+!     write(19,3101) j,s2(0:8,j)
+!3101 format(i2,9f8.2)
+  enddo
+
+  ss=0.
+  sig=0.
+  do j=1,69
+     smax=0.
+     do i=0,7
+        smax=max(smax,s3(i,j))
+        ss=ss+s3(i,j)
+     enddo
+     sig=sig+smax
+     ss=ss-smax
+  enddo
+  ave=ss/(69*7)                           !Baseline
+  call pctile(s2,9*85,35,xmed)
+  s3=s3/ave
+  sig=sig/69.                             !Signal
+  t=max(1.0,sig - 1.0)
+  snrdb=db(t)
+     
+  m0=3
+  k=0
+  do j=1,69
+        smax=0.
+        do i=0,7
+           if(s3(i,j).gt.smax) smax=s3(i,j)
+        enddo
+
+     do m=m0-1,0,-1                   !Get bit-wise soft symbols
+        if(m.eq.2) then
+           r1=max(s3(4,j),s3(5,j),s3(6,j),s3(7,j))
+           r0=max(s3(0,j),s3(1,j),s3(2,j),s3(3,j))
+        else if(m.eq.1) then
+           r1=max(s3(2,j),s3(3,j),s3(4,j),s3(5,j))
+           r0=max(s3(0,j),s3(1,j),s3(6,j),s3(7,j))
+        else
+           r1=max(s3(1,j),s3(2,j),s3(4,j),s3(7,j))
+           r0=max(s3(0,j),s3(3,j),s3(5,j),s3(6,j))
+        endif
+
+        k=k+1
+        i4=nint(10.0*(r1-r0))
+        if(i4.lt.-127) i4=-127
+        if(i4.gt.127) i4=127
+        i1SoftSymbolsScrambled(k)=i4
+     enddo
+  enddo
+
+! Remove interleaving
+  call interleave9(i1SoftSymbolsScrambled,-1,i1SoftSymbols)
+
+999  return
+end subroutine softsym9w

lib/sync9w.f90

@@ -0,0 +1,79 @@
+subroutine sync9w(ss,nzhsym,lag1,lag2,ia,ib,ccfred,ccfblue,ipkbest,lagpk,nadd)
+
+  include 'constants.f90'
+  real ss(184,NSMAX)
+  real ss1(184),ss1save(184)
+  real ccfred(NSMAX)
+  real ccfblue(-9:18)
+  real sa(NSMAX),sb(NSMAX)
+  include 'jt9sync.f90'
+
+! Smooth the symbol spectra (by an amount consistent with measured width??)
+  do j=1,nzhsym
+     sa=ss(j,1:NSMAX)
+     call smo(sa,NSMAX,sb,nadd)
+     call smo(sb,NSMAX,sa,nadd)
+     ss(j,1:NSMAX)=sa
+  enddo
+
+  ipk=0
+  ipkbest=0
+  sbest=0.
+  ccfred=0.
+  df=12000.0/16384.0
+
+  do i=ia,ib                                 !Loop over specified freq range
+     ss1=ss(1:184,i)                         !Symbol amplitudes at this freq
+     call pctile(ss1,nzhsym,50,xmed)         !Median level at this freq
+     ss1=ss1/xmed - 1.0
+
+     smax=0.                                 !Find DT in specified range
+     do lag=lag1,lag2
+        sum1=0.
+        nsum=nzhsym
+        do j=1,16                            !Sum over 16 sync symbols
+           k=ii2(j) + lag
+           if(k.ge.1 .and. k.le.nzhsym) then
+              sum1=sum1 + ss1(k)
+              nsum=nsum-1
+           endif
+        enddo
+        if(sum1.gt.smax) then
+           smax=sum1
+           ipk=i 
+        endif
+     enddo
+
+     ccfred(i)=smax                        !Best at this freq, over all lags
+     if(smax.gt.sbest) then
+        sbest=smax
+        ipkbest=ipk
+        ss1save=ss1
+     endif
+  enddo
+
+  call pctile(ccfred(ia),ib-ia+1,50,xmed)
+  if(xmed.le.0.0) xmed=1.0
+  ccfred=ccfred/xmed
+
+  ss1=ss1save
+  smax=0.                                 !Find DT in specified range
+  do lag=lag1,lag2
+     sum1=0.
+     nsum=nzhsym
+     do j=1,16                            !Sum over 16 sync symbols
+        k=ii2(j) + lag
+        if(k.ge.1 .and. k.le.nzhsym) then
+           sum1=sum1 + ss1(k)
+           nsum=nsum-1
+        endif
+     enddo
+     ccfblue(lag)=sum1
+     if(sum1.gt.smax) then
+        smax=sum1
+        lagpk=lag
+     endif
+  enddo
+
+  return
+end subroutine sync9w