Changes needed to accommodate channel equalization implemented in the new msk144signalquality.f90. Warning - msk144signalquality writes to the executable directory.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7412 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Steven Franke 2016-12-20 03:01:43 +00:00
parent f9e6c89847
commit e4c5811ac3
4 changed files with 53 additions and 37 deletions

View File

@ -1,29 +1,44 @@
subroutine analytic(d,npts,nfft,c)
subroutine analytic(d,npts,nfft,c,a,equalize)
! Convert real data to analytic signal
parameter (NFFTMAX=1024*1024)
real d(npts)
real h(NFFTMAX/2)
complex h(NFFTMAX/2)
real*8 a(5),alast(5),fp
real ac(5),a0(5)
complex c(NFFTMAX)
logical*1 equalize
data nfft0/0/
save nfft0,h
data alast/0.0,0.0,0.0,0.0,0.0/
data ac/1.0,0.05532,0.11438,0.12918,0.09274/ ! amp coeffs for TS2000
data a0/0.0,0.0,-0.952,0.768,-0.565/ ! baseline phase coeffs for TS2000
save nfft0,h,alast,a0,ac
! disable baseline phase correction for commit - should look for a file with coeffs
a0(1:5)=0.0
df=12000.0/nfft
nh=nfft/2
if(nfft.ne.nfft0) then
if( nfft.ne.nfft0 .or. any(alast .ne. a) ) then
t=1.0/2000.0
beta=0.1
pi=4.0*atan(1.0)
do i=1,nh+1
ff=(i-1)*df
f=ff-1500.0
h(i)=0.
if(abs(f).le.(1-beta)/(2*t)) h(i)=1.0
if(abs(f).gt.(1-beta)/(2*t) .and. abs(f).le.(1+beta)/(2*t)) then
h(i)=0.5*(1+cos((pi*t/beta )*(abs(f)-(1-beta)/(2*t))))
fp=f/1000.0
h(i)=cmplx(1.0,0.0)
if( equalize ) then
phase0=a0(1)+fp*(a0(2)+fp*(a0(3)+fp*(a0(4)+fp*a0(5))))
phase=a(1)+fp*(a(2)+fp*(a(3)+fp*(a(4)+fp*a(5))))
! amp=ac(1)+fp*(ac(2)+fp*(ac(3)+fp*(ac(4)+fp*ac(5))))
amp=1.0 ! no amplitude correction for now
h(i)=amp*cmplx(cos(phase),sin(phase))*cmplx(cos(phase0),sin(phase0))
endif
if(abs(f).gt.(1-beta)/(2*t) .and. abs(f).le.(1+beta)/(2*t)) then
h(i)=h(i)*0.5*(1+cos((pi*t/beta )*(abs(f)-(1-beta)/(2*t))))
endif
! h(i)=sqrt(h(i))
enddo
nfft0=nfft
endif
@ -33,22 +48,10 @@ subroutine analytic(d,npts,nfft,c)
c(npts+1:nfft)=0.
call four2a(c,nfft,1,-1,1) !Forward c2c FFT
! do i=1,nh
! f=(i-1)*df
! s(i)=real(c(i))**2 + aimag(c(i))**2
! write(12,3001) f,s(i),db(s(i))
!3001 format(3f12.3)
! enddo
! ia=700.0/df
! c(1:ia)=0.
! ib=2300.0/df
! c(ib:nfft)=0.
c(1:nh+1)=h(1:nh+1)*c(1:nh+1)
c(1)=0.5*c(1) !Half of DC term
c(nh+2:nfft)=0. !Zero the negative frequencies
call four2a(c,nfft,1,1,1) !Inverse c2c FFT
alast=a
return
end subroutine analytic

View File

@ -1,4 +1,4 @@
subroutine msk144decodeframe(c,msgreceived,nsuccess)
subroutine msk144decodeframe(c,softbits,msgreceived,nsuccess)
! use timer_module, only: timer
parameter (NSPM=864)

View File

@ -1,4 +1,4 @@
subroutine msk144spd(cbig,n,ntol,nsuccess,msgreceived,fc,fret,tret,navg)
subroutine msk144spd(cbig,n,ntol,nsuccess,msgreceived,fc,fret,tret,navg,ct,softbits)
! MSK144 short-ping-decoder
use timer_module, only: timer
@ -23,6 +23,7 @@ subroutine msk144spd(cbig,n,ntol,nsuccess,msgreceived,fc,fret,tret,navg)
real rcw(12)
real ferrs(MAXCAND)
real snrs(MAXCAND)
real softbits(144)
real tonespec(NFFT)
real tpat(NPATTERNS)
real*8 dt, df, fs, pi, twopi
@ -177,7 +178,7 @@ subroutine msk144spd(cbig,n,ntol,nsuccess,msgreceived,fc,fret,tret,navg)
if( is.eq.2) ic0=max(1,ic0-1)
if( is.eq.3) ic0=min(NSPM,ic0+1)
ct=cshift(c,ic0-1)
call msk144decodeframe(ct,msgreceived,ndecodesuccess)
call msk144decodeframe(ct,softbits,msgreceived,ndecodesuccess)
if( ndecodesuccess .gt. 0 ) then
tret=(nstart(icand)+NSPM/2)/fs

