From f45c6174226032596d7e4e3a3893bae30f64c3e8 Mon Sep 17 00:00:00 2001
From: Joe Taylor <joe@princeton.edu>
Date: Sat, 1 Aug 2020 09:24:59 -0400
Subject: [PATCH] First working QRA66 decoder.
---
CMakeLists.txt | 1 +
lib/qra/qra66/qra66sim.f90 | 10 ++-
lib/qra66_decode.f90 | 157 +++++++++++++++++++++++++++++++++++++
lib/spec66.f90 | 44 +++++++++++
lib/sync66.f90 | 60 ++++++++++++++
5 files changed, 269 insertions(+), 3 deletions(-)
create mode 100644 lib/qra66_decode.f90
create mode 100644 lib/spec66.f90
create mode 100644 lib/sync66.f90
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9726caee0..49cb6ad57 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -565,6 +565,7 @@ set (wsjt_FSRCS
lib/softsym9w.f90
lib/shell.f90
lib/spec64.f90
+ lib/spec66.f90
lib/spec9f.f90
lib/stdmsg.f90
lib/subtract65.f90
diff --git a/lib/qra/qra66/qra66sim.f90 b/lib/qra/qra66/qra66sim.f90
index 934239189..a482e4cc3 100644
--- a/lib/qra/qra66/qra66sim.f90
+++ b/lib/qra/qra66/qra66sim.f90
@@ -47,12 +47,17 @@ program qra66sim
nsps=960
nsym=169
endif
+
+ ichk=66 !Flag sent to genqra64
+ call genqra64(msg,ichk,msgsent,itone,itype)
+ write(*,1001) itone
+1001 format('Channel symbols:'/(20i3))
+
baud=12000.d0/nsps !Keying rate = 6.25 baud
h=default_header(12000,npts)
- ichk=66 !Flag sent to genqra64
write(*,1000)
-1000 format('File Freq A|B S/N DT Dop Message'/60('-'))
+1000 format('File Freq A|B S/N DT Dop Message'/60('-'))
do ifile=1,nfiles !Loop over requested number of files
write(fname,1002) ifile !Output filename
@@ -70,7 +75,6 @@ program qra66sim
bandwidth_ratio=2500.0/6000.0
sig=sqrt(2*bandwidth_ratio)*10.0**(0.05*snrdb)
if(snrdb.gt.90.0) sig=1.0
- call genqra64(msg,ichk,msgsent,itone,itype)
write(*,1020) ifile,f0,csubmode,xsnr,xdt,fspread,msgsent
1020 format(i4,f10.3,2x,a1,2x,f5.1,f6.2,f6.1,1x,a22)
phi=0.d0
diff --git a/lib/qra66_decode.f90 b/lib/qra66_decode.f90
new file mode 100644
index 000000000..53243e78c
--- /dev/null
+++ b/lib/qra66_decode.f90
@@ -0,0 +1,157 @@
+module qra66_decode
+
+ type :: qra66_decoder
+ procedure(qra66_decode_callback), pointer :: callback
+ contains
+ procedure :: decode
+ end type qra66_decoder
+
+ abstract interface
+ subroutine qra66_decode_callback (this,nutc,sync,nsnr,dt,freq, &
+ decoded,nap,qual,ntrperiod,fmid,w50)
+ import qra66_decoder
+ implicit none
+ class(qra66_decoder), intent(inout) :: this
+ integer, intent(in) :: nutc
+ real, intent(in) :: sync
+ integer, intent(in) :: nsnr
+ real, intent(in) :: dt
+ real, intent(in) :: freq
+ character(len=37), intent(in) :: decoded
+ integer, intent(in) :: nap
+ real, intent(in) :: qual
+ integer, intent(in) :: ntrperiod
+ real, intent(in) :: fmid
+ real, intent(in) :: w50
+ end subroutine qra66_decode_callback
+ end interface
+
+contains
+
+ subroutine decode(this,callback,iwave,nutc,nfa,nfb,nfqso,ndepth,lapdx, &
+ mycall,hiscall,hisgrid)
+
+ use timer_module, only: timer
+ use packjt
+ use, intrinsic :: iso_c_binding
+ parameter (NFFT1=15*12000,NFFT2=15*6000)
+ class(qra66_decoder), intent(inout) :: this
+ procedure(qra66_decode_callback) :: callback
+ character(len=12) :: mycall, hiscall
+ character(len=6) :: hisgrid
+ character*37 decoded
+ integer*2 iwave(NFFT1) !Raw data
+ integer ipk(1)
+ integer dat4(12)
+ logical lapdx,ltext
+ complex c0(0:NFFT1-1) !Spectrum, then analytic signal
+ real s(900)
+ real s3a(-64:127,63)
+ real s3(-64:127,63)
+ real a(5)
+ data nc1z/-1/,nc2z/-1/,ng2z/-1/,maxaptypez/-1/
+ save
+
+ this%callback => callback
+ nsps=1920
+ baud=12000.0/nsps
+ df1=12000.0/NFFT1
+
+ if(nutc.eq.-999) print*,mycall,hiscall,hisgrid,lapdx,ndepth,nfa,nfb,nfqso
+
+! Prime the QRA decoder for possible use of AP
+ call packcall(mycall(1:6),nc1,ltext)
+ call packcall(hiscall(1:6),nc2,ltext)
+ call packgrid(hisgrid(1:4),ng2,ltext)
+ nSubmode=0
+ b90=1.0
+ nFadingModel=1
+ maxaptype=4
+ if(iand(ndepth,64).ne.0) maxaptype=5
+ if(nc1.ne.nc1z .or. nc2.ne.nc2z .or. ng2.ne.ng2z .or. &
+ maxaptype.ne.maxaptypez) then
+ do naptype=0,maxaptype
+ if(naptype.eq.2 .and. maxaptype.eq.4) cycle
+ call qra64_dec(s3,nc1,nc2,ng2,naptype,1,nSubmode,b90, &
+ nFadingModel,dat4,snr2,irc)
+ enddo
+ nc1z=nc1
+ nc2z=nc2
+ ng2z=ng2
+ maxaptypez=maxaptype
+ endif
+ naptype=maxaptype
+
+! Compute the full-length spectrum
+ fac=2.0/NFFT1
+ c0=fac*iwave
+ call four2a(c0,NFFT1,1,-1,1) !Forward c2c FFT
+
+ nadd=101
+ nh=nadd/2
+ df2=nh*df1
+ iz=3000.0/df2
+ do i=1,iz !Compute smoothed spectrum
+ s(i)=0.
+ j0=nh*i
+ do j=j0-nh,j0+nh
+ s(i)=s(i) + real(c0(j))**2 + aimag(c0(j))**2
+ enddo
+ enddo
+ call smo121(s,iz)
+ ipk=maxloc(s)
+ f0=df2*ipk(1) - 0.5*baud !Candidate sync frequency
+
+! do i=1,iz
+! write(51,3051) i*df2,s(i)
+!3051 format(f12.6,e12.3)
+! enddo
+
+ c0(NFFT2/2+1:NFFT2)=0. !Zero the top half
+ c0(0)=0.5*c0(0)
+ call four2a(c0,nfft2,1,1,1) !Inverse c2c FFT
+ call sync66(c0,f0,jpk,sync) !c0 is analytic signal at 6000 S/s
+ xdt=jpk/6000.0 - 0.5
+
+ a=0.
+ a(1)=-(f0 + 1.5*baud)
+ call twkfreq(c0,c0,85*NSPS,6000.0,a)
+ call spec66(c0(jpk:jpk+85*NSPS-1),s3a)
+ s3=s3a/maxval(s3a)
+! do j=1,63
+! ipk=maxloc(s3(-64:127,j))
+! write(54,3054) j,ipk(1)-65
+!3054 format(2i8)
+! do i=-64,127
+! write(53,3053) i,2*s3(i,j)+j-1
+!3053 format(i8,f12.6)
+! enddo
+! enddo
+
+ call qra64_dec(s3a,nc1,nc2,ng2,naptype,0,nSubmode,b90, &
+ nFadingModel,dat4,snr2,irc)
+ if(irc.gt.0) call badmsg(irc,dat4,nc1,nc2,ng2)
+
+ decoded=' '
+ if(irc.ge.0) then
+ call unpackmsg(dat4,decoded) !Unpack the user message
+ call fmtmsg(decoded,iz)
+ if(index(decoded,"000AAA ").ge.1) then
+! Suppress a certain type of garbage decode.
+ decoded=' '
+ irc=-1
+ endif
+ nft=100 + irc
+ nsnr=nint(snr2)
+ else
+ snr2=0.
+ call this%callback(nutc,sYNC,nsnr,xdt,fsig,decoded, &
+ iaptype,qual,ntrperiod,fmid,w50)
+ endif
+! print*,'QRA66 decoder',nutc,jpk,xdt,f0,sync,snr2,irc,decoded
+ print*,decoded
+
+ return
+ end subroutine decode
+
+end module qra66_decode
diff --git a/lib/spec66.f90 b/lib/spec66.f90
new file mode 100644
index 000000000..9242f6a59
--- /dev/null
+++ b/lib/spec66.f90
@@ -0,0 +1,44 @@
+subroutine spec66(c0,s3)
+
+ parameter (LL=3*64) !Frequency channels
+ parameter (NN=63) !Data symbols
+ parameter (NSPS=960) !Samples per symbol at 6000 Hz
+ parameter (NMAX=85*NSPS)
+ complex c0(0:NMAX-1) !Synchrinized complex data
+ complex cs(0:NSPS-1) !Complex symbol spectrum
+ real s3(LL,NN) !Synchronized symbol spectra
+ real xbase0(LL),xbase(LL)
+
+ fac=1.0/NSPS
+ ja=-NSPS
+ do j=1,NN
+ ja=ja+NSPS
+ if(mod(ja/NSPS,4).eq.0) ja=ja+NSPS
+ jb=ja+NSPS-1
+ cs=fac*c0(ja:jb)
+ call four2a(cs,NSPS,1,-1,1) !c2c FFT to frequency
+ do ii=1,LL
+ i=ii-65
+ if(i.lt.0) i=i+NSPS
+ s3(ii,j)=real(cs(i))**2 + aimag(cs(i))**2
+ enddo
+ enddo
+
+ df=6000.0/NSPS
+ do i=1,LL
+ call pctile(s3(i,1:NN),NN,45,xbase0(i)) !Get baseline for passband shape
+ enddo
+
+ nh=9
+ xbase(1:nh-1)=sum(xbase0(1:nh-1))/(nh-1.0)
+ xbase(LL-nh+1:LL)=sum(xbase0(LL-nh+1:LL))/(nh-1.0)
+ do i=nh,LL-nh
+ xbase(i)=sum(xbase0(i-nh+1:i+nh))/(2*nh+1) !Smoothed passband shape
+ enddo
+
+ do i=1,LL
+ s3(i,1:NN)=s3(i,1:NN)/(xbase(i)+0.001) !Apply frequency equalization
+ enddo
+
+ return
+end subroutine spec66
diff --git a/lib/sync66.f90 b/lib/sync66.f90
new file mode 100644
index 000000000..bc8b41805
--- /dev/null
+++ b/lib/sync66.f90
@@ -0,0 +1,60 @@
+subroutine sync66(c0,f0,jpk,sync)
+
+! Use the sync vector to find xdt (and improved f0 ?)
+
+ PARAMETER (NMAX=15*6000)
+ parameter (NSPS=960) !Samples per symbol at 6000 Hz
+ complex c0(0:NMAX-1) !Complex data at 6000 S/s
+ complex csync0(0:NSPS-1)
+ complex csync1(0:NSPS-1)
+ complex z
+ integer b11(11) !Barker 11 code
+ data b11/1,1,1,0,0,0,1,0,0,1,0/ !Barker 11 definition
+ data mode66z/-1/
+ save
+
+ twopi=8.0*atan(1.0)
+ baud=6000.0/NSPS
+ dt=1.0/6000.0
+
+ dphi0=twopi*f0**dt
+ dphi1=twopi*(f0+baud)*dt
+ phi0=0.
+ phi1=0.
+
+ do i=0,NSPS-1 !Compute csync for f0 and f0+baud
+ csync0(i)=cmplx(cos(phi0),sin(phi0))
+ csync1(i)=cmplx(cos(phi1),sin(phi1))
+ phi0=phi0 + dphi0
+ phi1=phi1 + dphi1
+ if(phi0.gt.twopi) phi0=phi0-twopi
+ if(phi1.gt.twopi) phi1=phi1-twopi
+ enddo
+
+ sqmax=0.
+ jstep=NSPS/8
+ do j0=0,6000,jstep
+ sq=0.
+ do k=1,22
+ i=k
+ if(i.gt.11) i=i-11
+ j1=j0 + 4*(k-1)*NSPS
+ if(b11(i).eq.0) then
+ z=dot_product(c0(j1:j1+NSPS-1),csync0)
+ else
+ z=dot_product(c0(j1:j1+NSPS-1),csync1)
+ endif
+ z=0.001*z
+ sq=sq + real(z)**2 + aimag(z)**2
+ enddo
+ write(52,3052) j0/6000.0,sq,j0,j1
+3052 format(f10.6,f12.3,2i8)
+ if(sq.gt.smax) then
+ smax=sq
+ jpk=j0
+ endif
+ enddo
+ sync=smax
+
+ return
+end subroutine sync66