c
      program mtmcohere
c
c     (c) Michael Mann
c
c     MTM code of J. Park adapted to perform significance of spectral
c     peaks relative to "robust red noise" background as described
c     by Mann and Lees (1996) and to calculate MTM spectral coherence
c     of two series with jacknife uncertainties and automatic determination
c     of confidence levels.
c
      integer maxwin,maxdof,nlim,nmed
      parameter (maxwin=8,maxdof=2*maxwin,nlim=32768)
      parameter (nmed=2000)
      parameter (big=1.0e+36,small=1.0e-12,tol=1e-4)
      real rho,rho0
      real FUNC
      real*8 avar,el
      real val(nmed)
      real tm
      real rmean(nlim),specmed(nlim)
      real harmonic(nlim),background(nlim),psd(nlim)
      real rednoise(nlim),rednoise0(nlim)
      real dat(nlim,2)
      real base(nlim)
      real amp(nlim),phase(nlim),amphigh(nlim),amplow(nlim)
      real coh(maxwin)
      real takeep(nlim/2,maxwin,2),taikeep(nlim/2,maxwin,2)
      integer iharm(nlim)
      real ftest(6),fcoh(6),fcoh0(6),fconf(6)
      real adum
      character*1 k0,k1,k2,k3,k4,chan*8,ista*4
      character*80 ifmt,name,title
      complex zed
      common/npiprol/anpi
      common/data/a(50000)
      common/work/b(32800)
      common/work2/fspec(16400)
      COMMON /BLK1/white,specmed,fny,ddf,f0
      COMMON /BLK2/ilog,nfreq
c     program is limited to 8 windows
c     program is limited to 16400 frequency-ftest estimation pairs
c     #pts in padded series must be < 32768 '
      common/pltopt/key,k1,k2,k3,k4,nbin,npad,nst,npts,nmin,nmax,t(3),
     $     fmin,fmax,nf
      common/taperz/ta(16400,8),tai(16400,8),dcf(16400,8),amu(16400,8)
      common/staple/tad(32800),tad1(32800)
      common/stap2/tas(32800,16)
      common/idiot/qcyc,rfrq,ftfrq,fr(3),k0
      common/udhedz/iy,id,ih,im,iss,ktime(5),dt,ista,chan,nscan,commen
      dimension dum(35)
      dimension reim(2),el(10)
      equivalence (zed,reim),(iy,dum)
      pi=3.14159265358979
      radian=180./pi
      tiny = 1e-6
      iflag=1
c
6666  write (6,*) 'options: (quit=0)'
      write (6,*) '(1) MTM spectrum of a time series'
      write (6,*) '(2) MTM time-domain signal reconstruction'
      write (6,*) '(3) MTM spectral coherence of two time series'
      read (5,*) ioption
      ioption = ioption-1
      if (ioption.eq.1) then
         write (6,*) 'sorry--that section is still in progress'
         write (6,*) 'check back soon for updated version!'
         goto 6666
      else
       if (ioption.lt.0) goto 8888
      endif
c
c
c     read in data file/s
c
c     first 3 lines contain
c        1. title
c        2. #samples, sample interval
c        3. fortran format statement
c
101   format(80a)
c
      if (ioption.lt.2) then
c
        nrep = 1
        write (6,*) 'input file name : '
        read (5,101) name
        open(7,file=name,form='formatted')
        rewind(7)
        read(7,101) title
        read(7,*) nscan,dt
        read(7,101) ifmt
        read(7,ifmt)(dat(i,nrep),i=1,nscan)
        close(7)
        open (unit=44,file='data.out',status='unknown')
        do i=1,nscan
           write (44,*) i,dat(i,1)
        end do
        close (unit=44)
c
      else
        nrep = 2
        do irep = 1,nrep
          write (6,*) 'input file name for series #',
     $      irep,':'
          read (5,101) name
          open(7,file=name,form='formatted')
          rewind(7)
          read(7,101) title
          read(7,*) nscan1,dt1
          read(7,101) ifmt
          read(7,ifmt) (dat(i,irep),i=1,nscan1)
          close(7)
          if (irep.eq.1) then
             nscan=nscan1
             dt = dt1
          else
             if (abs(dt1-dt).gt.tol) then
                write (6,*) 'sampling rates not equal'
                write (6,*) 'cannot compute coherences...'
                goto 6666
             else
              if (nscan1.ne.nscan) then
                write (6,*) 'lengths not equal'
                write (6,*) 'series1: ',nscan,' samples'
                write (6,*) 'series2: ',nscan1,' samples'
                write (6,*) 'taking shorter of the two:'
                nscan = min(nscan,nscan1)
                write (6,*) 'using first ',nscan,' data points...'
              endif
             endif
          endif
        end do
      endif
