WSJT-X/lib/msk144signalquality.f90
Steven Franke 9582ccd5d5 Lower SWL SNR thresholds by one notch. Note that the xsnr that is compared to the threshold is the nonlinear snr from squared signal spectrum, not 'real' snr.
git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7441 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
2017-01-02 03:56:29 +00:00

229 lines
6.3 KiB
Fortran

subroutine msk144signalquality(cframe,snr,freq,t0,softbits,msg,dxcall, &
nbiterrors,eyeopening,trained,pcoeffs)
character*22 msg,msgsent
character*12 dxcall
character*12 training_dxcall
character*12 trained_dxcall
character*6 mygrid
character*40 pcoeff_filename
character*8 date
character*10 time
character*5 zone
complex cframe(864)
complex cross(864)
complex cross_avg(864)
complex canalytic(1024)
complex cmodel(1024)
integer i4tone(144)
integer hardbits(144)
integer msgbits(144)
integer values(8)
logical*1 bcontest
logical*1 first
logical*1 currently_training
logical*1 trained
logical*1 msg_has_dxcall
logical*1 is_training_frame
real abscross_avg(864)
real softbits(144)
real waveform(0:863)
real d(1024)
real phase(864)
real twopi,freq,phi,dphi0,dphi1,dphi
real*8 x(145),y(145),pp(145),sigmay(145),a(5),chisqr
real pcoeffs(3)
data first/.true./
save cross_avg,abscross_avg,wt_avg,first,currently_training, &
navg,tlast,training_dxcall,trained_dxcall
if (first) then
navg=0
cross=cmplx(0.0,0.0)
cross_avg=cmplx(0.0,0.0)
wt_avg=0.0
tlast=0.0
trained_dxcall(1:12)=' '
training_dxcall(1:12)=' '
trained=.false.
currently_training=.false.
pcoeffs(1:3)=0.0
first=.false.
endif
if( (trained .and. (dxcall .ne. trained_dxcall)) .or. &
(currently_training .and. (dxcall .ne. training_dxcall)) .or. &
(navg .gt. 10 )) then !reset and retrain
navg=0
cross=cmplx(0.0,0.0)
cross_avg=cmplx(0.0,0.0)
wt_avg=0.0
tlast=0.0
trained_dxcall(1:12)=' '
trained=.false.
currently_training=.false.
training_dxcall(1:12)=' '
trained_dxcall(1:12)=' '
pcoeffs(1:3)=0.0
write(*,*) 'reset to untrained state '
endif
indx_dxcall=index(msg,trim(dxcall))
msg_has_dxcall = indx_dxcall .ge. 4
if( msg_has_dxcall .and. (.not. trained) .and. .not. currently_training ) then
currently_training=.true.
training_dxcall=trim(dxcall)
trained_dxcall(1:12)=' '
pcoeffs(1:3)=0.0
write(*,*) 'start training on call ',training_dxcall
endif
if( msg_has_dxcall .and. currently_training ) then
trained=.false. ! just to be sure
trained_dxcall(1:12)=' '
training_dxcall=dxcall
pcoeffs(1:3)=0.0
endif
! use decoded message to figure out how many bit errors in the frame
do i=1, 144
hardbits(i)=0
if(softbits(i) .gt. 0 ) hardbits(i)=1
enddo
! generate tones from decoded message
mygrid="EN50"
ichk=0
bcontest=.false.
call genmsk144(msg,mygrid,ichk,bcontest,msgsent,i4tone,itype)
! reconstruct message bits from tones
msgbits(1)=0
do i=1,143
if( i4tone(i) .eq. 0 ) then
if( mod(i,2) .eq. 1 ) then
msgbits(i+1)=msgbits(i)
else
msgbits(i+1)=mod(msgbits(i)+1,2)
endif
else
if( mod(i,2) .eq. 1 ) then
msgbits(i+1)=mod(msgbits(i)+1,2)
else
msgbits(i+1)=msgbits(i)
endif
endif
enddo
nbiterrors=0
do i=1,144
if( hardbits(i) .ne. msgbits(i) ) nbiterrors=nbiterrors+1
enddo
nplus=0
nminus=0
eyetop=1
eyebot=-1
do i=1,144
if( msgbits(i) .eq. 1 ) then
if( softbits(i) .lt. eyetop ) eyetop=softbits(i)
else
if( softbits(i) .gt. eyebot ) eyebot=softbits(i)
endif
enddo
eyeopening=eyetop-eyebot
is_training_frame = &
(snr.gt.5.0 .and.(nbiterrors.lt.7)) .and. &
(abs(t0-tlast) .gt. 0.072) .and. &
msg_has_dxcall
if( currently_training .and. is_training_frame ) then
twopi=8.0*atan(1.0)
nsym=144
if( i4tone(41) .lt. 0 ) nsym=40
dphi0=twopi*(freq-500)/12000.0
dphi1=twopi*(freq+500)/12000.0
phi=-twopi/8
indx=0
do i=1,nsym
if( i4tone(i) .eq. 0 ) then
dphi=dphi0
else
dphi=dphi1
endif
do j=1,6
waveform(indx)=cos(phi);
indx=indx+1
phi=mod(phi+dphi,twopi)
enddo
enddo
! convert the passband waveform to complex baseband
npts=864
nfft=1024
d=0
d(1:864)=waveform(0:863)
call analytic(d,npts,nfft,canalytic,pcoeffs,.false.,.false.) ! don't equalize the model
call tweak1(canalytic,nfft,-freq,cmodel)
call four2a(cframe(1:864),864,1,-1,1)
call four2a(cmodel(1:864),864,1,-1,1)
! Cross spectra from different messages can be averaged
! as long as all messages originate from dxcall.
cross=cmodel(1:864)*conjg(cframe)/1000.0
cross=cshift(cross,864/2)
cross_avg=cross_avg+10**(snr/20.0)*cross
wt_avg=wt_avg+10**(snr/20.0)
navg=navg+1
tlast=t0
phase=atan2(imag(cross_avg),real(cross_avg))
df=12000.0/864.0
nm=145
do i=1,145
x(i)=(i-73)*df/1000.0
enddo
y=phase((864/2-nm/2):(864/2+nm/2))
sigmay=wt_avg/abs(cross_avg((864/2-nm/2):(864/2+nm/2)))
mode=1
npts=145
nterms=5
call polyfit(x,y,sigmay,npts,nterms,mode,a,chisqr)
pp=a(1)+x*(a(2)+x*(a(3)+x*(a(4)+x*a(5))))
rmsdiff=sum( (pp-phase((864/2-nm/2):(864/2+nm/2)))**2 )/145.0
write(*,*) 'training ',navg,sqrt(chisqr),rmsdiff
! if( (sqrt(chisqr).lt.1.5) .and. (rmsdiff.lt.1.1) .and. (navg.ge.5) ) then
if( (sqrt(chisqr).lt.1.8) .and. (rmsdiff.lt.0.5) .and. (navg.ge.2) ) then
open(18,file='phasefit.dat',status='unknown')
do i=1, 145
write(18,*) x(i),pp(i),y(i),sigmay(i)
enddo
close(18)
open(19,file='cframe.dat',status='unknown')
do i=1,864
write(19,*) i,real(cframe(i)),imag(cframe(i)),real(cross_avg(i)),imag(cross_avg(i))
enddo
close(19)
pcoeffs=a(3:5)
trained_dxcall=dxcall
training_dxcall(1:12)=' '
currently_training=.false.
trained=.true.
call date_and_time(date,time,zone,values)
write(pcoeff_filename,'(i2.2,i2.2,i2.2,"_",i2.2,i2.2,i2.2)') &
values(1)-2000,values(2),values(3),values(5),values(6),values(7)
pcoeff_filename=trim(trained_dxcall)//"_"//pcoeff_filename//".pcoeff"
write(*,*) pcoeff_filename
open(17,file=pcoeff_filename,status='new')
write(17,'(i4,2f10.2,5f10.4)') navg,sqrt(chisqr),rmsdiff,a(1),a(2),a(3),a(4),a(5)
close(17)
endif
endif
return
end subroutine msk144signalquality