program mtm
c
c     driver for MTM spectral analysis/signal reconstruction code
c
c     (c) Michael E. Mann 12/96 
c
c
      parameter (maxlen=50000,zero=0.0,maxsignal=100)
      parameter (nlim=32768)
      character*80 ifmt,name
      real a(maxlen),rcon(maxlen),signal(maxlen)
      real fsignal(maxsignal),confsignal(maxsignal)
      integer irec(maxsignal)
      real dt,anpi,f1,f2,fthresh,f0
      integer nscan,nwin,npi,ioption,nsignals,ithresh
      integer nf0,nbnd,iwhich
      real fconf(6),conf(4)
      real specraw0(nlim),specmed0(nlim),specresh0(nlim),
     $   harmonic0(nlim),specback0(nlim),ftest(nlim)
      real whiteraw0,whiterob0,rhoraw0,rhorob0,tauraw0,taurob0
c
      write (6,*)
      write (6,*) 'MTM SPEC/RECON'
      write (6,*)
      write (6,*) '(c) 1996 Michael E. Mann'
      write (6,*)
c
c     read in data
c
      write (6,*) 'Time series file : '
      read (5,101) name
101   format(80a)
      open (7,file=name,form='formatted')
      rewind(7)
      read (7,*) 
      read (7,*) nscan,dt  
      read (7,101) ifmt
      read(7,ifmt)(a(i),i=1,nscan)
      close(7)
      do i=1,nscan
         demn = demn+a(i)
      end do
      demn = demn/float(nscan)
      do i=nscan+1,maxlen
         a(i)=zero
      end do
c
c     calculate Rayleigh and Nyquist frequencies
c
      fray=1.0/(nscan*dt)
      fny = 0.5/dt
c  
c
c     SET DEFAULTS
c
c     indicate that no reconstructions/plotting possible until
c     spectrum has been calculated...
c
      iflag = 0
c
c     resolution/variance tradeoff
c
      npi = 2
      nwin = 3
      anpi=float(npi)
      fray=1.0/(nscan*dt)
      fny = 0.5/dt
      bndwdth = 2.0*anpi*fray
c
c     spectrum
c
      f1 = 0.0
      f2 = fny
      frange = f2-f1
      inorm = 0
      ispec = 1
      iresh = 1
      ithresh = 3
      fthresh = 95.0
c
c     null hypothesis
c
      inoise = 0
      ismooth = 1
      fsmooth = frange/5.0
      if (fsmooth.lt.bndwdth) fsmooth=bndwdth
      ilog = 1
c
c     signal assumption
c   
      isignal = 0
c
c     reconstruction
c
      irecon = 0
      do i=1,maxsignal
         irec(i)=0
         fsignal(i)=0.0
         confsignal(i)=0.0
      end do
      nsignals = 0
c
c     display
c
      iplotresh = 1
      iplotftest = 1
      iplotsmoo = 1
      iplotraw = 1
      iplotconf = 1
c
c     menu
c
5555  continue
      write (6,*)
      write (6,*) 'MAIN MENU - MTM SPEC/RECON'
      write (6,*)
      write (6,*) '1) variance/resolution tradeoff:'
      write (6,*) '   * resolution = ',npi,' f_R'
      write (6,*) '   * number of tapers = ',nwin
      if (inoise.eq.0) then
        write (6,*) '2) null hypothesis: red noise'
      else
         if (inoise.eq.1) then
            write (6,*) '2) null hypothesis: white noise'
         else
            write (6,*) '2) null hypothesis: locally-white noise'
         endif
      endif
c
      if (isignal.ne.2) then
        if (ismooth.eq.0) then
           write (6,*) '   * raw noise background estimation'
        else
           write (6,*) '   * robust noise background estimation'
           write (6,*) '         & fsmooth = ',fsmooth
           if (inoise.eq.0)  then
              if (ilog.eq.0) then
                write (6,*) '    (min misfit noise background)'
              else
                write (6,*) '    (min log-misfit noise background)'
              endif
           endif
        endif
      endif
        
c
      if (isignal.eq.0) then
        write (6,*) '3) signal assumption: harmonic or narrowband'
      else
         if (isignal.eq.1) then
            write (6,*) '3) signal assumption: narrowband'
         else
            write (6,*) '3) signal assumption: harmonic'
         endif
      endif
c
      write (6,*) '4) spectrum: '
      if (ispec.eq.1) then
         write (6,*) '    * adaptive estimate'
      else
         write (6,*) '    * high-resolution estimate'
      endif
      if (inorm.eq.0) then
         write (6,*)    '    * no normalization'
      else
         if (inorm.eq.1) then
           write (6,*)  '    * normalize by  N '
         else
           write (6,*)  '    * normalize by 1/dt'
         endif
      endif
      if (iresh.eq.0) then
           write (6,*)  '    * unreshaped'
      else
           write (6,*)  '    * reshaped at threshold =',
     $       fthresh,'%'
      endif
      write (6,*)       '    * frequency range: ',f1,f2