c
      write (6,*) 'read in: '
      write (6,*) nrep,' time series of length ',nscan
      write (6,*) float(int(100*nscan*dt+0.5))/100.0, ' years long...'
      write (6,*) 'sampling rate ',dt
c
c
c     determine MTM parameters
c
c     default settings
c
      npi = 2
      anpi=float(npi)
      nwin=3
      imode = 0
      nmode = 1
      iseg = 0
      iformat = 0
      iref = 1
      jsmoo = 0
      iaverage = 0
      irecon = 0
      imoderecon = 0
      ioutput = 0
      itype = 0
      nbeg=1
      nend = nscan
      n1 = nbeg
      n2 = nend
      npts = n2-n1+1
      neval = nscan
      nmove = nscan
      ninterv = nmove
      fsmooth = 0.0
      nsmooth = 0
      inoise = 1
      ismooth = 1
      ilog = 1
      tm0 = 0.0
c
444   write (6,*) 'present settings'
      write (6,*) '1) use',nwin, ' x ',npi,' pi tapers'
      write (6,*) '2) data interval (',nbeg,' to ',nend,' )'
      if (itype.eq.0) then
        write (6,*) '3) fixed window analysis'
      else
         write (6,*) '3) evolutive analysis'
         write (6,*) '    --',nmove,' time unit moving window'
         write (6,*) '    --',ninterv,' time unit step'
      endif
      if (ioption.gt.1) goto 555
c
      if (inoise.eq.0) then
         write (6,*) '4) white noise assumption'
      else
        if (inoise.eq.1) then
           write (6,*) '4) red noise assumption'
        else
           if (inoise.eq.2) then
              write (6,*) '4) median-estimated "local white" noise'
           endif
        endif
      endif
555   write (6,*) 'item to change (0=continue)'
      read (5,*) icont
      if (icont.eq.1) then
       goto 1111
      else
       if (icont.eq.2) then
         goto 2222
       else
         if (icont.eq.3) then
          goto 3333
         else
            if (icont.eq.4) then
               goto 4444
            else
              goto 9999
            endif
         endif
       endif
      endif
c
c     anpi = desired time/frequency bandwith "p" (ie, "p pi tapers")
c     nwin = number of tapers used (maximum is 2p-1)
c
1111  write (6,*) 'time-frequency bandwidth product "p"'
      read (5,*) npinew
      if (npinew.gt.maxwin) then
         write (6,*) 'sorry--exceeds maximum bandwidth'
      else
         write (6,*) 'number of tapers'
         read (5,*) nwinnew
         if (nwinnew.gt.2*npinew-1) then
           write (6,*) 'sorry--cannot exceed 2*p-1'
         else
           npi = npinew
           anpi = float(npi)
           nwin=nwinnew
         endif
      endif
      goto 444
c
2222  write (6,*) 'beginning time index'
      read (5,*) nbeg
      write (6,*) 'ending time index'
      read (5,*) nend
      nmove = nend-nbeg+1
      neval = nend-nbeg+1
      n1 = nbeg
      n2 = nend
      goto 444
c
3333  write (6,*) 'choose'
      write (6,*) '(0) fixed window'
      write (6,*) '(1) evolutive analysis'
      read (5,*) itype
      if (itype.eq.1) then
         write (6,*) 'window width (in time units)'
         read (5,*) nmove
         write (6,*) 'interval between evaluations (in time units)'
         read (5,*) ninterv
         write (6,*) 'initial time (year,month, etc)'
         read (5,*) tm0
      endif
      goto 444
c
c
4444  write (6,*) 'choose'
      write (6,*) '(0) white noise background assumption'
      write (6,*) '(1) red noise background assumption'
      write (6,*) '(2) median estimated "local white" assumption'
      read (5,*) inoise
      if (inoise.eq.0) then
         rho = 0.0
         rho0 = rho
      else
        fsmooth = 0.0
        nsmooth = 0
c
        demn = 0.0
        do i=n1,n2
         demn=demn+a(i)
        end do
        demn = demn/float(npts)
