c
c $Revision: 1.2 $
c $Modtime:   9/30/96 10:58:32  $
cC **********************************************************************
C *                    H E A D E R  B L O C K                          *
C **********************************************************************
C
C ######################################################################
C #                                                                    #
C #                                                                    #
C #   ##############                                                   #
C #   #            #                                                   #
C #   #  NEWRATE   #                                                   #
C #   #            #                                                   #
C #   ##############                                                   #
C #                                                                    #
C #   NAME:  NEWRATE                                                   #
C #                                                                    #
C #   PURPOSE:                                                         #
C #      Resample an input file at a new sampling rate                 #
C #                                                                    #
C #   DESCRIPTION:                                                     #
C #      Get command line parameters                                   #
C #      Set-up interpolation tables                                   #
C #      Perform interpolation calculations                            #
C #      Write buffer of data                                          #
C #      Post headers, close files and return to MIDAS                 #
C #                                                                    #
C #   CALLING SEQUENCE:                                                #
C #      NEWRATE <input> <output> <new sampling rate> <sync file>      #
C #      <filter file>                                                 #
C #                                                                    #
C #   COMMON:                                                          #
C #      PRIMITIVE.INC (which contains the necessary MIDAS commons)    #
C #                                                                    #
C #   INPUT FILES:                                                     #
C #      MIDAS signal input file                                       #
C #      MIDAS sync file file (optional for variable sample rate)      #
C #      MIDAS type 2000 filter file (optional)                        #
C #                                                                    #
C #   OUTPUT FILES:                                                    #
C #      MIDAS file of resampled signal                                #
C #                                                                    #
C #   SUBROUTINES CALLED:                                              #
C #      MIDAS system utility routines                                 #
C #                                                                    #
C #   RESTRICTIONS:                                                    #
C #      All type 1000 real data types are permitted                   #
C #      LONG_INTEGER(I*4), INTEGER(I*2), or REAL.                     #
C #                                                                    #
C #   MODIFICATIONS:                                                   #
C #      REV 0.0 - CRAY MIDAS release                    DSL   01/92   #
C #      REV 1.0 - X-MIDAS release                       DSL           #
C #      REV 2.0 Add filter interpolation                DSL   07/93   #
C #      REV 2.1 Corrected problem with buffer sizes     DSL   09/93   #
C #      REV 2.2 Changed next_clock algorithm for                      #
C #              improved accuracy with long files.      DSL   10/93   #
C #      REV 2.3 Changed looping for faster execution                  #
C #              Also added implicit none and declared                 #
C #              more variables as real*8 for accuracy   DSL   11/93   #
C #                                                                    #
C #   DR XI0143  Design enhancement for the inclusion                  #
C #              of this primitive into the SYS option   SRC   1/24/95 #
C #              1) Added M$FORCE1000                                  #
C #              2) Created NEWRATE.INC for block data   SRC   8/30/95 #
C #                                                                    #
C #      REV 2.4 Support complex filter                  DSL    1/96   #
C #      REV 2.5 Added comments                                        #
C #              Fixed potential bug preventing exit     dsl    9/96   #
C #      REV 2.6 Added /TUNE=dval switch                               #
C #              Added tune widget                                     #
C #              Changed BLOCKDATA subprogram to simple                #
C #              data include to avoid bug in Promula                  #
C #              Fortran.                                dsl    7/97   #
C #      REV 2.7 Moved widget init statements ahead of                 #
C #              m$sync.                                 dsl   1/27/99 #
C #                                                                    #
C ######################################################################
c Principal of operation
c
c Newrate steps through an input with steps that need not be an integral
c multiple of input elements. The steps size is determined by the desired
c output sample rate. This step can be modified by an input sync file or
c by changing the widget that controls the sample rate. The variable
c aindex determines the fractional index of the input file for the desired
c output point. The output sample as actually aligned at aindex+(ntaps/2)-1.
c The fractional part of the index is used to select a filter from a bank
c that provides a sequence of fine time steps. when generating the output.
c The default filter uses 28 taps to accurately generate fractional sample
c time displacements for frequencies up to 0.475 of the sample rate. Other
c filters can be generated by the RESAMPF primitive for a variety of
c purposes such as low-pass filtering (necessary for rate reduction),
c or real to complex data conversion.
c
c Using a sync file.
c
c The sync file is intended to support baud synchronous resampling or
c alignment of signals from different moving platforms. The sync file
c presents displacements from the nominal output sample positions. The
c displacements are supplied as a fraction of the output clock interval. The
c displacements can cover the range of -1/2 to +1/2 of a clock cycle, relative
c to an extrapolation from the last sync point at the specified output rate.
c Values outside of this range are adjusted by an integral number of output
c cycles. The xdelta for the sync file must be larger than both the xdelta
c of both the input file and the nominal rate of the output file. In a typical
c application, it may be a thousand times larger. The sync file may be
c obtained from a complex baud clock by using POLAR to get phase in degrees,
c and then multiplying by the samples_per_baud/360.
c
c Program design
c
c The program produces samples at a constant output rate. Sync data is
c used to calculate an adjusted sample rate between sync points so that
c the clock phase matches the sync phase at each sync point. Output data
c values are effectively interpolations of the input data samples. The
c number of interpolated positions is determined by the filter size, and 
c the size of the filter table. The default filter uses 2048 interpolation
c positions. Internal buffer management routines retain needed input data.
c Output is accumulated to reduce the overhead write operations. The actual
c filtering loop uses the maximin possible number of iterations, subjet to
c input data available, room in the output buffer, and the possible need
c to set a new sample rate at the next phase point.
c
c Care has been taken to minimize cumulative errors in abscissa values when
c processing large files, and to reduce the number of operations in the
c inner processing loop. This has resulted in some increase in program
c complexity. Comments in the data declaration section indicate the
c intended function of some of the control variables.
C
C **********************************************************************
C *               D E C L A R A T I O N   B L O C K                    *
C **********************************************************************
C
C
      program newrate
      implicit none
      INCLUDE 'primitive.inc'
      integer blksize,blksize2,filtwrds,filtsize
      PARAMETER (blksize=131072,filtwrds=262144)
      record/header/ HCB1,hcb2,hcb5,hcb4
      CHARACTER*80 NAME1,name2,name4,name5,dummy_str
