diff --git a/lib/refspectrum.f90 b/lib/refspectrum.f90
index 0f234a46f..91de08910 100644
--- a/lib/refspectrum.f90
+++ b/lib/refspectrum.f90
@@ -4,11 +4,51 @@ subroutine refspectrum(id2,brefspec)
 !  id2       i*2        Raw 16-bit integer data, 12000 Hz sample rate
 !  brefspec  logical    True when accumulating a reference spectrum
 
-  integer*2 id2(3456)
-  logical brefspec
+  parameter (NFFT=6912,NH=NFFT/2)
+  integer*2 id2(NH)
+  logical brefspec,brefspec0
+  real x(NFFT)
+  real s(0:NH)
+  complex cx(0:NH)
+  equivalence(x,cx)
+  data nsave/0/,brefspec0/.false./
+  save brefspec0,nsave
 
-!  write(*,3001) id2(1:10),brefspec
-!3001 format(10i5,L8)
+  if(brefspec) then
+     if(.not.brefspec0) then
+        nsave=0
+        s=0.
+        brefspec0=.true.
+     endif
+     
+     x(1:NH)=0.001*id2
+     x(NH+1:)=0.0
+     call four2a(x,NFFT,1,-1,0)                 !r2c FFT
+
+     do i=1,NH
+        s(i)=s(i) + real(cx(i))**2 + aimag(cx(i))**2
+     enddo
+     nsave=nsave+1
+
+     if(mod(nsave,34).eq.0) then                   !About 9.8 sec 
+        df=12000.0/NFFT
+        ia=nint(500.0/df)
+        ib=nint(2500.0/df)
+        call pctile(s(ia),ib-ia+1,50,xmed)
+        db0=db(xmed)
+        nhadd=10
+        open(16,file='refspec.dat',status='unknown')
+        do i=1,NH
+           freq=i*df
+           ia=max(1,i-nhadd)
+           ib=min(NH,i+nhadd)
+           smo=sum(s(ia:ib))/(ib-ia+1)
+           write(16,1000) freq,db(smo)-db0
+1000       format(2f10.3)
+        enddo
+        close(16)
+     endif
+  endif
 
   return
 end subroutine refspectrum