WSJT-X/lib/msk144signalquality.f90

  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.5) .and. (rmsdiff.lt.1.1) .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