diff --git a/lib/Makefile.jt9w b/lib/Makefile.jt9w new file mode 100644 index 000000000..3d44edb7c --- /dev/null +++ b/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 diff --git a/lib/jt9w.f90 b/lib/jt9w.f90 new file mode 100644 index 000000000..b3fb62281 --- /dev/null +++ b/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 diff --git a/lib/lorentzian.f90 b/lib/lorentzian.f90 index 6b88a4d94..33d6ff661 100644 --- a/lib/lorentzian.f90 +++ b/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) diff --git a/lib/softsym9w.f90 b/lib/softsym9w.f90 new file mode 100644 index 000000000..b3b08eb73 --- /dev/null +++ b/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 diff --git a/lib/sync9w.f90 b/lib/sync9w.f90 new file mode 100644 index 000000000..f6797535c --- /dev/null +++ b/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