From 3541d12a6735c1581d6e8838ed8d3e8a24764021 Mon Sep 17 00:00:00 2001
From: Joe Taylor <k1jt@arrl.org>
Date: Mon, 11 Apr 2016 15:24:55 +0000
Subject: [PATCH] Test program

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@6607 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
---
 lib/t3.f90 | 59 +++++++++++++++++++++++++++---------------------------
 1 file changed, 30 insertions(+), 29 deletions(-)

diff --git a/lib/t3.f90 b/lib/t3.f90
index af818061d..82050f932 100644
--- a/lib/t3.f90
+++ b/lib/t3.f90
@@ -19,30 +19,32 @@ program t3
      call filter(x0(ia:ib),x1(ia:ib))
   enddo
 
-  do i=1,NZ
+  x1(1:NZ-NBLK)=x1(NBLK+1:NZ)
+  do i=1,NZ-NBLK
      write(13,1001) i,x0(i),x1(i),x1(i)-x0(i)
-1001 format(i6,3e12.3)
+1001 format(i6,3f13.9)
   enddo
 
 end program t3
 
 subroutine filter(x0,x1)
 
-! Process time-domain data sequentially, optionally using 'refspec.dat' 
-! to flatten the spectrum.
+! Process time-domain data sequentially, optionally using a frequency-domain
+! filter to alter the spectrum.
 
-! NB: sin^2 window with 50% overlap; sin^2 + cos^2 = 1.0.
+! NB: uses a sin^2 window with 50% overlap.
 
   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)
+  real x0(0:NH-1)                         !Input samples
+  real x1(0:NH-1)                         !Output samples (delayed by one block)
+  real x0s(0:NH-1)                        !Saved upper half of input samples
+  real x1s(0:NH-1)                        !Saved upper half of output samples
+  real x(0:NFFT-1)                        !Work array
+  real*4 w(0:NFFT-1)                      !Window function
+  real f(0:NH)                            !Filter to be applied
+  real*4 s(0:NH)                          !Average spectrum
   logical first
+  complex cx(0:NH)                        !Complex frequency-domain work array
   equivalence (x,cx)
   data first/.true./
   save
@@ -50,26 +52,25 @@ subroutine filter(x0,x1)
   if(first) then
      pi=4.0*atan(1.0)
      do i=0,NFFT-1
-        w(i)=(sin(i*pi/NFFT))**2
+        ww=sin(i*pi/NFFT)
+        w(i)=ww*ww/NFFT
      enddo
-     s=0.
-     fac=1.0/NFFT
+     s=0.0
+     f=1.0
+     x0s=0.
+     x1s=0.
      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
+  x(0:NH-1)=x0s                             !Previous 2nd half to new 1st half
+  x(NH:NFFT-1)=x0                           !New 2nd half
+  x0s=x0                                    !Save the new 2nd half
+  x=w*x                                     !Apply window
+  call four2a(x,NFFT,1,-1,0)                !r2c FFT (to frequency domain)
+  cx=f*cx
+  call four2a(cx,NFFT,1,1,-1)               !c2r FFT (back to time domain)
+  x1=x1s + x(0:NH-1)                        !Add previous segment's 2nd half
+  x1s=x(NH:NFFT-1)                          !Save the new 2nd half
 
   return
 end subroutine filter