Add an option on advanced time to enable MSK144 Rx filter equalizer. Other changes necessary to accommodate coexistence of static (Rx Filter) and dynamic (QSO partner-dependent) phase corrections.

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

CMakeLists.txt

@@ -421,6 +421,7 @@ set (wsjt_FSRCS
   lib/msk40_freq_search.f90
   lib/msk144_freq_search.f90
   lib/mskrtd.f90
+  lib/msk144signalquality.f90
   lib/msk144sim.f90
   lib/mskrtd.f90
   lib/options.f90

Configuration.cpp

@@ -540,6 +540,7 @@ private:
   bool twoPass_;
   bool x2ToneSpacing_;
   bool contestMode_;
+  bool rxEqualize_;
   bool realTimeDecode_;
   QString udp_server_name_;
   port_type udp_server_port_;
@@ -628,6 +629,7 @@ bool Configuration::single_decode () const {return m_->single_decode_;}
 bool Configuration::twoPass() const {return m_->twoPass_;}
 bool Configuration::x2ToneSpacing() const {return m_->x2ToneSpacing_;}
 bool Configuration::contestMode() const {return m_->contestMode_;}
+bool Configuration::rxEqualize() const {return m_->rxEqualize_;}
 bool Configuration::realTimeDecode() const {return m_->realTimeDecode_;}
 bool Configuration::split_mode () const {return m_->split_mode ();}
 QString Configuration::udp_server_name () const {return m_->udp_server_name_;}
@@ -1057,6 +1059,7 @@ void Configuration::impl::initialize_models ()
   ui_->cbTwoPass->setChecked(twoPass_);
   ui_->cbx2ToneSpacing->setChecked(x2ToneSpacing_);
   ui_->cbContestMode->setChecked(contestMode_);
+  ui_->cbRxEqualize->setChecked(rxEqualize_);
   ui_->cbRealTime->setChecked(realTimeDecode_);
   ui_->cbRealTime->setVisible(false);                    //Tempoary -- probably will remove this control
   ui_->type_2_msg_gen_combo_box->setCurrentIndex (type_2_msg_gen_);
@@ -1285,6 +1288,7 @@ void Configuration::impl::read_settings ()
   twoPass_ = settings_->value("TwoPass",true).toBool ();
   x2ToneSpacing_ = settings_->value("x2ToneSpacing",false).toBool ();
   contestMode_ = settings_->value("ContestMode",false).toBool ();
+  rxEqualize_ = settings_->value("RxEqualize",false).toBool ();
   realTimeDecode_ = settings_->value("RealTimeDecode",false).toBool ();
   rig_params_.poll_interval = settings_->value ("Polling", 0).toInt ();
   rig_params_.split_mode = settings_->value ("SplitMode", QVariant::fromValue (TransceiverFactory::split_mode_none)).value<TransceiverFactory::SplitMode> ();
@@ -1384,6 +1388,7 @@ void Configuration::impl::write_settings ()
   settings_->setValue ("TwoPass", twoPass_);
   settings_->setValue ("x2ToneSpacing", x2ToneSpacing_);
   settings_->setValue ("ContestMode", contestMode_);
+  settings_->setValue ("RxEqualize", rxEqualize_);
   settings_->setValue ("RealTimeDecode", realTimeDecode_);
   settings_->setValue ("UDPServer", udp_server_name_);
   settings_->setValue ("UDPServerPort", udp_server_port_);
@@ -1779,6 +1784,7 @@ void Configuration::impl::accept ()
   twoPass_ = ui_->cbTwoPass->isChecked ();
   x2ToneSpacing_ = ui_->cbx2ToneSpacing->isChecked ();
   contestMode_ = ui_->cbContestMode->isChecked ();
+  rxEqualize_ = ui_->cbRxEqualize->isChecked ();
   realTimeDecode_ = ui_->cbRealTime->isChecked ();
   frequency_calibration_intercept_ = ui_->calibration_intercept_spin_box->value ();
   frequency_calibration_slope_ppm_ = ui_->calibration_slope_ppm_spin_box->value ();

Configuration.hpp

@@ -123,6 +123,7 @@ public:
   bool twoPass() const;
   bool x2ToneSpacing() const;
   bool contestMode() const;
+  bool rxEqualize() const;
   bool realTimeDecode() const;
   bool MyDx() const;
   bool CQMyN() const;

Configuration.ui

@@ -2295,6 +2295,16 @@ Right click for insert and delete options.</string>
            </widget>
           </item>
           <item row="8" column="0">
+           <widget class="QCheckBox" name="cbRxEqualize">
+            <property name="toolTip">
+             <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;Correct for Rx filter group-delay variation.&lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
+            </property>
+            <property name="text">
+             <string>MSK144 Rx Equalization</string>
+            </property>
+           </widget>
+          </item>
+          <item row="9" column="0">
            <widget class="QCheckBox" name="cbx2ToneSpacing">
             <property name="toolTip">
              <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;Generate Tx audio with twice the normal tone spacing.  Intended for special LF/MF transmitters that use a divide-by-2 before generating RF.&lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
@@ -2511,6 +2521,7 @@ soundcard changes</string>
   <tabstop>sbBandwidth</tabstop>
   <tabstop>sbTxDelay</tabstop>
   <tabstop>cbContestMode</tabstop>
+  <tabstop>cbRxEqualize</tabstop>
   <tabstop>cbx2ToneSpacing</tabstop>
   <tabstop>cbRealTime</tabstop>
  </tabstops>

lib/analytic.f90

@@ -1,41 +1,42 @@
-subroutine analytic(d,npts,nfft,c,a,equalize)
+subroutine analytic(d,npts,nfft,c,dpc,bseq,bdeq)
 
 ! Convert real data to analytic signal
 
   parameter (NFFTMAX=1024*1024)
-  real d(npts)
-  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/
-  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
+  real d(npts)              ! passband signal
-  a0(1:5)=0.0
+  real h(NFFTMAX/2)         ! real BPF magnitude
+  real dpc(3),dpclast(3)    ! dynamic phase coeffs
+  real spc(3),spclast(3)    ! static phase coeffs
+  real ac(5)                ! currently unused static amp coeffs
+  real fp     
+
+  complex corrs(NFFTMAX/2)  ! allpass static phase correction 
+  complex corrd(NFFTMAX/2)   ! allpass overall phase correction
+  complex c(NFFTMAX)        ! analytic signal
+
+  logical*1 bseq            ! boolean static equalizer flag
+  logical*1 bdeq            ! boolean dynamic equalizer flag
+  logical*1 bseqlast      
+
+  data nfft0/0/
+  data bseqlast/.false./
+  data spclast/0.0,0.0,0.0/
+  data spc/-0.952,0.768,-0.565/  ! baseline phase coeffs for TS2000
+  data ac/1.0,0.05532,0.11438,0.12918,0.09274/ ! amp coeffs for TS2000
+
+  save nfft0,h,spc,spclast,dpclast,ac,pi,t,beta
 
   df=12000.0/nfft
   nh=nfft/2
-  if( nfft.ne.nfft0 .or. any(alast .ne. a) ) then
+  if( nfft.ne.nfft0 ) then
+     pi=4.0*atan(1.0)
      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
-        fp=f/1000.0
+        h(i)=1.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
@@ -43,15 +44,39 @@ subroutine analytic(d,npts,nfft,c,a,equalize)
      nfft0=nfft
   endif
 
+  if( any(spclast .ne. spc) .or. any(dpclast .ne. dpc) ) then
+     spclast=spc
+     dpclast=dpc
+     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)))
+        corrs(i)=cmplx(cos(ps),sin(ps))
+        pd=fp*fp*(dpc(1)+fp*(dpc(2)+fp*dpc(3)))
+        corrd(i)=cmplx(cos(pd),sin(pd))
+     enddo
+  endif
+
   fac=2.0/nfft
   c(1:npts)=fac*d(1:npts)
   c(npts+1:nfft)=0.
   call four2a(c,nfft,1,-1,1)               !Forward c2c FFT
 