c
c
      write (6,*) '5) signals to reconstruct: ',nsignals
c
c
      iplot = iplotresh+iplotsmoo+iplotraw+iplotconf+iplotftest
      if (iplot.eq.0) then
         write (6,*) '6) no display:'
      else
         write (6,*) '6) display:'
         if (iplotftest.eq.1) then
            write (6,*) '   *  F-test spectrum'
            write (6,*) '   *   & 50% 90% 95% 99% sig levels'
         endif
         if (iplotraw.eq.1)  write (6,*) '  * raw spectrum'
         if (iplotresh.eq.1) 
     $      write (6,*) '  * reshaped & harmonic spectra'
         if (iplotsmoo.eq.1) 
     $      write (6,*) '  * median smoothed spectrum'
         if (iplotconf.eq.1) 
     $      write (6,*) '  * 50% 90% 95% 99% sig levels'
      endif
c
      write (6,*) '7) exit'
c
      write (6,*) 'option (continue=0)'
      read (5,*) ioption
c
c     proceed w/ calculations or quit?
c
      if (ioption.eq.0) goto 8888
      if (ioption.eq.7) goto 9999
c
c
c     if any of the attributes of the spectrum are selected to
c     be changed, we nullify the existing set of selected signals
c     since they might not be "signals" under changed assumptions,etc
c
c     also, flag that spectrum, peak detection has to be redone...
c
      if (ioption.lt.5) then
          nsignals=0
          do i=1,maxsignal
             irec(i)=0
          end do
          iflag=0
      endif
c
      if (ioption.eq.1) then
       write (6,*) 'resolution in "p" multiples of Rayleigh frequency'
       read (5,*) npi
       write (6,*) 'number of tapers (suggested = ',2*npi-1,' ) '
       read (5,*) nwin
       anpi=float(npi)
       bndwdth = 2.0*anpi*fray
       if (fsmooth.lt.bndwdth) fsmooth=bndwdth
       goto 5555
      endif
c      
      if (ioption.eq.2) then
        if (isignal.eq.2) then
           write (6,*) 'harmonic signal-only detection requires'
           write (6,*) 'a "locally-white noise" assumption...'
           goto 5555
        endif 
        write (6,*) 'null hypothesis:'
        write (6,*) '(0) red noise'
        write (6,*) '(1) white noise'
        write (6,*) '(2) locally-white noise'
        read (5,*) inoise
        write (6,*) 'background estimation:'
        write (6,*) '(0) raw'
        write (6,*) '(1) robust'
        read (5,*) ismooth
        if (ismooth.eq.1) then
333         write (6,*) 'median smoothing width'
            write (6,*) 'choose: ',bndwdth,' < f <',frange/2.0
            read (5,*) fsmooth
            if ((fsmooth.gt.frange/2.0).or.(fsmooth.lt.bndwdth)) then
              write (6,*) 'your frequency width is out of bounds!'
              goto 333
            endif
            if (inoise.eq.0) then
               write (6,*) 'red noise background fit:'
               write (6,*) '(0) linear'
               write (6,*) '(1) log'
               read (5,*) ilog
            endif
        endif
        goto 5555
      endif

      if (ioption.eq.3) then
        write (6,*) 'signal assumption'
        write (6,*) '(0) harmonic or narrowband'
        write (6,*) '(1) narrowband'
        write (6,*) '(2) harmonic'
        read (5,*) isignal
        if (isignal.eq.2) then
          inoise=2
          iplotsmoo = 0
          iplotconf = 0
        endif
      endif

      if (ioption.eq.4) then
c
         write (6,*) '(0) high-resolution or (1) adaptive estimate'
         read (5,*) ispec
c
         write (6,*) 'normalization:'
         write (6,*) '(0) none (1) N  (2) 1/dt'
         read (5,*) inorm
c
         write (6,*) 'reshaping (yes=1)?'
         read (5,*) iresh
         if (iresh.eq.1) then
          write 
     $      (6,*) 'sig. level for harmonic peak detection/reshaping'
          write 
     $      (6,*) '(0) 50% (1) 90% (2) 95% (3) 99% (4) 99.5% (5) 99.9%'
          read (5,*) ithresh
          ithresh = ithresh + 1  
          if (ithresh.eq.1) fthresh=50.0
          if (ithresh.eq.2) fthresh=90.0
          if (ithresh.eq.3) fthresh=95.0
          if (ithresh.eq.4) fthresh=99.0
          if (ithresh.eq.5) fthresh=99.5
          if (ithresh.eq.6) fthresh=99.9
         endif
c
         write (6,*) 'fmin:'
         read (5,*) f1
         write (6,*) 'fmax:'
         read (5,*) f2
         frange = f2-f1
         fsmooth = frange/5.0
         if (fsmooth.lt.bndwdth) fsmooth=bndwdth
         goto 5555
      endif


      if (ioption.eq.6) then