c
c       determine raw AR(1) coefficient 
c
        var = 0.0
        do i=n1,n2
         a(i)=dat(i,1)
         rmean(i)=a(i)-demn 
         var = var + rmean(i)**2
        end do
        var = var/float(npts)
        sd = sqrt(var)
        write (6,*) 'series mean: ',demn
        write (6,*) 'series standard deviation: ',sd
        c1 = 0.0
        icount = 0
        do i=2,n2
         icount = icount + 1
         c1 = c1 + rmean(i-1)*rmean(i)
        end do
        c1 = c1/float(icount)
        rho = c1/var
        write (6,*) 'rho (raw series): ',rho
        rho0 = rho
        if (inoise.eq.1) then
         write (6,*) 'robust noise background option'
         write (6,*) '(0) no'
         write (6,*) '(1) yes'
         read (5,*) ismooth
        else
          ismooth = 1
        endif
      endif
      goto 444
c
c
c
9999  npts = nmove
c
c     determine chi-square values for confidence level
c     determination of spectrum
c
      open (unit=21,file='/big2/MTM-COHERE/chisquare.dat',
     $ status='old')
         do i=2,maxdof
            idofs = 2*nwin
            if (i.eq.idofs) then
               read (21,*) idum,conf50,conf90,conf95,conf99
            else
               read (21,*) 
            endif
         end do
      close (unit=21)
c
      ick1 = 0
      ick2 = 1
      ick3 = 1
      jper = 3
      bfact = conf95
      fconf(1)=50.0
      fconf(2)=90.0
      fconf(3)=95.0
      fconf(4)=99.0
      fconf(5)=99.5
      fconf(6)=99.9
      if (nwin.eq.1) ick3 = 0
c
c     read in f-test data
c
      if (ick3.eq.1) then
       iper = 3
       open (unit=17,file='/big2/MTM-COHERE/ftest.dat',
     $    status='unknown')
       do i=2,14,2
         if (i.eq.2*nwin-2) then
            read (17,*) idum,ftest(1),ftest(2),ftest(3),ftest(4),
     $            ftest(5),ftest(6)
         else
            read (17,*)
         endif
       end do
       close (unit=17)
      endif
      thresh = ftest(iper)   
c
      if (ioption.eq.2) then
       do i=1,6
          fcoh(i)=0.0
          fcoh0(i)=0.0
       end do
       if (nwin.eq.1) then
          write (6,*) 'warning: '
          write (6,*) 'confidence intervals are indeterminate'
          write (6,*) 'for the case K=1...'
       endif
c
c      read in coherence signficance levels for non-secular
c
       open (unit=17,file='/big2/MTM-COHERE/fcoher.dat',
     $ status='unknown')
       do i=2,14,2
         if (i.eq.2*nwin-2) then
            read (17,*) idum,fcoh(1),fcoh(2),fcoh(3),fcoh(4),
     $            fcoh(5),fcoh(6)
         else
            read (17,*)
         endif
       end do
       close (unit=17)
c
c      read in coherence signficance levels for secular band
c
c      note we are comparing 1 and nwin-2 dof since dft
c      is real at zero frequency, and the mean counts as
c      one dof (which we have removed)
c
       open (unit=17,file='/big2/MTM-COHERE/fcoher-sec.dat',
     $    status='unknown')
       do i=0,6
         if (i.eq.nwin-2) then
            read (17,*) idum,fcoh0(1),fcoh0(2),fcoh0(3),fcoh0(4),
     $            fcoh0(5),fcoh0(6)
         else
            read (17,*)
         endif
       end do
       close (unit=17)
       bfact = fcoh(jper)
       bfact0 = fcoh0(jper)
c
      endif
c
      if(iflag.eq.0) go to 60
      key=0
      k2='n'
      k1='y'
      nbin=2
      inorm = 0
      k4='y'
      k3='y'
      qcyc=0.
      k0='n'