+  if( (.not. bseq) .and. (.not. bdeq) ) then
     c(1:nh+1)=h(1:nh+1)*c(1:nh+1)
+  else if( bseq .and. (.not. bdeq) ) then
+    c(1:nh+1)=h(1:nh+1)*corrs(1:nh+1)*c(1:nh+1)
+  else if( (.not. bseq) .and. bdeq ) then
+    c(1:nh+1)=h(1:nh+1)*corrd(1:nh+1)*c(1:nh+1)
+  else if( bseq .and. bdeq ) then
+    c(1:nh+1)=h(1:nh+1)*corrs(1:nh+1)*corrd(1:nh+1)*c(1:nh+1)
+  endif
+
   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
+  spclast=spc
+  dpclast=dpc
   return
 end subroutine analytic

lib/hspec.f90

@@ -1,5 +1,5 @@
 subroutine hspec(id2,k,nutc0,ntrpdepth,nrxfreq,ntol,bmsk144,bcontest,  &
-     ingain,mycall,hiscall,bshmsg,green,s,jh,line1,mygrid)
+     brxequal,ingain,mycall,hiscall,bshmsg,green,s,jh,line1,mygrid)
 
 ! Input:
 !  k         pointer to the most recent new data
@@ -8,6 +8,7 @@ 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
 !  ingain    Relative gain for spectra
 
 ! Output:
@@ -20,7 +21,7 @@ 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
+  logical*1 bmsk144,bcontest,bshmsg,brxequal
   real green(0:JZ-1)
   real s(0:63,0:JZ-1)
   real x(512)
@@ -83,7 +84,7 @@ 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,line1)
+                mycall,mygrid,hiscall,bshmsg,bcontest,brxequal,line1)
         endif
      endif
   endif

lib/msk144signalquality.f90

@@ -1,4 +1,5 @@
-subroutine msk144signalquality(cframe,snr,freq,t0,softbits,msg,dxcall,nbiterrors,eyeopening,trained,pcoeffs)
+  subroutine msk144signalquality(cframe,snr,freq,t0,softbits,msg,dxcall,  &
+       nbiterrors,eyeopening,trained,pcoeffs)
 
   character*22 msg,msgsent
   character*12 dxcall
@@ -35,7 +36,8 @@ subroutine msk144signalquality(cframe,snr,freq,t0,softbits,msg,dxcall,nbiterrors
   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,pcoeffs(5)
+  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,   &
@@ -51,7 +53,7 @@ subroutine msk144signalquality(cframe,snr,freq,t0,softbits,msg,dxcall,nbiterrors
     training_dxcall(1:12)=' '
     trained=.false.
     currently_training=.false.
-    pcoeffs(1:5)=0.0
+    pcoeffs(1:3)=0.0
     first=.false.
   endif
 
@@ -68,7 +70,7 @@ subroutine msk144signalquality(cframe,snr,freq,t0,softbits,msg,dxcall,nbiterrors
     currently_training=.false.
     training_dxcall(1:12)=' '
     trained_dxcall(1:12)=' '
-    pcoeffs(1:5)=0.0
+    pcoeffs(1:3)=0.0
 write(*,*) 'reset to untrained state '
   endif
 
@@ -79,7 +81,7 @@ write(*,*) 'reset to untrained state '
     currently_training=.true.
     training_dxcall=trim(dxcall)
     trained_dxcall(1:12)=' '
-    pcoeffs(1:5)=0.0
+    pcoeffs(1:3)=0.0
 write(*,*) 'start training on call ',training_dxcall
   endif
 
@@ -87,7 +89,7 @@ write(*,*) 'start training on call ',training_dxcall
     trained=.false. ! just to be sure
     trained_dxcall(1:12)=' '
     training_dxcall=dxcall
-    pcoeffs(1:5)=0.0
+    pcoeffs(1:3)=0.0
   endif
 
 ! use decoded message to figure out how many bit errors in the frame 
@@ -167,8 +169,7 @@ write(*,*) 'start training on call ',training_dxcall
       nfft=1024
       d=0
       d(1:864)=waveform(0:863)
-      a=0
+      call analytic(d,npts,nfft,canalytic,pcoeffs,.false.,.false.) ! don't equalize the model
-      call analytic(d,npts,nfft,canalytic,a,.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)
@@ -207,7 +208,7 @@ write(*,*) 'training ',navg,sqrt(chisqr),rmsdiff
           write(19,*) i,real(cframe(i)),imag(cframe(i))
         enddo
         close(19)
-        pcoeffs=a
+        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,line)
+     bshmsg,bcontest,brxequal,line)
 
 ! Real-time decoder for MSK144.  
 ! Analysis block size = NZ = 7168 samples, t_block = 0.597333 s 
@@ -30,11 +30,11 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall,   &
   real pow(8)
   real softbits(144)
   real xmc(NPATTERNS)
-  real*8 pcoeffs(5)
+  real pcoeffs(3)
 
-  logical*1 bshmsg,bcontest
+  logical*1 bshmsg,bcontest,brxequal
   logical first
-  logical*1 equalized
+  logical*1 trained 
 
   data first/.true./
   data iavpatterns/ &
@@ -43,14 +43,14 @@ 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,pcoeffs,equalized
+  save first,tsec0,nutc00,pnoise,nsnrlast,msglast,cdat,pcoeffs,trained
 
   if(first) then
      tsec0=tsec
      nutc00=nutc0
      pnoise=-1.0
-     pcoeffs(1:5)=0.0
+     pcoeffs(1:3)=0.0
-     equalized=.false.
+     trained=.false.
      first=.false.
   endif
 
@@ -74,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,pcoeffs,.true.) !Convert to analytic signal and filter
+  call analytic(d,NZ,NFFT1,cdat,pcoeffs,brxequal,.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.
@@ -101,8 +101,6 @@ 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
@@ -125,7 +123,7 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall,   &
      if( nsyncsuccess .eq. 0 ) cycle
 
      do ipk=1,npeaks
-        do is=1,3    ! With equalization, this loop may not be necessary
+        do is=1,3   
            ic0=npkloc(ipk)
            if(is.eq.2) ic0=max(1,ic0-1)
            if(is.eq.3) ic0=min(NSPM,ic0+1)
@@ -133,9 +131,6 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall,   &
            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
@@ -165,7 +160,7 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall,   &
   nsnr=nint(snr0)
 
   call msk144signalquality(ct,snr0,fest,tdec,softbits,msgreceived,hiscall,   &
-                           ncorrected,eyeopening,equalized,pcoeffs)
+                           ncorrected,eyeopening,trained,pcoeffs)
 
 ! Dupe check. Only print if new message, or higher snr.
   if(msgreceived.ne.msglast .or. nsnr.gt.nsnrlast .or. tsec.lt.tsec0) then