222      write (6,*) 'display options (toggle)'
         write (6,*) '(1) raw spectrum                : ',iplotraw
         write (6,*) '(2) reshaped & harmonic spectra : ',iplotresh
         write (6,*) '(3) median smoothed background  : ',iplotsmoo
         write (6,*) '(4) 50% 90% 95% 99% conf levels : ',iplotconf
         write (6,*) '(5) F-test spectrum             : ',iplotftest
         write (6,*) 'option to change: (0=continue)'
         read (5,*) ichoice
         if (ichoice.eq.0) goto 5555
         if (ichoice.eq.1) iplotraw=1-iplotraw
         if (ichoice.eq.2) iplotresh=1-iplotresh
         if (ichoice.eq.3) iplotsmoo=1-iplotsmoo
         if (ichoice.eq.4) iplotconf=1-iplotconf
         if (ichoice.eq.5) iplotftest=1-iplotftest
         if (isignal.eq.2) then
            iplotsmoo = 0
            iplotconf = 0
         endif
         if (isignal.eq.1) then
            iplotftest = 0
            iplotresh=0
         endif
         goto 222
      endif

      irecon=1
      write (6,*) '4) signals to reconstruct: ',nsignals
      if (nsignals.gt.0) then
c
c       constraint option
c
c
c       signal selection
c
444     write (6,*) 'signal frequency signif% reconstruct'
        do i=1,nsignals
          if (fsignal(i).lt.bndwdth) then
           write (6,778) '(',i,')   TREND  ' ,
     $       int(confsignal(i)+0.5),'%',irec(i)
          else
           write (6,777) '(',i,') ',fsignal(i),
     $       int(confsignal(i)+0.5),'%',irec(i)
          endif
        end do
        write (6,779) '(',nsignals+1,')    ADD NEW frequency'
        write (6,*) 'select (toggle) 0=continue'
        read (5,*) iwhich
        if (iwhich.eq.0) goto 8888
        if (iwhich.eq.nsignals+1) then
           write (6,*) 'frequency?'
           read (5,*) fnew
           nsignals=nsignals+1
           fsignal(nsignals)=fnew
           goto 444
        else
         irec(iwhich)=1-irec(iwhich)
         goto 444
        endif
      else
        write (6,*) 'no signals to reconstruct'
        write (6,*) '(',nsignals+1,') ',
     $    'enter your own frequency?'
        write (6,*) 'select (toggle) 0=continue'
        read (5,*) iwhich
        if (iwhich.eq.0) then
           if (nsignals.eq.0) goto 5555
           if (nsignals.gt.0) goto 8888
        endif
        write (6,*) 'frequency?'
        read (5,*) fnew
        nsignals=nsignals+1
        fsignal(nsignals)=fnew
        goto 444
      endif
      write (6,*) 'constraint option:'
      write (6,*) 
     $ '(0) min misfit (1) min norm (2) min slope (3) min rough'
      read (5,*) icon

777   format (a1,i2,a2,f9.4,i9,a1,i6)
778   format (a1,i2,a11,i9,a1,i6)
779   format (a1,i2,a22)
c
c
c
8888  continue

      if (iflag.eq.0) then
       call spec(a,nwin,npi,dt,nscan,f1,f2,
     $     inoise,ismooth,ilog,fsmooth,
     $     isignal,ispec,inorm,iresh,ithresh,
     $     nf0,df,specraw0,specmed0,specresh0,harmonic0,
     $     specback0,ftest,
     $     conf,fconf,nbnd,
     $     whiteraw0,whiterob0,rhoraw0,rhorob0,tauraw0,taurob0)

       call findpeaks(flow,fhigh,isignal,
     $      nf0,df,nbnd,specraw0,specmed0,specresh0,harmonic0,
     $      specback0,ftest,
     $      conf,fconf,
     $      fsignal,confsignal,nsignals)

        call display(nf0,f1,df,
     $     specraw0,specmed0,specresh0,harmonic0,
     $     specback0,ftest,conf,fconf,
     $     whiteraw0,whiterob0,rhoraw0,rhorob0,
     $     tauraw0,taurob0,
     $     iplotresh,iplotftest,iplotsmoo,
     $     iplotraw,iplotconf)
        iflag=1
        goto 5555
      endif
  
      if (irecon.eq.1) then
        do i=1,maxlen
           signal(i)=0.0
           rcon(i)=0.0
        end do
        do j=1,nsignals
         if (irec(j).eq.1) then
           f0 = fsignal(j)
           call recon(a,nwin,npi,dt,nscan,bndwdth,f0,icon,rcon)
           do i=1,nscan
              signal(i)=signal(i)+rcon(i)
           end do
         endif
        end do
        open (unit=3,file='recon.out',status='unknown')
        do i=1,nscan
            write (3,*) float(i-1)*dt,signal(i),a(i)-demn
        end do
        close (unit=3)
        endif
        goto 5555
c
9999  continue
      write (6,*) 'quitting...'
      stop
      end