C
      INTEGER*4 blkend2
C      
      INTEGER*4 TYPE_CODE1, FORMAT_CODE1, X_UNITS1
      REAL*8    X1_START,x1_delta, SIZE1, START1,PTIME,DX,tstart
      real*8 time_now,aindex,astart,t2,sigfirst,siglast
      logical syncfile,cxfilt,tuning
C
      REAL*4    sigbuf(blksize),r4buf2(blksize),r4buf5(filtwrds)
      complex csigbuf(blksize/2),c8buf5(filtwrds/2),c8buf2(blksize/2)
      equivalence (sigbuf,csigbuf),(r4buf5,c8buf5),(r4buf2,c8buf2)
      real*8 index_now,astep
      real*8 p_start,pp_delta,new_rate,x2_delta,x2_start,count
      integer pbufsize,pbufsize_max
      parameter (pbufsize_max=16384)
      real*4 phasebuf(pbufsize_max+1),rsum,qsum,halffilt,round
      complex csum
      integer*4 nbpe,ntaps,nread,kintr,ngot,ninterp,idr,nstep,ncleft
      integer idt
      integer*4 nhave,i,idum,nrep,irep,ngot5,ifirst,ioff,itoff
      real*8 tstep, tune, tune_dph
      complex*16 tune_cx,tune_twiddle
      common/syninfo/hcb4,p_start,pp_delta,x2_delta,pbufsize,phasebuf 
      common /filter/ntaps,ninterp,halffilt(28*17)
C
C
C
C **********************************************************************
C *               I N I T I A L I Z A T I O N   B L O C K              *
C **********************************************************************
C
c      call lib$init_timer()
      CALL m$CHECKIN
C
C     Pick-up parameters from command line
C
      new_rate = m$DPICK(3)
      name4=m$apick(4)
      name5=m$apick(5)

c     open input file
      hcb1.file_name=m$apick(1)
      hcb1.type_allow='1000'
      hcb1.format_allow='S#,C#'
      call m$inok(hcb1)
      x1_start=hcb1.xstart
      x1_delta=hcb1.xdelta
      size1=hcb1.data_size/hcb1.bpe
      hcb1.buf_type='F'
      nbpe=hcb1.spa
      blksize2=blksize/nbpe