@@ -177,6 +172,10 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall,   &
      if(msgreceived(1:1).ne.'<') then
         call fix_contest_msg(mycall(1:6),mygrid,hiscall(1:6),msgreceived)
      endif
+     decsym=' & '
+     if( brxequal .and. (.not. trained) ) decsym=' ^ '
+     if( brxequal .and. trained ) decsym=' $ '
+     if( (.not. brxequal) .and. trained ) decsym=' @ '
      write(line,1020) nutc0,nsnr,tdec,nint(fest),decsym,msgreceived,    &
           navg,ncorrected,eyeopening,char(0)
 1020 format(i6.6,i4,f5.1,i5,a3,a22,i2,i3,f5.1,a1)
@@ -187,5 +186,3 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall,   &
 end subroutine mskrtd
 
 include 'fix_contest_msg.f90'
-
-include 'msk144signalquality.f90'

mainwindow.cpp

@@ -67,9 +67,9 @@ extern "C" {
                 int* minw, float* px, float s[], float* df3, int* nhsym, int* npts8);
 
   void hspec_(short int d2[], int* k, int* nutc0, int* ntrperiod, int* nrxfreq, int* ntol,
-              bool* bmsk144, bool* bcontest, int* ingain, char mycall[], char hiscall[],
+              bool* bmsk144, bool* bcontest, bool* brxequalize, int* ingain, char mycall[], 
-              bool* bshmsg, float green[], float s[], int* jh, char line[],
+              char hiscall[], bool* bshmsg, float green[], float s[], int* jh,
-              char mygrid[], int len1, int len2, int len3, int len4);
+              char line[], char mygrid[], int len1, int len2, int len3, int len4);
 
   void gen4_(char* msg, int* ichk, char* msgsent, int itone[],
                int* itext, int len1, int len2);
@@ -1284,9 +1284,10 @@ void MainWindow::fastSink(qint64 frames)
   QString hisCall {ui->dxCallEntry->text ()};
   bool bshmsg=ui->cbShMsgs->isChecked();
   bool bcontest=m_config.contestMode();
+  bool brxequalize=m_config.rxEqualize();
   strncpy(dec_data.params.hiscall,(hisCall + "            ").toLatin1 ().constData (), 12);
   strncpy(dec_data.params.mygrid, (m_config.my_grid()+"      ").toLatin1(),6);
-  hspec_(dec_data.d2,&k,&nutc0,&nTRpDepth,&m_RxFreq,&m_Ftol,&bmsk144,&bcontest,
+  hspec_(dec_data.d2,&k,&nutc0,&nTRpDepth,&m_RxFreq,&m_Ftol,&bmsk144,&bcontest,&brxequalize,
          &m_inGain,&dec_data.params.mycall[0],&dec_data.params.hiscall[0],&bshmsg,
          fast_green,fast_s,&fast_jh,&line[0],&dec_data.params.mygrid[0],12,12,80,6);
   float px = fast_green[fast_jh];