c
c  find the rayleigh freq
      rfrq=1./(npts*dt)
      if(npts.gt.32768) type *,'too large - cut it down'
      npad=npts-1
      ij=0
 1001 npad=npad/2
      ij=ij+1
      if(npad.ne.0) go to 1001
      npad=2**ij
      write (6,*) ,npad
      fmin=0.
      f0 = fmin
      fmax=.5/dt
      nst=0
      hstart=nst*dt
      thrs=npts*dt
   60 ddf=1.0/(npad*dt)
      eps=.1*ddf
      nn1=(fmin+eps)/ddf
      fmin=nn1*ddf
      f0 = fmin
      nn2=((fmax+eps)/ddf) 
      fmax=nn2*ddf
      nf=nn2-nn1+1
      nfreq = nf
      write (6,*) ' number of freqs = ',nf
      nmin=nn1*2+1
      nmax=nn2*2+2
      t(1)=fmin
      t(2)=ddf
      t(3)=ddf
      ftfrq=ddf
      lwt=2
      
   50 write (6,*) ' 1) freq spacing(can be decreased by padding): ',ddf
      write (6,*) ' 2) min and max freq for analysis: ',fmin,fmax
      if (ioption.eq.2) then
           write (6,*) ' 3) filter spectral coherence at ',fconf(jper),
     $          ' significance level'
      else
      write (6,*) ' 3) spectrum type:'
      if (ick1.eq.1) then
        write (6,*) '  * high resolution MTM spectrum'
      else
        write (6,*) '  * adaptive MTM spectrum'
      endif
        write (6,*) '  * filter spectrum at ',fconf(jper),
     $        ' signficance level'
      if (ick3.eq.1) then
        write (6,*) '  * harmonic peak detection/reshaping'
        write (6,*) '  * ',fconf(iper),
     $      '% sig. level for harm. peak detection'
      else
        write (6,*) '  * no harmonic peak detection/reshaping'
      endif
      write (6,*) ' 4) units - amp spectrum per unit time '
      endif
666   write (6,*) 'which number do you want to change (0=continue) '
      read (5,*) ijk
      if(ijk.eq.0) goto 9876
c
c     set iflag if user selects 1 or 2 
c     iflag=1 causes eigentapers to be reinterpolated
c
      go to (1,2,3,4), ijk
c
    1 write (6,*) 'current padded length',npad
      write (6,*) 'pad by how many powers of two? '
      read (5,*) ij
      ip=2**(iabs(ij))
      nd=npad*ip
      if(ij.lt.0) nd=npad/ip
      if(nd.gt.32768) then
        type *,'too many - max padded length=32768'
        go to 1
      endif
      npad=nd
      write (6,*) ' padded fft length = ',npad
      go to 60
c
    4 write (6,*) '(0) no normalization'
      write (6,*) '(1) normalization = dft*npts'
      write (6,*) '(2) normalization = dft/dt'
      read (5,*) inorm
      go to 50
c
    2 write (6,*) 'freq units are in cycles per time unit'
      write (6,*) 'enter fmin and fmax '
      read (5,*) fmin,fmax
      go to 60
c
    3 if (ioption.eq.2) then
c
        write (6,*) 'significance level (%) for filtering'
        write (6,*) '(0) 50%'
        write (6,*) '(1) 90%'
        write (6,*) '(2) 95%'
        write (6,*) '(3) 99%'
        write (6,*) '(4) 99.5%'
        write (6,*) '(5) 99.9%'
        read (5,*) jper
        jper = jper+1
        bfact = fcoh(jper)
        bfact0 = fcoh0(jper)
c
      else
c
      write (6,*) 'choose:'
      write (6,*) '(0) high-resolution MTM spectrum'
      write (6,*) '(1) adaptive MTM spectrum'
      read (5,*) ick2
      if ((ick2.gt.1).or.(ick2.lt.0)) then
         write (6,*) 'illegal entry'
         goto 4
      endif
      ick1 =  1-ick2
        write (6,*) 'significance level (%) for filtering'
        write (6,*) '(0) 50%'
        write (6,*) '(1) 90%'
        write (6,*) '(2) 95%'
        write (6,*) '(3) 99%'
        write (6,*) '(4) 99.5%'
        write (6,*) '(5) 99.9%'
        read (5,*) jper
        jper = jper+1
        bfact = conf50
        if (jper.eq.2) bfact = conf90
        if (jper.eq.3) bfact = conf95
        if (jper.gt.3) bfact = conf99
      write (6,*) 'F-test for harmonic lines (0=no,1=yes)'
      read (5,*) ick3
      if (ick3.eq.1) then
       if (nwin.eq.1) then
           write (6,*) 'sorry--F test is indeterminate for K=1'
           goto 60
       else
        write (6,*) 'significance level (%) for harmonic peak detection'
        write (6,*) '(0) 50%'
        write (6,*) '(1) 90%'
        write (6,*) '(2) 95%'
        write (6,*) '(3) 99%'
        write (6,*) '(4) 99.5%'
        write (6,*) '(5) 99.9%'
        read (5,*) iper
        iper = iper+1
        thresh = ftest(iper)
       endif
      endif