c      write(6,*)'nbpe,blksize2 ',nbpe,blksize2
      x2_delta=1./new_rate
c      write(6,*)'x2_delta ',x2_delta
c
c     open sync file (if named)
c
      if(name4(1:1).eq.' ')then
       syncfile=.false.
      else
       syncfile=.true.
       hcb4.file_name=name4
       hcb4.type_allow='1000'
       hcb4.format_allow='S#'
       call m$inok(hcb4)
       if(hcb4.xdelta .le. x2_delta .or. hcb4.xdelta .le. x1_delta)
     1  call m$error('Syncfile xdelta too small. Try thinning it.') 
c     rate widget only available if not using sync file
      end if
c     The rate widget is ignored if a synd file is used. It is always
c     created so the the TUNE widget will be winget 2.
      idr=M$DWINIT('RATE',new_rate,.9*new_rate,1.1*new_rate,
     1     .001*new_rate)
c
c     open filter file or use default filter
c
      if(name5(1:1).eq.' ')then
       cxfilt=.false.
       filtsize=filtwrds
       call default_filter(r4buf5)
c       write(6,9000)(r4buf5(i),i=1,ntaps)
c9000   format(' first row of default filter'/(5f12.6))
      else
       hcb5.file_name=m$apick(5)
       hcb5.type_allow='2000'
       hcb5.format_allow='SF,CF'
       call m$inok(hcb5)
       ntaps=hcb5.subsize
       nread=hcb5.size
       cxfilt=hcb5.format.eq.'CF'
       if(cxfilt)then
        filtsize=filtwrds/2
        blksize2=blksize/2
       else
        filtsize=filtwrds
       end if
      end if
c      write(6,*)'cxfilt,filtsize ',cxfilt,filtsize
c       write(6,*)'ntaps,nread ',ntaps,nread
       x2_start=x1_start+(real(ntaps)/2.-1.)*x1_delta
c       write(6,*)'x2_start ',x2_start
      index_now=1
      kintr=0
c
c     open output file, now that we have all required info.
c
      call m$propagate(hcb1,hcb2)
      hcb2.file_name=m$apick(2)
      hcb2.xstart=x2_start
      hcb2.xdelta=x2_delta
      hcb2.size=size1*x1_delta/x2_delta
      hcb2.format=hcb1.format
      if(cxfilt)hcb2.format='CF'
c      write(6,*)'output format ',hcb2.format
c
c     The TUNE widget is ignored if the output is not complex
c
      idt=M$DWINIT('TUNE',tune,1.d0,-1.d0,1d0)
      call m$allocate(hcb2)
      hcb2.buf_type='F'
      siglast = -1e6
      indexb = 0
      ngot = -1
      count = 0
c force 1000 type file handling. This may not be a good idea, but it is
c a MCCB convention.
       call M$FORCE1000 (HCB1,HCB2)

c
c     Get /tune switch.
c
      if(m$get_uswitch('TUNE',dummy_str).gt.0)then
       tune=m$get_dswitch('TUNE')
       if(hcb2.format.ne.'CF' .and. tune.ne.0d0)then
        call m$warning('TUNE switch ignored for real output.')
        tune = 0d0
        tuning = .FALSE.
       endif
       tune_dph = dmod(tune * hcb2.xdelta*dtwopi,dtwopi)
       tune_cx = dcmplx(dcos(tune_dph),dsin(tune_dph))
       if(tune.ne.0)then
        tuning = .TRUE.
       else
        tuning = .FALSE.
       endif
      else
       tune = 0d0
       tuning = .FALSE.
      endif
c	write(6,*)'tune_cx',tune_cx,tuning
      tune_twiddle = dcmplx(1d0,0d0)
c
      call m$sync
c
c     read in filter before processing any data
c
      if(name5(1:1).ne.' ')then
       call m$grab(hcb5,r4buf5,1.0d0,nread,ngot5)
c       if(ngot5.ne.ninterp)call m$error('Error reading filter file')
c       write(6,*)'ntaps,ninterp ',ntaps,ninterp,ngot5
       call fixfilter(r4buf5,ntaps,nread,ninterp,filtsize,cxfilt)