View File

@ -28,10 +28,14 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall, &
real d(NFFT1)
real pow(8)
real softbits(144)
real xmc(NPATTERNS)
real*8 pcoeffs(5)
logical*1 bshmsg,bcontest
logical first
logical*1 equalized
data first/.true./
data iavpatterns/ &
1,1,1,1,0,0,0,0, &
@ -39,18 +43,20 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall, &
1,1,1,1,1,0,0,0, &
1,1,1,1,1,1,1,0/
data xmc/2.0,4.5,2.5,3.5/ !Used to set time at center of averaging mask
save first,tsec0,nutc00,pnoise,nsnrlast,msglast,cdat
save first,tsec0,nutc00,pnoise,nsnrlast,msglast,cdat,pcoeffs,equalized
if(first) then
tsec0=tsec
nutc00=nutc0
pnoise=-1.0
pcoeffs(1:5)=0.0
equalized=.false.
first=.false.
endif
fc=nrxfreq
!!! Dupe checking should probaby be moved to mainwindow.cpp
! Dupe checking setup
if(nutc00.ne.nutc0 .or. tsec.lt.tsec0) then ! reset dupe checker
msglast=' '
nsnrlast=-99
@ -68,7 +74,7 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall, &
fac=1.0/rms
d(1:NZ)=fac*d(1:NZ)
d(NZ+1:NFFT1)=0.
call analytic(d,NZ,NFFT1,cdat) !Convert to analytic signal and filter
call analytic(d,NZ,NFFT1,cdat,pcoeffs,.true.) !Convert to analytic signal and filter
! Calculate average power for each frame and for the entire block.
! If decode is successful, largest power will be taken as signal+noise.
@ -86,10 +92,8 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall, &
! center a 3-frame analysis window and attempts to decode each of the
! 3 frames along with 2- and 3-frame averages.
np=8*NSPM
call msk144spd(cdat,np,ntol,nsuccess,msgreceived,fc,fest,tdec,navg)
call msk144spd(cdat,np,ntol,nsuccess,msgreceived,fc,fest,tdec,navg,ct,softbits)
!############################################################
!##### hardwired for testing - need to bring in Sh box status
if(nsuccess.eq.0 .and. bshmsg) then
call msk40spd(cdat,np,ntol,mycall(1:6),hiscall(1:6),nsuccess, &
msgreceived,fc,fest,tdec,navg)
@ -98,6 +102,7 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall, &
if( nsuccess .eq. 1 ) then
tdec=tsec+tdec
decsym=' & '
if( equalized ) decsym=' ^ '
ipk=0
is=0
goto 900
@ -120,15 +125,17 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall, &
if( nsyncsuccess .eq. 0 ) cycle
do ipk=1,npeaks
do is=1,3
do is=1,3 ! With equalization, this loop may not be necessary
ic0=npkloc(ipk)
if(is.eq.2) ic0=max(1,ic0-1)
if(is.eq.3) ic0=min(NSPM,ic0+1)
ct=cshift(c,ic0-1)
call msk144decodeframe(ct,msgreceived,ndecodesuccess)
call msk144decodeframe(ct,softbits,msgreceived,ndecodesuccess)
if(ndecodesuccess .gt. 0) then
tdec=tsec+xmc(iavg)*tframe
decsym=' & '
if( equalized ) decsym=' ^ '
if( equalized .and. is .ne. 1 ) decsym=' ! ' !help decide if is dither is needed
goto 900
endif
enddo !Slicer dither
@ -149,7 +156,7 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall, &
go to 999
900 continue
! Successful decode - estimate snr !!! noise estimate needs work
! Successful decode - estimate snr
if( pnoise .gt. 0.0 ) then
snr0=10.0*log10(pmax/pnoise-1.0)
else
@ -157,7 +164,10 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall, &
endif
nsnr=nint(snr0)
!!!! Temporary - dupe check. Only print if new message, or higher snr.
call msk144signalquality(ct,snr0,fest,tdec,softbits,msgreceived,hiscall, &
ncorrected,eyeopening,equalized,pcoeffs)
! Dupe check. Only print if new message, or higher snr.
if(msgreceived.ne.msglast .or. nsnr.gt.nsnrlast .or. tsec.lt.tsec0) then
msglast=msgreceived
nsnrlast=nsnr
@ -168,8 +178,8 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall, &
call fix_contest_msg(mycall(1:6),mygrid,hiscall(1:6),msgreceived)
endif
write(line,1020) nutc0,nsnr,tdec,nint(fest),decsym,msgreceived, &
navg,char(0)
1020 format(i6.6,i4,f5.1,i5,a3,a22,i2,a1)
navg,ncorrected,eyeopening,char(0)
1020 format(i6.6,i4,f5.1,i5,a3,a22,i2,i3,f5.1,a1)
endif
999 tsec0=tsec
@ -177,3 +187,5 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall, &
end subroutine mskrtd
include 'fix_contest_msg.f90'
include 'msk144signalquality.f90'