First attempt at a UI phase compensation tool for MSK144

This builds on the static phase compensation in the MSK144 decoder and
the phase  analysis and  polynomial fitting  capabilities also  in teh
MSK144 decoder,  by allowing  captured data to  be selected  for phase
equalization from the WSJT-X UI.

Reads  captured phase  compensation  curve  estimate files  containing
fitted  polynomial coefficients  and measured  phase data  from MSK144
receptions. Intent  is to select a  compensation curve that is  from a
known transmitter like  an SDR which have good  phase linearity. Phase
plots and compensation polynomials may be viewed and compared with the
current compensation polynomial. A  suitable polynomial can be applied
to be use in all further decoding of MSK144 signals.

Plots of  the currently  selected polynomial  and its  modified higher
order terms  polynomial which is  actually used in  equalization (this
plot may  be dropped - it  is just for  kicks at the moment).   When a
captured phase analysis file is loaded plots of the measured phase and
the proposed best fit polynomial are shown.

Basic  maintenance  is also  included  allowing  clearing and  loading
captured  plots and  an option  to revert  to a  flat no  equalization
curve.

More to come on this as  amplitude equalization is also possible, this
will probably  be similar, maybe even  plotted on the same  graph with
dual  axes  for phase  and  amplitude.   Amplitude correction  from  a
measured  reference   spectrum  could  be  viewed   and  selected  for
equalization for all modes. TBC...

This   change    also   introduces    the   QCustomPlot    3rd   party
widget. Currently this  is statically linked from a  qcp library built
by the WSJT-X CMake script. This will probably be migrated to a shared
object (DLL) build as a CMake external project, once some CMake script
re-factoring  has been  completed,  which  is more  in  line with  the
QCustomPlot author's  intentions. This  will allow efficient  reuse in
other tools shipped with WSJT-X.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7570 ab8295b8-cf94-4d9e-aec4-7959e3be5d79

lib/analytic.f90

@@ -6,13 +6,12 @@ subroutine analytic(d,npts,nfft,c,dpc,bseq,bdeq)
 
   real d(npts)              ! passband signal
   real h(NFFTMAX/2)         ! real BPF magnitude
-  real dpc(3),dpclast(3)    ! dynamic phase coeffs
-  real spc(3),spclast(3)    ! static phase coeffs
+  real dpc(5)               ! dynamic phase coeffs
   real sac(5),saclast(5)    ! amp coeffs
   real fp     
 
   complex corrs(NFFTMAX/2)  ! allpass static phase correction 
-  complex corrd(NFFTMAX/2)   ! allpass overall phase correction
+  complex corrd(NFFTMAX/2)  ! allpass overall phase correction
   complex c(NFFTMAX)        ! analytic signal
 
   logical*1 bseq            ! boolean static equalizer flag
@@ -21,12 +20,10 @@ subroutine analytic(d,npts,nfft,c,dpc,bseq,bdeq)
 
   data nfft0/0/
   data bseqlast/.false./
-  data spclast/0.0,0.0,0.0/
   data saclast/1.0,0.0,0.0,0.0,0.0/
-  data spc/-0.952,0.768,-0.565/  ! baseline phase coeffs for TS2000
   data sac/1.0,0.05532,0.11438,0.12918,0.09274/ ! amp coeffs for TS2000
 
-  save corrs,corrd,nfft0,h,sac,spc,saclast,spclast,dpclast,pi,t,beta
+  save corrs,corrd,nfft0,h,sac,saclast,pi,t,beta
 
   df=12000.0/nfft
   nh=nfft/2
@@ -47,18 +44,14 @@ subroutine analytic(d,npts,nfft,c,dpc,bseq,bdeq)
      nfft0=nfft
   endif
 
-  if( any(spclast .ne. spc) .or. any(saclast .ne. sac) .or. any(dpclast .ne. dpc) ) then
-     spclast=spc
-     dpclast=dpc
+  if( any(saclast .ne. sac) ) then
      saclast=sac
      do i=1,nh+1
         ff=(i-1)*df
         f=ff-1500.0
         fp=f/1000.0
-        ps=fp*fp*(spc(1)+fp*(spc(2)+fp*spc(3)))
-        amp=sac(1)+fp*(sac(2)+fp*(sac(3)+fp*(sac(4)+fp*sac(5))))
-        corrs(i)=amp*cmplx(cos(ps),sin(ps))
-        pd=fp*fp*(dpc(1)+fp*(dpc(2)+fp*dpc(3)))
+        corrs(i)=sac(1)+fp*(sac(2)+fp*(sac(3)+fp*(sac(4)+fp*sac(5))))
+        pd=fp*fp*(dpc(3)+fp*(dpc(4)+fp*dpc(5))) ! ignore 1st two terms
         corrd(i)=cmplx(cos(pd),sin(pd))
      enddo
   endif

lib/hspec.f90

@@ -1,5 +1,6 @@
 subroutine hspec(id2,k,nutc0,ntrpdepth,nrxfreq,ntol,bmsk144,bcontest,  &
-     brxequal,btrain,ingain,mycall,hiscall,bshmsg,bswl,datadir,green,s,jh,line1,mygrid)
+     btrain,pcoeffs,ingain,mycall,hiscall,bshmsg,bswl,datadir,green,s,jh,line1, &
+     mygrid)
 
 ! Input:
 !  k         pointer to the most recent new data
