c************
c  Description:	Host version of QFSHIFT routine
c		FFT based frequency domain shift and block filter operation
c
c  Author:	JGS   9/20/90
c
c  $Revision: 1.1 $
c  $Modtime:   04/04/95 16:56:54  $
c  
c************
	program RQFSHIFT
c
	implicit none
c
	include 'primitive.inc'
	integer*4 bin,bout
	character*20 label
	character*4 formout
	logical*4 cplex
	real*8    delx,delf,fshiftin,fshift,fcenter,fwidth,fnyquist
c
	record /HEADER/ HIN,HOUT
c
	integer*4 nfft, nffth, nblk, nbins, nbin1, nbin2, idfs, idfc, idfw
	integer*4 ngot, noff, nput
	real*4	scale
c
	call M$CHECKIN

	call M$WARNING('RQFSHIFT has serious signal processing errors')
	call M$WARNING('Recommend using SBTUNER instead (see Explains)')
c
c open and test signal input file #1
	HIN.file_name = M$APICK(1)
	HIN.type_allow = '1000'	
	HIN.format_allow = 'S#,C#'
	call M$INOK (HIN)
c
c get the fftsize and blocking factor and complex mode
	nfft = M$GET_SWITCH('FFT')
	if (nfft.lt.2) nfft = 8192
	nffth = nfft/2
	delx = HIN.xdelta
	delf = 1.d0/(delx*nfft)
	fnyquist = .5/delx
	nblk = max(1,M$LPICK(6))
	if (nblk.gt.nfft) call M$ERROR 
     &	('Cannot block average by more than the FFT size')
	scale = 1./float(nfft*nblk)
	cplex = HIN.format(1:1) .eq. 'C'
c
c setup output file
	call M$PROPAGATE (HIN,HOUT)
	HOUT.file_name = M$APICK(2)
	formout = M$APICK(7)
	if (formout.ne.' ') HOUT.format = M$TEST_FORMAT(formout,'S#,C#')
	HOUT.xstart = HIN.xstart + INT(nblk/2)*HIN.xdelta
	HOUT.xdelta = HIN.xdelta*nblk
	HOUT.size = HIN.size/nblk
	call M$ALLOCATE (HOUT)
c
c get freq shift
	fshiftin = M$DPICK(3)	
	nbins = NINT( fshiftin/delf )
	fshift = nbins*delf
	label = M$UPICK(8)
	call M$DRSLT (label,fshift)
c
c get filter center and width
	fwidth = M$DPICK(4)
	fcenter = M$DPICK(5)
	if (fwidth.le.0.d0) fwidth = fnyquist*2
	nbin1 = max(-nffth, NINT((fcenter-.5*fwidth)/delf) )
	nbin2 = min( nffth, NINT((fcenter+.5*fwidth)/delf) )
c
c setup real time controls
	idfs = M$DWINIT ('FREQ SHIFT',  fshift, -fnyquist, fnyquist, delf)
	idfc = M$DWINIT ('FILT CENTER', fcenter, -fnyquist, fnyquist, delf)
	idfw = M$DWINIT ('FILT WIDTH',  fwidth, delf, fnyquist*2, delf)
c
c set up the data flow parameters
	HIN.buf_type = 'F'  ! for GRABX will reformat if needed
	HOUT.buf_type = 'F' ! for FILAD will reformat if needed
c
c get buffers for the data
	call M$MALLOC (nfft*8, bin)
	call M$MALLOC (nfft*8, bout)
c
c signal start of processing block
	call M$SYNC()
c
c loop through buffers
	do while ( M$DO(nfft,HIN.xfer_len) )
	  HIN.cons_len = HIN.xfer_len/nblk * nblk
c
c grab the data and reformat
	  call M$GRABX (HIN,%val(bin),ngot)
c
c get real time freq shift and filter parameters
	  if (M$DWGET(idfs,fshiftin) .ne. 0) then
	    nbins = NINT( fshiftin/delf )
	    fshift = nbins*delf
	    call M$DRSLT (label,fshift)
	  endif
	  if (M$DWGET(idfw,fwidth) + M$DWGET(idfc,fcenter) .ne. 0) then
	    nbin1 = max(-nffth, NINT((fcenter-.5*fwidth)/delf) )
	    nbin2 = min( nffth, NINT((fcenter+.5*fwidth)/delf) )
	  endif
c
c see if we can take the quick way out
	  if (nbins.eq.0 .and. nblk.eq.1) then
	    if (nbin1.le.-nffth .and. nbin2.ge.nffth) then
	      if (HOUT.format(1:1) .eq. HIN.format(1:1)) then
	        call M$FILAD (HOUT,%val(bin),ngot)
	        goto 10
	      endif
	    endif
	  endif
c
c zero pad (if necessary) the buffer and FFT
	  if (cplex) then
	    call M$VFILL (%val(bin+ngot*8), 0.0, (nfft-ngot)*2)
	    call M$CFFT (%val(bin),nfft,1)
	  else
	    call M$VFILL (%val(bin+ngot*4), 0.0, nfft-ngot)
	    call M$RFFT (%val(bin),nfft,1)
	    ! freq domain hilbert transform approximation
	    call M$VFILL (%val(bin+nfft*4+8), 0.0, nfft-2) 
	    call M$VSMUL (%val(bin),0.5,%val(bin),1)
	  endif
c
c rotate the spectra and scale
	  if (nbins.lt.0) nbins = nfft + nbins
	  noff = nfft-nbins
	  call M$VSMUL (%val(bin+noff*8), scale, %val(bout), nbins*2)
	  call M$VSMUL (%val(bin), scale, %val(bout+nbins*8), noff*2)
c
c apply the block filter
	  if (nbin1.gt.-nffth) then		! low stop
	    if (nbin1.ge.0) then
	      call M$VFILL (%val(bout+nfft*4+4), 0.0, nfft-1)
	      call M$VFILL (%val(bout), 0.0, nbin1*2 )
	    else
	      call M$VFILL (%val(bout+nfft*4+4), 0.0, nfft+nbin1*2-1)
	    endif
	  endif
	  if (nbin2.lt.nffth) then		! high stop
	    if (nbin2.ge.0) then
	      call M$VFILL (%val(bout+nbin2*8), 0.0, nfft-nbin2*2-1)
	    else
	      call M$VFILL (%val(bout), 0.0, nfft+1)
	      call M$VFILL (%val(bout+nfft*8+nbin2*8+8), 
     &                       0.0, -(nbin2+1)*2)
	    endif
	  endif
c
c inverse FFT
	  call M$CFFT (%val(bout),nfft,-1)
c
c block average the data
	  if (nblk.gt.1) then
	    call M$CVBLK (%val(bout), %val(bout), nblk, ngot)	    
	  endif
c
c write to file
	  nput = ngot/nblk
	  if (HOUT.format(1:1).eq.'C') then
	    call M$FILAD (HOUT,%val(bout),nput)
	  else
	    call M$VMOV (%val(bout),2,%val(bout),1,nput)
	    call M$FILAD (HOUT,%val(bout),nput)
	  endif
c
  10	enddo
c
c close files
	call M$CLOSE (HIN)
	call M$CLOSE (HOUT)
	call M$RETURN
c
	end