c       if(cxfilt)write(6,'(8f10.6)')(c8buf5(i),i=1,2*ntaps)
       round = 0.5
      else
       round = 0.0
      endif
      if(.not. syncfile)then
c      astep is the number of input samples to step for each output sample
c      this is NOT an integer.
       astep = x2_delta/x1_delta
       aindex = 1
      else
       pbufsize=pbufsize_max
c
c      call for initialization of the sync data
c      returns: t2= the time of the first output samples, tstep= the
c      actual-nominal time delta between output samples, nstep=then number
c      of output samples for which tstep is constant. Count= the number of
c      ateps that we have used.
c
       call first_clock(x2_start,0.0,t2,tstep,nstep,count)
       ncleft=nstep
c       write(6,*)'t2,tstep,nstep,count ',t2,tstep,nstep,count
       aindex=(t2+count*x2_delta+p_start-x2_start)/x1_delta+1
       astep=(tstep+x2_delta)/x1_delta
c       kintr=int(ninterp*(aindex-index_now)+round)
      end if
C
C
C **********************************************************************
C *                 P R O C E S S I N G   B L O C K                    *
C **********************************************************************
C
c      write(6,*)'start of process loop '
c      call lib$show_timer()
10    if(Mc.break)goto 991
      if(index_now.gt.siglast)then
c      We need more data from input file
       if(ngot.eq.blksize2)then
c       last read was fully satisfied so we probably can read more
c       but we need to move old data to start of input buffer
        ifirst=index_now-sigfirst+1
        nhave=max(ngot-ifirst+1,0)
c        write(6,*)'ifirst,index_now,nhave ',ifirst,index_now,nhave
        if(nbpe.eq.2)then
         do i=ifirst,ngot
          csigbuf(i-ifirst+1)=csigbuf(i)
         end do
        else
         do i=ifirst,ngot
          sigbuf(i-ifirst+1)=sigbuf(i)
         end do
        end if
       else
c       we get here before first grab and at end of file        
        if(ngot.ge.0)goto 991       !at end of file so exit
c       At start of file. Nothing to save.
        nhave=0
	ngot=0
       end if
c      Read new data to follow old data.
       call m$grab(hcb1,sigbuf(nhave*nbpe+1),index_now+nhave,
     1  blksize2-nhave,ngot)
       ngot=ngot+nhave
c       write(6,*)'nhave,ngot,index_now ',nhave,ngot,index_now
c       write(6,'(5f12.6)')(sigbuf(i),i=ngot-29,ngot)
       if(ngot.lt.ntaps)goto 991
c      file index values for data in input buffer.
       sigfirst=index_now
       siglast=index_now+ngot-ntaps
c       write(6,*)'siglast ',siglast
c      test parameter change
       if(.not. syncfile)then
        idum=m$dwget(idr,new_rate)
        if(idum > 0.5)then
         x2_delta=1.d0/new_rate
         astep=x2_delta/x1_delta
        end if
       end if
       if(hcb2.format.eq.'Cf')then
        idum=m$dwget(idt,tune)
        if(idum > 0.5)then
c        Do not turn off tuning even if tune goes to 0, as this would
c        cause a phase discontuinity.
         tuning = .true.        
         tune_dph = dmod(tune * hcb2.xdelta*dtwopi,dtwopi)
         tune_cx = dcmplx(dcos(tune_dph),dsin(tune_dph))
        end if
       end if
      end if
c
c     Almost ready to filter data. Do a few more initializations.
c     We can now iterate the interpolation process until the earliest
c     of (1) the number of points left in the output buffer, (2) the
c     number of points for which we have input data, or (3) the number
c     of steps left for the current syncfile interval.
c
      nrep = min(blksize2-indexb,int((siglast-aindex+.999)/astep)+1)
      if(syncfile)nrep=min(nrep,ncleft)
c     write(6,*)'nrep,siglast,aindex,astep ',nrep,siglast,aindex,astep
c      call lib$show_timer()
c
      astart=aindex
c
c     This is the actual filter loop
c
      do irep = 0,nrep-1
       aindex=astart+irep*astep
       index_now=dint(aindex)
c      calculate index of the filter to be used for inperpolation.
       kintr = ninterp*(aindex - index_now) + round
       itoff = kintr*ntaps