@@ -8,7 +9,6 @@ subroutine hspec(id2,k,nutc0,ntrpdepth,nrxfreq,ntol,bmsk144,bcontest,  &
 !  nrxfreq   Rx audio center frequency
 !  ntol      Decoding range is +/- ntol
 !  bmsk144   Boolean, true if in MSK144 mode
-!  brxequal  Boolean, turns on equalization in MSK144 mode
 !  btrain    Boolean, turns on training in MSK144 mode
 !  ingain    Relative gain for spectra
 
@@ -23,10 +23,11 @@ subroutine hspec(id2,k,nutc0,ntrpdepth,nrxfreq,ntol,bmsk144,bcontest,  &
   character*12 mycall,hiscall
   character*6 mygrid
   integer*2 id2(0:120*12000-1)
-  logical*1 bmsk144,bcontest,bshmsg,brxequal,btrain,bswl
+  logical*1 bmsk144,bcontest,bshmsg,btrain,bswl
   real green(0:JZ-1)
   real s(0:63,0:JZ-1)
   real x(512)
+  real pcoeffs(5)
   complex cx(0:256)
   data rms/999.0/,k0/99999999/
   equivalence (x,cx)
@@ -86,7 +87,8 @@ subroutine hspec(id2,k,nutc0,ntrpdepth,nrxfreq,ntol,bmsk144,bcontest,  &
         tt2=sum(float(abs(id2(k0:k0+3583))))
         if(tt1.ne.0.0 .and. tt2.ne.0) then
            call mskrtd(id2(k-7168+1:k),nutc0,tsec,ntol,nrxfreq,ndepth,   &
-                mycall,mygrid,hiscall,bshmsg,bcontest,brxequal,btrain,bswl,datadir,line1)
+                mycall,mygrid,hiscall,bshmsg,bcontest,btrain,pcoeffs,bswl,&
+                datadir,line1)
         endif
      endif
   endif

lib/msk144signalquality.f90

@@ -1,4 +1,4 @@
-  subroutine msk144signalquality(cframe,snr,freq,t0,softbits,msg,dxcall,       &
+subroutine msk144signalquality(cframe,snr,freq,t0,softbits,msg,dxcall,       &
    btrain,datadir,nbiterrors,eyeopening,pcoeffs)
 
   character*22 msg,msgsent
@@ -36,7 +36,7 @@
   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)
+  real pcoeffs(5)
 
   data first/.true./
   save cross_avg,wt_avg,first,currently_training,   &
@@ -51,7 +51,6 @@
     trained_dxcall(1:12)=' '
     training_dxcall(1:12)=' '
     currently_training=.false.
-    pcoeffs(1:3)=0.0
     first=.false.
   endif
 
@@ -66,7 +65,6 @@
     currently_training=.false.
     training_dxcall(1:12)=' '
     trained_dxcall(1:12)=' '
-    pcoeffs(1:3)=0.0
 write(*,*) 'reset to untrained state '
   endif
 
@@ -77,14 +75,12 @@ write(*,*) 'reset to untrained state '
     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_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 
@@ -193,7 +189,6 @@ write(*,*) 'start training on call ',training_dxcall
       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.8) .and. (rmsdiff.lt.0.5) .and. (navg.ge.2) ) then 
-!        pcoeffs=a(3:5)
         trained_dxcall=dxcall
         training_dxcall(1:12)=' '
         currently_training=.false.

lib/mskrtd.f90

@@ -1,5 +1,5 @@
 subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall,   &
-     bshmsg,bcontest,brxequal,btrain,bswl,datadir,line)
+     bshmsg,bcontest,btrain,pcoeffs,bswl,datadir,line)
 
 ! Real-time decoder for MSK144.  
 ! Analysis block size = NZ = 7168 samples, t_block = 0.597333 s 
@@ -37,9 +37,9 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall,   &
   real pow(8)
   real softbits(144)
   real xmc(NPATTERNS)
-  real pcoeffs(3)
+  real pcoeffs(5)
 
-  logical*1 bshmsg,bcontest,brxequal,btrain,bswl
+  logical*1 bshmsg,bcontest,btrain,bswl
   logical*1 first
   logical*1 bshdecode
   logical*1 seenb4
@@ -53,14 +53,13 @@ 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,cdat,pcoeffs,msglast,msglastswl,     &
+  save first,tsec0,nutc00,pnoise,cdat,msglast,msglastswl,     &
        nsnrlast,nsnrlastswl,recent_calls,nhasharray,recent_shmsgs
 
   if(first) then
      tsec0=tsec
      nutc00=nutc0
      pnoise=-1.0
-     pcoeffs(1:3)=0.0
      do i=1,nrecent
        recent_calls(i)(1:12)=' '
      enddo
@@ -96,7 +95,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.
-  bvar=brxequal
+  bvar=.true.
   if( btrain ) bvar=.false.   ! if training, turn off rx eq
   call analytic(d,NZ,NFFT1,cdat,pcoeffs,bvar,.false.)  ! never apply dynamic coeffs
 
@@ -191,8 +190,8 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall,   &
                           btrain,datadir,ncorrected,eyeopening,pcoeffs)
   endif
 
-  decsym=' & '
-  if( brxequal ) decsym=' ^ '
+  decsym=' ^ '
+  if( btrain ) decsym=' & '
   if( msgreceived(1:1).eq.'<') then
     ncorrected=0
     eyeopening=0.0