c
      endif
      goto 50
c
c     we're happy with the settings. Continue...
c
9876  ddf = 1.0/(dt*npad)
      fny = 0.5/dt
      if(nf.ge.16400) then
        write (6,*) 'too many DFT points'
        go to 444
      endif
c
      afact = conf50
      if (iper.eq.2) afact = conf90
      if (iper.eq.3) afact = conf95
      if (iper.gt.3) afact = conf99
c
      bndwdth = 2.0*anpi*rfrq  
      halfwdth = bndwdth/2.0
      nbnd = bndwdth/ddf
      timefreq = 2.0*anpi
      write (6,*) 'Rayleigh frequency: ',rfrq
      write (6,*) 'time-frequency bandwidth: ',timefreq,' N'
      write (6,*) 'spectral bandwidth: ',bndwdth
      write (6,*) 'Nyquist frequency: ',fny
      fsugg = fny/5.0
      if (fsugg.lt.bndwdth) fsugg=bndwdth
      if ((inoise.gt.0).and.(ioption.eq.0)) then
828      write (6,*) 'smoothing width (cycles/yr):'
         write (6,*) 'suggested value = ',fsugg
         read (5,*) fsmooth
         if ((fsmooth.gt.fny).or.(fsmooth.lt.bndwdth)) then
          write (6,*) 'spectral bdwdth < width < Nyquist freq!'
           write (6,*) 'your frequency width is out of bounds'
           goto 828
         endif
         nsmooth = fsmooth/ddf
         if (inoise.eq.1) then
          write (6,*)
          write (6,*) 'fit red noise to smoothed'
          write (6,*) '(0) spectrum'
          write (6,*) '(1) log spectrum'
          read (5,*) ilog
         endif
      endif
c
c  if we have not changed anything, we need not recalculate the windowed
c  spectra (e.g. if we wish to re-interpolate spectra)
c  get slepian windows - unless we can skip this step
c
      if(iflag.eq.1) then
        call taper2(npts,nwin,el)
      endif
c
c     refresh iflag
c
      iflag=0
c
c
c    normalization: mult by dt if we wish absolute spectral estimate
c    e.g. analysis of time-limited signal
c    divide by npts if we wish amplitude spectrum per unit time
c    or divide by sqrt(npts)  if we are adhering to conventions of 
c    inverse theory paper (since tapers self-dotted = n)
c
      if(inorm.eq.2) then  !  amp spec (time limited)
        anrm=1/dt
      elseif(inorm.eq.1) then  ! amp spect per unit time
        anrm=float(npts)
      else
        anrm=1.
      endif
c
c
c     begin the analysis....
c
c     create output files
c
      open (unit=31,file='coher-amp.out',status='unknown')
      open (unit=32,file='coher-phase.out',status='unknown')
      open (unit=33,file='coher-sig.out',status='unknown')
      open (unit=21,file='spec-harm.out',status='unknown')
      open (unit=22,file='spec-robst.out',status='unknown')
      open (unit=23,file='spec-raw.out',status='unknown')
      open (unit=24,file='spec-med.out',status='unknown')
c
      open (unit=29,file='spec-inf-evol.out',status='unknown')
c
c     loop in time (if non-evolutive, reduces to a single window)
c
      write (6,*) 'nmove ',nmove
      write (6,*) 'ninterv ',ninterv
      niter = (neval-nmove)/ninterv + 1
      write (6,*) 'beginining 1st of ',niter, ' iteration/s...'
      amax = -99999
c
c     begin moving time window loop
c
      do it = 1,niter
c
c
      n1 = 1+(it-1)*ninterv
      n2 = n1+npts-1
      ncent = (n1+n2-1)/2
      tm = ncent*dt
c
      write (6,*) 'interval: ',n1,' to ',n2,' time: ',tm
c
c     demean series in local window -- loop over two series
c     if spectral coherence option
c
      do k=1,nrep
c
      demn=0.
      do i=n1,n2
         a(i)=dat(i,k)
         demn=demn+a(i)
      end do
      demn=demn/(npts)
