More work toward a usable reference spectrum.

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

lib/Makefile.jt9w

@@ -48,6 +48,10 @@ OBJS3 = t2.o four2a.o db.o
 t2:   $(OBJS3)
 	$(FC) -o t2 $(OBJS3) -lfftw3f
 
+OBJS4 = t3.o four2a.o db.o
+t3:   $(OBJS4)
+	$(FC) -o t3 $(OBJS4) -lfftw3f
+
 .PHONY : clean
 
 clean:

lib/Makefile.jt9w_Win

@@ -0,0 +1,58 @@
+
+# Set paths
+EXE_DIR = ..\\..\\wsjtx_install
+QT_DIR = C:/wsjt-env/Qt5/5.2.1/mingw48_32
+FFTW3_DIR = ..
+
+INCPATH = -I${QT_DIR}/include/QtCore -I${QT_DIR}/include 
+
+# Compilers
+CC = gcc
+CXX = g++
+FC = gfortran
+AR = ar cr
+RANLIB = ranlib
+MKDIR = mkdir -p
+CP = cp
+RM = rm -f
+
+FFLAGS = -O2 -fbounds-check -Wall -Wno-conversion
+CFLAGS = -O2 -I. 
+
+# Default rules
+%.o: %.c
+	${CC} ${CFLAGS} -c $<
+%.o: %.f
+	${FC} ${FFLAGS} -c $<
+%.o: %.F
+	${FC} ${FFLAGS} -c $<
+%.o: %.f90
+	${FC} ${FFLAGS} -c $<
+%.o: %.F90
+	${FC} ${FFLAGS} -c $<
+
+all:	 jt9w 
+
+OBJS1 = jt9w.o smo.o sync9w.o pctile.o shell.o lorentzian.o fchisq0.o \
+      softsym9w.o four2a.o interleave9.o jt9fano.o fano232.o packjt.o \
+      deg2grid.o grid2deg.o fmtmsg.o db.o decode9w.o
+      
+jt9w: $(OBJS1) 
+	$(FC) -o jt9w $(OBJS1) -lfftw3f
+
+OBJS2 = t1.o four2a.o db.o
+t1:   $(OBJS2)
+	$(FC) -o t1 $(OBJS2) -lfftw3f
+
+OBJS3 = t2.o four2a.o db.o
+t2:   $(OBJS3)
+	$(FC) -o t2 $(OBJS3) -lfftw3f
+
+OBJS4 = t3.o 
+t3:   $(OBJS4)
+	$(FC) -o t3 $(OBJS4) -L. -ljt9 C:\JTSDK\fftw3f\libfftw3f-3.dll
+
+.PHONY : clean
+
+clean:
+	$(RM) *.o JTMSKcode JTMSKcode.exe

lib/refspectrum.f90

@@ -12,7 +12,7 @@ subroutine refspectrum(id2,brefspec)
   complex cx(0:NH)
   equivalence(x,cx)
   data nsave/0/,brefspec0/.false./
-  save brefspec0,nsave
+  save brefspec0,nsave,s
 
   if(brefspec) then
      if(.not.brefspec0) then
@@ -32,19 +32,19 @@ subroutine refspectrum(id2,brefspec)
 
      if(mod(nsave,34).eq.0) then                   !About 9.8 sec 
         df=12000.0/NFFT
-        ia=nint(500.0/df)
+!        ia=nint(500.0/df)
-        ib=nint(2500.0/df)
+!        ib=nint(2500.0/df)
-        call pctile(s(ia),ib-ia+1,50,xmed)
+!        call pctile(s(ia),ib-ia+1,50,xmed)
-        db0=db(xmed)
+!        db0=db(xmed)
-        nhadd=10
+!        nhadd=10
         open(16,file='refspec.dat',status='unknown')
         do i=1,NH
            freq=i*df
-           ia=max(1,i-nhadd)
+!           ia=max(1,i-nhadd)
-           ib=min(NH,i+nhadd)
+!           ib=min(NH,i+nhadd)
-           smo=sum(s(ia:ib))/(ib-ia+1)
+!           smo=sum(s(ia:ib))/(ib-ia+1)
-           write(16,1000) freq,db(smo)-db0
+           write(16,1000) freq,s(i),db(s(i))
-1000       format(2f10.3)
+1000       format(f10.3,e12.3,f12.6)
         enddo
         close(16)
      endif

lib/t3.f90

@@ -0,0 +1,75 @@
+program t3
+
+  parameter (NBLK=3456,NZ=10*NBLK)
+  real x0(NZ)
+  real x1(NZ)
+
+  twopi=8.0*atan(1.0)
+  dphi=twopi*1000.0/12000.0
+  phi=0.
+  do i=1,NZ
+     phi=phi+dphi
+     x0(i)=sin(phi)
+     if(mod(i,10007).eq.100) x0(i)=2.0
+  enddo
+
+  do j=1,10
+     ib=j*NBLK
+     ia=ib-NBLK+1
+     call filter(x0(ia:ib),x1(ia:ib))
+  enddo
+
+  do i=1,NZ
+     write(13,1001) i,x0(i),x1(i),x1(i)-x0(i)
+1001 format(i6,3e12.3)
+  enddo
+
+end program t3
+
+subroutine filter(x0,x1)
+
+! Process time-domain data sequentially, optionally using 'refspec.dat' 
+! to flatten the spectrum.
+
+! NB: sin^2 window with 50% overlap; sin^2 + cos^2 = 1.0.
+
+  parameter (NFFT=6912,NH=NFFT/2)
+  real x0(0:NH-1)                         !Real input data
+  real x1(0:NH-1)                         !Real output data
+  real xov(0:NH-1)
+
+  real x(0:NFFT-1)
+  complex cx(0:NH)
+  real*4 w(0:NFFT-1)
+  real*4 s(0:NH)
+  logical first
+  equivalence (x,cx)
+  data first/.true./
+  save
+
+  if(first) then
+     pi=4.0*atan(1.0)
+     do i=0,NFFT-1
+        w(i)=(sin(i*pi/NFFT))**2
+     enddo
+     s=0.
+     fac=1.0/NFFT
+     first=.false.
+     xov=0.
+  endif
+
+  x(:NH-1)=xov                              !Previous 2nd half to new 1st half
+  x(NH:)=x0                                 !New 2nd half
+  x=x*w                                     !Apply window
+  call four2a(x,NFFT,1,-1,0)                !r2c FFT: to freq domain
+
+! Apply filter to cx()
+
+  call four2a(cx,NFFT,1,1,-1)               !c2r FFT: back to time domain
+
+  x(0:NH-1)=x(0:NH-1)+xov(0:NH-1)           !Add previous segment's 2nd half
+
+  x1=x
+
+  return
+end subroutine filter

mainwindow.cpp

@@ -874,7 +874,7 @@ void MainWindow::dataSink(qint64 frames)
   static float df3;
 
   int k (frames);
-  refspectrum_(&dec_data.d2[k-3456],&m_bRefSpec);
+  if(m_bRefSpec) refspectrum_(&dec_data.d2[k-3456],&m_bRefSpec);
 
   if(m_diskData) {
     dec_data.params.ndiskdat=1;