c       if(kintr.lt.0 .or. kintr.ge.ninterp)
c     1  write(6,*)'kintr,indexb ',kintr,indexb
c       write(6,*)'index_now,ioff,kintr ',index_now,ioff,kintr
c
c      calculate which point in the input buffer is the the first to be
c      filtered
c
       ioff=index_now-sigfirst
c      increment pointer to output buffer
       indexb=indexb+1
       if(cxfilt)then
        if(nbpe.eq.2)then
c        complex input file      
         csum=(0,0)
         do i=1,ntaps
          csum=csum+csigbuf(i+ioff)*c8buf5(i+itoff)
         end do
         if(tuning)then
          c8buf2(indexb)=csum*cmplx(tune_twiddle)
          tune_twiddle=tune_twiddle*tune_cx
         else
          c8buf2(indexb)=csum
         endif
        else
         rsum=0
         qsum=0
         do i=1,ntaps
          rsum=rsum+sigbuf(i+ioff)*r4buf5(2*(i+itoff)-1)
          qsum=qsum+sigbuf(i+ioff)*r4buf5(2*(i+itoff))
         end do
         if(tuning)then
          c8buf2(indexb)=cmplx(rsum,qsum)*cmplx(tune_twiddle)
          tune_twiddle=tune_twiddle*tune_cx
         else
          c8buf2(indexb)=cmplx(rsum,qsum)
         endif
        end if
       else
        if(nbpe.eq.2)then
c        complex input file      
         rsum=0
         qsum=0
         do i=1,ntaps
          rsum=rsum+sigbuf(2*(i+ioff)-1)*r4buf5(i+itoff)
          qsum=qsum+sigbuf(2*(i+ioff))*r4buf5(i+itoff)
         end do
         if(tuning)then
          c8buf2(indexb)=cmplx(rsum,qsum)*cmplx(tune_twiddle)
          tune_twiddle=tune_twiddle*tune_cx
         else
          c8buf2(indexb)=cmplx(rsum,qsum)
         endif
        else
         rsum=0
         do i=1,ntaps
          rsum=rsum+sigbuf(i+ioff)*r4buf5(i+itoff)
         end do
         r4buf2(indexb)=rsum
        end if
       end if
      end do
c      write(6,*)'after irep loop '
c      call lib$show_timer()
c
c     We could output filtered points now, but we may wait to avoid writing
c     small blocks.
c
      if(indexb.ge.blksize2*0.9)then
c       write(6,*)'writing ',indexb, r4buf2(indexb) 
       call m$filad(hcb2,r4buf2,indexb)
       indexb=0
      end if
c
c     Some preparation for the next loop
c
c     recalculate aindex for improved numeric accuracy
      aindex=astart+nrep*astep
      count=count+nrep
      if(syncfile)then
       ncleft=ncleft-nrep
       if(ncleft.lt.1)then
        t2=t2+nstep*tstep
        call next_clock(t2,t2,tstep,nstep,count)
        ncleft=nstep
c        write(6,*)'t2,tstep,Nstep,count ',t2,tstep,nstep,count
        aindex=(t2+count*x2_delta+p_start-x2_start)/x1_delta+1
        astep=(tstep+x2_delta)/x1_delta
       end if
c      else
c       aindex=1+count*astep
      end if
      index_now=dint(aindex)
c      write(6,*)'aindex,index_now ',aindex,index_now
      goto 10
C
C **********************************************************************
C *                  W R A P - U P    B L O C K                        *
C **********************************************************************
C
991   continue
C     Close all files, post headers, and pass control back to MIDAS 
C
      if(indexb.gt.0)THEN
       call m$filad(hcb2,r4buf2,indexb)
      END IF
      call m$close (hcb1)
      CALL m$CLOSE (HCB2)
      call m$close (hcb4)
      call m$close (hcb5)
      CALL m$return
C
      END
c
      subroutine first_clock(time_start,tce,t_now,tstep,nstep,count)