c
      do iwin=1,nwin
        do i=n1,n2
          b(i-n1+1)=(a(i)-demn)*tas(i+1-n1,iwin) 
        end do
        j=npts+1
        do i=j,npad
          b(i)=0.
        end do
        call realft(b,npad/2,1)
        b(npad-1)=b(2)
        b(npad)=0.
        b(2)=0.
        sum=0.
        j=0
        do i=nmin,nmax,2
          j=j+1
          ta(j,iwin)=b(i)/anrm
          tai(j,iwin)=b(i+1)/anrm
          takeep(j,iwin,k)=ta(j,iwin)
          taikeep(j,iwin,k)=tai(j,iwin)
          sumi=(b(i)*b(i)+b(i+1)*b(i+1))
          sum=sum+sumi
        end do
        avamp=sqrt(sum/nf)/anrm
      end do
c
c
      end do
      if (ioption.eq.2) then
c
c     calculate the spectral coherence and jacknife uncertainties
c
       do i=1,nf-1
        abr=0.
        abi=0.
        aa=0.
        bb=0.
        cmean = 0.
        do k=1,nwin
          coh(k)=0.
          aa=aa+takeep(i,k,1)*takeep(i,k,1)
     $      +taikeep(i,k,1)*taikeep(i,k,1)
          bb=bb+takeep(i,k,2)*takeep(i,k,2)
     $      +taikeep(i,k,2)*taikeep(i,k,2)
          abr=abr+takeep(i,k,1)*takeep(i,k,2)
     $      +taikeep(i,k,1)*taikeep(i,k,2)
          abi=abi+takeep(i,k,1)*taikeep(i,k,2)
     $      -taikeep(i,k,1)*takeep(i,k,2)
          abr2=0.
          abi2=0.
          aa2=0.
          bb2=0.
          do j=1,nwin
             if (j.ne.k) then
               aa2=aa2+takeep(i,j,1)*takeep(i,j,1)
     $           +taikeep(i,j,1)*taikeep(i,j,1)
               bb2=bb2+takeep(i,j,2)*takeep(i,j,2)
     $           +taikeep(i,j,2)*taikeep(i,j,2)
               abr2=abr2+takeep(i,j,1)*takeep(i,j,2)
     $           +taikeep(i,j,1)*taikeep(i,j,2)
               abi2=abi2+takeep(i,j,1)*taikeep(i,j,2)
     $           -taikeep(i,j,1)*takeep(i,j,2)
             endif
          enddo
          coh(k)=(abr2*abr2+abi2*abi2)/(aa2*bb2)
          cmean = cmean + coh(k)
        end do
        cmean = cmean/float(nwin)
        var = 0.
        do k=1,nwin
         var = var + (coh(k)-cmean)**2
        end do
        unc = sqrt(var/float(nwin))
        amp(i)=(abr*abr+abi*abi)/(aa*bb)
        phase(i)=180.*atan2(abi,abr)/pi
        amphigh(i)=amp(i)+unc
        amplow(i)=amp(i)-unc
        ff = (i-1)*ddf+fmin
        if (itype.eq.0) then
           write (31,854) ff,amp(i),amplow(i),amphigh(i)
           write (32,855) ff,phase(i)
           if (ff.ge.halfwdth) then
            write (33,856) ff,amp(i),fcoh(1),fcoh(2),fcoh(3),fcoh(4)
           else
            write (33,856) ff,amp(i),fcoh0(1),fcoh0(2),fcoh0(3),
     $                        fcoh0(4)
           endif
        else
           write (31,954) tm,ff,amp(i)
           write (32,955) tm,ff,phase(i)
           if (ff.ge.halfwdth) then
             if (amp(i).gt.bfact) then
              write (33,955) tm,ff,amp(i),phase(i)
             else
              write (33,955) tm,ff,0.0,-200.0
             endif
           else
             if (amp(i).gt.bfact0) then
              write (33,955) tm,ff,amp(i),phase(i)
             else
              write (33,955) tm,ff,0.0,-200.0
             endif
           endif
        endif
       end do
       goto 777
      endif
c
c     high-resolution mtm spectrum
c
      if(ick1.eq.1) then
        call hires(ta,tai,el,nwin,nf-1,amu)
c
c       normalized psd, disallow small negative values
c
        do i=1,nf-1
          amu(i,1)=abs(amu(i,1))**2
        end do
      endif
c
c     adaptively weighted mtm spectrum 
c
      if(ick2.eq.1) then
        avar=0.d0
        do i=n1,n2
          avar=avar+(a(i)-demn)*(a(i)-demn)
        end do
c
c       avar is a factor that scales the bias factor 
c       in adaptive weighting
c       scheme. it will depend on the transform normalization thusly:
c
        if (inorm.eq.1) then ! amp spect per unit time
          avar=avar/(npts*npts)
        elseif(inorm.eq.2) then  ! absolute amp spect
          avar=avar*dt*dt
        endif
        call adwait(ta,tai,dcf,el,nwin,nf-1,amu,amu(1,2),avar)
