Add reference spectrum to equalization plots and more plotting enhancements

Use  a header  format for  polynomial coefficients  that includes  the
valid  X  range  in  scaled  terms  and  a  count  of  the  number  of
coefficients.

Use double  precision consistently  for polynomial  coefficients. This
includes formatting with sufficient DPs when writing to files.

Many changes to the equalization plots, more to come.

Add  error   handling  for   reading  coefficient,  plot   and  filter
files.  This  includes  being  backward  compatible  for  old  format
refspec.dat files with no header.

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

lib/analytic.f90

@@ -6,9 +6,9 @@ subroutine analytic(d,npts,nfft,c,pc,beq)
 
   real d(npts)              ! passband signal
   real h(NFFTMAX/2)         ! real BPF magnitude
-  real pc(5),pclast(5)      ! static phase coeffs
-  real ac(5),aclast(5)      ! amp coeffs
-  real fp     
+  real*8 pc(5),pclast(5)    ! static phase coeffs
+  real*8 ac(5),aclast(5)    ! amp coeffs
+  real*8 fp
 
   complex corr(NFFTMAX/2)  ! complex frequency-dependent correction 
   complex c(NFFTMAX)        ! analytic signal

lib/hspec.f90

@@ -27,7 +27,7 @@ subroutine hspec(id2,k,nutc0,ntrpdepth,nrxfreq,ntol,bmsk144,bcontest,  &
   real green(0:JZ-1)
   real s(0:63,0:JZ-1)
   real x(512)
-  real pcoeffs(5)
+  real*8 pcoeffs(5)
   complex cx(0:256)
   data rms/999.0/,k0/99999999/
   equivalence (x,cx)

lib/msk144signalquality.f90

@@ -36,8 +36,9 @@ subroutine msk144signalquality(cframe,snr,freq,t0,softbits,msg,dxcall,       &
   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(5)
+  real*8 pcoeffs(5)
 
+  parameter (NFREQLOW=500,NFREQHIGH=2500)
   data first/.true./
   save cross_avg,wt_avg,first,currently_training,   &
        navg,tlast,training_dxcall,trained_dxcall
@@ -202,7 +203,7 @@ write(*,*) 'training ',navg,sqrt(chisqr),rmsdiff
         pcoeff_filename=datadir(1:l1+1)//trim(pcoeff_filename)
 write(*,*) 'trained - writing coefficients to: ',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)
+        write(17,'(i4,2f10.2,3i5,5e25.16)') navg,sqrt(chisqr),rmsdiff,NFREQLOW,NFREQHIGH,nterms,a
         do i=1, 145
           write(17,*) x(i),pp(i),y(i),sigmay(i)
         enddo

lib/mskrtd.f90

@@ -37,7 +37,7 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,mycall,mygrid,hiscall,   &
   real pow(8)
   real softbits(144)
   real xmc(NPATTERNS)
-  real pcoeffs(5)
+  real*8 pcoeffs(5)
 
   logical*1 bshmsg,bcontest,btrain,bswl
   logical*1 first

lib/refspectrum.f90

@@ -4,9 +4,9 @@ subroutine refspectrum(id2,bclear,brefspec,buseref,fname)
 !  id2       i*2        Raw 16-bit integer data, 12000 Hz sample rate
 !  brefspec  logical    True when accumulating a reference spectrum
 
-  parameter (NFFT=6912,NH=NFFT/2)
+  parameter (NFFT=6912,NH=NFFT/2,NPOLYLOW=400,NPOLYHIGH=2600)
   integer*2 id2(NFFT)
-  logical*1 bclear,brefspec,buseref
+  logical*1 bclear,brefspec,buseref,blastuse
   
   real x0(0:NH-1)                         !Input samples
   real x1(0:NH-1)                         !Output samples (delayed by one block)
@@ -17,12 +17,12 @@ subroutine refspectrum(id2,bclear,brefspec,buseref,fname)
   real*4 s(0:NH)                          !Average spectrum
   real*4 fil(0:NH)
   real*8 xfit(1500),yfit(1500),sigmay(1500),a(5),chisqr !Polyfit arrays
-  logical first,firstuse
+  logical first
   complex cx(0:NH)                        !Complex frequency-domain work array
   character*(*) fname
   common/spectra/syellow(6827),ref(0:NH),filter(0:NH)
   equivalence(x,cx)
-  data first/.true./,firstuse/.true./
+  data first/.true./,blastuse/.false./
   save
 
   if(first) then
@@ -98,8 +98,8 @@ subroutine refspectrum(id2,bclear,brefspec,buseref,fname)
            if(s(i).gt.0.0) filter(i)=20.0*log10(fil(i))
         enddo
 
-        il=nint(400.0/df)
-        ih=nint(2600.0/df)
+        il=nint(NPOLYLOW/df)
+        ih=nint(NPOLYHIGH/df)
         nfit=ih-il+1
         mode=0
         nterms=5
@@ -111,8 +111,8 @@ subroutine refspectrum(id2,bclear,brefspec,buseref,fname)
         call polyfit(xfit,yfit,sigmay,nfit,nterms,mode,a,chisqr)
 
         open(16,file=fname,status='unknown')
-        write(16,1003) a
-1003    format(5f10.4)
+        write(16,1003) NPOLYLOW,NPOLYHIGH,nterms,a
+1003    format(3i5,5e25.16)
         do i=1,NH
            freq=i*df
            ref(i)=db(s(i)/avemid)
@@ -125,17 +125,18 @@ subroutine refspectrum(id2,bclear,brefspec,buseref,fname)
   endif
 
   if(buseref) then
-     if(firstuse) then
+     if(blastuse.neqv.buseref) then !just enabled so read filter
         fil=1.0
-        open(16,file=fname,status='old',err=10)
-        read(16,1003,err=10,end=10) a
-        do i=1,NH
-           read(16,1005,err=10,end=10) freq,s(i),ref(i),fil(i),filter(i)
+        open(16,file=fname,status='old',err=110)
+        read(16,1003,err=20,end=100) ndummy,ndummy,nterms,a
+        goto 30
+20      rewind(16)              !allow for old style refspec.dat with no header
+30      do i=1,NH
+           read(16,1005,err=100,end=100) freq,s(i),ref(i),fil(i),filter(i)
         enddo
-10      close(16)
-      firstuse=.false.
+100     close(16)
+110     continue
      endif
-!     x0=id2(NH+1:NFFT)
      x0=id2(1:NH)
      x(0:NH-1)=x0s                         !Previous 2nd half to new 1st half
      x(NH:NFFT-1)=x0                       !New 2nd half
@@ -148,6 +149,6 @@ subroutine refspectrum(id2,bclear,brefspec,buseref,fname)
      id2(1:NH)=nint(x1)
      x1s=x(NH:NFFT-1)                      !Save the new 2nd half
   endif
-
+  blastuse=buseref
   return
 end subroutine refspectrum