c
c     Sync file HCB and sync period (pri) passed by /syninfo/ common
c     block. Pbufsize may be set to limit sync input for improved operation
c     in a piped macro.
c     Time_start is the earliest time that can be processed from the input
C     file. Abscissa registration relative to the sync file is assumed.
C     Time_last is the time offset of the last symbol. It is set to
C     0  by first_clock.
C     Tce is the offset between the sync time and the time for the 
C     selected sample from the input file. It is implemented by
C     adding tce to the start of the sync file. P_start is, therefore,
C     set equal to the abscissa start of the sync file + tce.
C     Time_now is the time offset of the current sample.
C     Tstep is the offset increment per symbol. The instantaneous 
C     symbol period is tstep+pri.
C     Nstep is the number of periods to be processed with the current 
C     value of tstep.
C     Count is both an input and a returned value. It is ininitialized 
C     by first clock and used as input by next_clock. The sync phase
C     file will be read at: time_last + pri*count + p_start.
C     The next sample to be taken from the input file is at
C     t_now + pri*count + p_start.
c
      implicit none
      include 'primitive.inc'
c
      real*8 time_last,t_last,t_now,time_now,index,bufstart,x_start
      integer pbufsize,bindex,pbufsize_max
      parameter (pbufsize_max=16384)
      real phasebuf(pbufsize_max+1)
      integer*4  ngot,npts,nstep,maxint,nfix
      real*8 tce,pstep_nominal,phase1,time1,phase2,pdif,phasint,tleft
      real*8 p_start,pp_delta,time_start,p_index,count,pri,ioff,tstep
      record /HEADER/ hcbp
      common/syninfo/hcbp,p_start,pp_delta,pri,pbufsize,phasebuf
      data bufstart,x_start/2*0.0d0/
      data pstep_nominal,npts/2*0/,maxint/2147483647/
      logical first

c     Set pbufsize within legal limits
      pbufsize=min(max(pbufsize,2),pbufsize_max)
c*    read one buffer at start
c     write(6,*)' pbufsize ',pbufsize
      call M$GRAB(hcbp,phasebuf(1),1.d0,pbufsize,ngot)
c        write(6,*)'Initial phase grab',pbufsize,ngot
      npts=ngot
      if(ngot.eq.0)call M$error(' Empty phase file')
      p_index=1
      p_start=HCBP.xstart+tce
      bufstart=p_start
      pp_delta=HCBP.xdelta
      pstep_nominal=pp_delta/pri
c      write(6,*)'pri,pstep_nominal ',pri,pstep_nominal
      x_start=time_start
c
c     The initial value of count is an estimate for the purpose of
c     reading the sync file. It may be corrected before exiting.
c
      count=nint((x_start - p_start)/pri-0.5)
      first=.true.
      t_last=0
      go to 10

      entry next_clock(time_last,t_now,tstep,nstep,count)

      first=.false.
      t_last=time_last
10    continue
c      write(6,'(6h count,f5.0,7h t_last,e14.6)')count,t_last
      time_now=t_last+pri*count+p_start
  
c*    test to see if time_now is within the phase buffer
      bindex=int((time_now-bufstart)/pp_delta+1)
c      write(6,9001)time_now,bindex,npts
9001  format(' time_now,bindex,npts ',e14.6,2i5)
c*    test if we are in leadin extrapolation region
      if(bindex.lt.1)then
c     adjust time_now so that it is aligned with symbol clock
       t_now=-amod(phasebuf(1),1.0)*pri
       ioff=nint((t_last-t_now)/pri)
       tstep=0
       t_now=t_now-ioff*pri
       nstep=1-count-nint((t_now)/pri+0.5)
c       write(6,*)'t_now,ioff ,time_now,nstep',t_now,ioff,time_now,nstep
       goto 99
      end if
c*    test if we need more data
      if(bindex.ge.npts)then
       if(ngot.lt.pbufsize)then
c       we are past end of phase file so extrapolate
        tstep=0
        nstep=maxint
        t_now=-phasebuf(ngot)*pri+bufstart-p_start+pp_delta*(ngot-1)
        t_now=t_now+nint((t_last-t_now)/pri)*pri
        goto 99
       endif
       p_index=nint((time_now-p_start)/pp_delta+0.5)
       bufstart=p_start+(p_index-1)*pp_delta
       call M$GRAB(hcbp,phasebuf(1),p_index,pbufsize,ngot)
c       write(6,*)' new phase',p_index,ngot
       npts=ngot
       bindex=int((time_now-bufstart)/pp_delta)+1
      end if

c     we are after first point and before last point so interpolate
      phase1=phasebuf(bindex)
      time1=bufstart+(bindex-1)*pp_delta
      phase2=phasebuf(bindex+1)
      pdif=phase2-phase1
      pdif=pdif+nint(pstep_nominal-pdif)