c
c       normalized psd, disallow small negative values
c
        do i=1,nf-1
          amu(i,1)=abs(amu(i,1))
        end do
      endif
c
c
      do i=1,nf-1
         psd(i)=amu(i,1)
      end do
c
      if (it.eq.1) then
c
c     determine white noise level
c     and  determine median smoothed spectrum
c      (if moving window, based on first window)
c
      white0 = 0.0
      do i=1,nf-1
         white0 = white0+psd(i)
      end do
      white0 = white0/float(nf-1)
      rho = rho0
      white = white0
      endif
c
c
c     estimate spectrum amplitude
c
      do i=1,nf-1
         psd(i)=psd(i)+white0*small
      end do
c
      if (it.eq.1) then
       do i=1,nf-1
         red = white0*(1.0-rho0**2)/
     $      (1.0-2.0*rho0*cos(arg)+rho0**2)
         if ((red.gt.amax).and.(it.eq.1)) amax = red
       end do
c
       white = 0.0
       do j=1,nf-1
         if1 = j-nsmooth/2
         if2 = j+nsmooth/2
         if (if1.lt.1) if1=1
         if (if2.gt.nf-1) if2=nf-1
         nblk = if2-if1+1
         do i=1,nblk
            val(i)=psd(if1+i-1)
         end do
c
c        sort spectrum in this block
c         
         kmid = (nblk+1)/2
         do kk1=1,nblk
            do kk2=kk1+1,nblk
               if (val(kk2).lt.val(kk1)) then
                  adum = val(kk2)
                  val(kk2)=val(kk1)
                  val(kk1)=adum
               endif
            end do
         end do
         specmed(j)=val(kmid)
         white = white + specmed(j)
       end do
       white = white/float(nf-1)
c
c      now determine best fit red noise spectrum of the form
c
c         rednoise = white*(1.0-rho**2)/
c     $      (1.0-2.0*rho*cos(arg)+rho**2)
c
c      to the median-smoother.
c
c      we take "white" as estimated above from median smoothed
c      spectrum, and do a global search of the interval [0,1)
c      to find the optimal rho as determined by minimum MSE
c
       if (ismooth.eq.1) then
        amin = big
        do rho1=0.0,0.999,0.001
         amiss = FUNC(rho1)
         if (amiss.lt.amin) then
            rho=rho1
            amin = amiss
         endif
        end do
       endif
c
c     amax is just a scaling factor to insure a standard
c     magnitude scale of the output spectrum
c
      amax = amax/100.0
c
      tau = 1e+16
      tau0 = 1e+16
      if (rho.gt.0.0) tau = -dt/log(rho)
      if (rho0.gt.0.0) tau0 = -dt/log(rho0)
      write (29,657) n1,n2,white,rho,tau
      write (6,*) n1,n2,white,rho,tau
c
      endif
c
      do i=1,nf-1
         ff = (i-1)*ddf+fmin
         freq_norm = ff/fny
         arg = freq_norm*pi
         rednoise0(i) = white0*(1.0-rho0**2)/
     $      (1.0-2.0*rho0*cos(arg)+rho0**2)
         rednoise(i) = white*(1.0-rho**2)/
     $      (1.0-2.0*rho*cos(arg)+rho**2)
         if (inoise.ne.2) then
             base(i) = rednoise(i)
         else
             base(i) = specmed(i)
         endif
c
         if (itype.eq.0) then
c
         write (22,654) ff,psd(i)/amax,base(i)/amax,
     $         base(i)*conf90/amax,
     $         base(i)*conf95/amax,base(i)*conf99/amax
         write (23,654) ff,psd(i)/amax,rednoise0(i)/amax,
     $         rednoise0(i)*conf90/amax,
     $         rednoise0(i)*conf95/amax,
     $         rednoise0(i)*conf99/amax
         write (24,655) ff,psd(i)/amax,specmed(i)/amax
c
         else
c
         if (psd(i).gt.bfact*base(i)) then
           write (22,754) tm+tm0,ff,log(psd(i)/amax)
         else
           write (22,754) tm+tm0,ff,0.0
         endif
         write (23,754) tm+tm0,ff,log(psd(i)/amax)
c
         endif
      end do
c
      if (ick3.eq.1) then
