First working QRA66 decoder.

This commit is contained in:
Joe Taylor 2020-08-01 09:24:59 -04:00
parent c4ef1e3e25
commit f45c617422
5 changed files with 269 additions and 3 deletions

View File

@ -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

View File

@ -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

157
lib/qra66_decode.f90 Normal file
View File

@ -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

44
lib/spec66.f90 Normal file
View File

@ -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

60
lib/sync66.f90 Normal file
View File

@ -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