c      write(6,*)'phase1,phase2,pdif ',phase1,phase2,pdif
      phasint=phase1+pdif*(time_now-time1)/pp_delta
c     from the interpolated phase adjust time_next to the nearest symbol
      t_now=t_last-(phasint-nint(phasint))*pri
c      write(6,*)'t_now,phasint ',t_now,phasint
      tleft=time1+pp_delta-time_now
c      write(6,*)'tleft ',tleft
      if(tleft.lt.0) tleft = tleft + pp_delta
      tstep=pp_delta/pdif-pri
      nstep=int(tleft/pri+1)
C     Test to account for file start problem. Adjust in pri steps at start.
99    if(first)then
       if(t_now+count*pri+p_start .lt. x_start)then
        nfix=int((x_start-(t_now+count*pri+p_start))/pri)+1
       else
        nfix=-int((t_now+count*pri+p_start-x_start)/pri)
       end if
       nstep=nstep-nfix
       count=count+nfix
c       write(6,*)'nfix ',nfix
      end if
c      write(6,9002)count,nstep,t_now,tstep,time_now
9002  format(' count,nstep,t_now,tstep ',f7.0,i11,2e14.6,e18.10)
      return
      end
c
      subroutine default_filter (fbuf)
      real fbuf(262144)
      common /filter/ntaps,ninterp,halffilt(28*17)
      include 'newrate.inc'
      ninterp=2048
c     copy 28*17 halffilt into 28*33, spaced 64
      do i=1,17
       do j=1,28
        fbuf(j+(i-1)*64*28)  =  halffilt(j+(i-1)*28)
        fbuf(29-j+(33-i)*64*28)=halffilt(j+(i-1)*28)
       end do
      end do
c     now use linear interpolation to craeate 28*2048 filters
c     there will be one extra filter row at end which will not be used
      do i=1,32
       do j=1,28
        do k=1,63
         fbuf(j+(k+64*(i-1))*28)=((64-k)*fbuf(j+(i-1)*64*28)+
     1       k*fbuf(j+i*64*28))/64
        end do
       end do
      end do
cc     print out filter for test
c      write(7,1000)((fbuf(i+j*28),i=1,28),j=0,2048,32)
c1000  format(1x,4e16.7)
      return
      end
      subroutine fixfilter(buf,ntaps,nread,ninterp,ibufsize,cxfilt)
      real buf(ibufsize)
      logical cxfilt
c     This subroutine takes an input filter set consisting of nread
c     filters exactly spanning the delay range from -0.5 to 0.5
c     and interpolates to find filters for as many delays as will
c     fit into the available memory.

      nspread=(ibufsize-ntaps)/(ntaps*(nread-1))
      if(nspread.ge.2)then

c      spread out data in buffer
       m=ntaps*nspread
       if(cxfilt)then
        do i=nread,2,-1
         do j=1,ntaps
          buf(2*(j+m*(i-1))-1)=buf(2*(j+ntaps*(i-1))-1)
          buf(2*(j+m*(i-1)))=buf(2*(j+ntaps*(i-1)))
         end do
        end do
       else
        do i=nread,2,-1
         do j=1,ntaps
          buf(j+m*(i-1))=buf(j+ntaps*(i-1))
         end do
        end do
       end if
c      Now interpolate between the initial delay values
       do j=1,nread-1
        jm1=j-1
        do k=2,nspread
         w1=float(k-1)/nspread
         w2=float(nspread-k+1)/nspread
         if(cxfilt)then
          do i=1,ntaps
           buf(2*(i+m*jm1+(k-1)*ntaps)-1)=w2*buf(2*(i+m*jm1)-1)+
     1       w1*buf(2*(i+m*j)-1)
           buf(2*(i+m*jm1+(k-1)*ntaps))=w2*buf(2*(i+m*jm1))+
     1       w1*buf(2*(i+m*j))
          end do
         else
          do i=1,ntaps
           buf(i+m*jm1+(k-1)*ntaps)=w2*buf(i+m*jm1)+w1*buf(i+m*j)
          end do
         end if
        end do
       end do
      end if
1      ninterp=(nread-1)*nspread
      return
      end