c
        call regre(ta,tai,nf-1,nwin,amu)
c
c       do a poor mans reshaping -detect harmonic peaks and then
c       interpolate the continuous spectrum across the effected
c       bandwidth  -only reshape if harmonic peak is greater than
c       the significance level in terms of overall power that was
c       indicated for the F-test harmonic detection procedure
c
c
        do i=1,nf-1
           background(i)=psd(i)
           harmonic(i)=background(i)
           iharm(i)=0
           if ((amu(i,3).gt.thresh).
     $         and.(psd(i).gt.afact*base(i))) iharm(i)=1
        end do
        do i=1,nf-1
           if (iharm(i).eq.1) then
c
c           determine frequency points at boarder of the bandwidth
c
            ipre = i-nbnd/2-1
            iaft = i+nbnd/2+1
            if (ipre.lt.1) ipre=1
            if (iaft.gt.nf-1) iaft=nf-1 
c
            ave = 0.0
            do j=ipre+1,iaft-1
              background(j)=0.5*(psd(ipre)+psd(iaft))
              harmonic(j)=psd(i)
            end do
           endif
        end do
        do i=1,nf-1
           ff = float(i-1)*ddf+fmin
           if (itype.eq.0) then
           write (21,656) ff,harmonic(i)/amax,background(i)/amax,
     $        amu(i,3),ftest(1),ftest(2),ftest(3),ftest(4)
           else
            alarge = max(afact,bfact)
            if (harmonic(i).gt.alarge*base(i)) then
              write (21,756) tm,ff,log(harmonic(i)/amax),
     $           amu(i,3)
            else
              write (21,756) tm,ff,0.0,amu(i,3)
            endif
           endif
        end do
      endif
c
c     finished with evolutive loop
c
777   end do
c
      open (unit=27,file='spec-inf.out',status='unknown')
      if (itype.eq.1) write (27,*) 'final window:'
      write (27,*) 'smoothing:'
      write (27,*) '# points: ',nsmooth
      write (27,*) '   bndwdth: ',fsmooth
      write (27,*) '   log(yes=1) ',ilog
      write (27,*) 'normalization option: ',inorm
      write (27,*) 'anrm: ',anrm
      write (27,*) 'normalizing factor for output: ',amax
      write (27,*) 'raw determination:'
      write (27,*) '    white noise variance: ',white0
      write (27,*) '    rho: ',rho0
      write (27,*) '    tau: ',tau0
      write (27,*) 'robust determination:'
      write (27,*) '    white noise variance: ',white
      write (27,*) '    rho: ',rho
      write (27,*) '    tau: ',tau
c
      close (unit=21)
      close (unit=22)
      close (unit=23)
      close (unit=24)
      close (unit=27)
      close (unit=29)
      close (unit=31)
      close (unit=32)
      close (unit=33)
      goto 6666
c
8888  continue
c
654   format (f7.4,1x,5f12.6)
655   format (f7.4,1x,2f12.6)
656   format (f7.4,1x,3f12.6,4f6.2)
657   format (2i5,6f9.3)
c
754   format (2f9.4,1x,1f12.6)
756   format (2f9.4,1x,2f12.6)
c
854   format (f7.4,1x,3f12.6)
855   format (f7.4,1x,1f12.6)
856   format (f7.4,1x,5f6.3)
c
954   format (2f9.4,1x,1f12.6)
955   format (2f9.4,1x,2f12.6)
956   format (2f9.4,1x,5f6.3)
      stop
      end
c
      real function FUNC(rho)
c
      parameter (nlim=32768)
      COMMON /BLK1/white,specmed(nlim),fny,ddf,f0
      COMMON /BLK2/ilog,nfreq
      real rho
      real pie,dff,freq_norm,ff,arg,small,rednoise,val1,val2
      pie = 3.1415926
      small = 1e-12
      dff = 0.0
      do j=1,nfreq-1
         ff = (j-1)*ddf+f0
         freq_norm = ff/fny
         arg = freq_norm*pie
         rednoise = white*(1.0-rho**2)/
     $      (1.0-2.0*rho*cos(arg)+rho**2)
         if (ilog.eq.0) then
           dff = dff + (specmed(j)-rednoise)**2
         else
           val1 = abs(specmed(j))
           val2 = abs(rednoise)
           if (val1.lt.small) val1=small
           if (val2.lt.small) val2=small
           dff = dff + 
     $     (log(val1)-log(val2))**2
         endif
      end do
      FUNC = dff
      return
      end