/*  AutoCorrEstimator.h

    Mark Woolrich, FMRIB Image Analysis Group

    Copyright (C) 1999-2000 University of Oxford  */

/*  CCOPYRIGHT  */

#if !defined(__AutoCorrEstimator_h)
#define __AutoCorrEstimator_h

#include <iostream>
#include <fstream>
#define WANT_STREAM
#define WANT_MATH

#include "miscmaths/miscmaths.h"
#include "newimage/newimageall.h"

using namespace NEWMAT;
using namespace MISCMATHS;
using namespace NEWIMAGE;

namespace FILM {
     
  class AutoCorrEstimator
    {
    public:
      AutoCorrEstimator(const Matrix& pxdata) : 
        sizeTS(pxdata.Nrows()),
	numTS(pxdata.Ncols()),
	xdata(pxdata),
	acEst(pxdata.Nrows(), pxdata.Ncols()),
	vrow(),
	xrow(),
	dummy(),
	realifft(), 
	dm_mn(),
	zeropad(0)
	{ 
	  zeropad = MISCMATHS::nextpow2(pxdata.Nrows());
	  vrow.ReSize(zeropad);
	  xrow.ReSize(zeropad);
	  dummy.ReSize(zeropad);
	  realifft.ReSize(zeropad);
	  countLargeE.ReSize(zeropad);
	  countLargeE = 0;
	  acEst = 0;
	  acEst.Row(1) = 1;
	}

      void calcRaw(int lag = 0);
      void spatiallySmooth(const string& usanfname, const ColumnVector& epivol, int masksize, const string& epifname, int usanthresh, int lag=0);
      void applyConstraints();
      void filter(const ColumnVector& filterFFT);
      Matrix fitAutoRegressiveModel();
      void pava();
      //void preWhiten(VolumeSeries& in, VolumeSeries& ret);
      void preWhiten(ColumnVector& in, ColumnVector& ret, int i, Matrix& dmret, bool highfreqremovalonly=false);
      void setDesignMatrix(const Matrix& dm);
      int establishUsanThresh(const ColumnVector& epivol);

      void getMeanEstimate(ColumnVector& ret);

      Matrix& getEstimates() { return acEst; }
      Matrix& getE() { return E; }
      ColumnVector& getCountLargeE(){ return countLargeE; }

      int getZeroPad() { return zeropad; }
      void tukey(int M);
      void multitaper(int M);
      int pacf(const ColumnVector& x, int minorder, int maxorder, ColumnVector& betas);
      
      volume<float> mask;

    private:
      const int sizeTS;
      const int numTS;
      AutoCorrEstimator();
      const AutoCorrEstimator& operator=(AutoCorrEstimator&);
      AutoCorrEstimator(AutoCorrEstimator&);
      void getSlepians(int M, int sizeTS, Matrix& slepians);

      const Matrix& xdata;
      Matrix acEst;
      Matrix E;
      ColumnVector countLargeE;

      Matrix dminFFTReal;
      Matrix dminFFTImag;

      ColumnVector vrow;
      ColumnVector xrow;     

      ColumnVector dm_fft_real, dm_fft_imag;
      ColumnVector x_fft_real, ac_fft_real;
      ColumnVector x_fft_im, ac_fft_im;
     
      ColumnVector dummy;
      ColumnVector realifft;
      ColumnVector dm_mn;

      int zeropad;
    };
 
}

#endif