/*  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/volumeseries.h"
#include "miscmaths/volume.h"
#include "newmatap.h"
#include "newmatio.h"
#include "miscmaths/miscmaths.h"

using namespace NEWMAT;
using namespace MISCMATHS;

namespace FILM {
     
  class AutoCorrEstimator
    {
    public:
      AutoCorrEstimator(const VolumeSeries& pxdata) : 
	xdata(pxdata),
	acEst(pxdata.getNumVolumes(), pxdata.getNumSeries(), pxdata.getInfo(), pxdata.getPreThresholdPositions()),
	vrow(),
	xrow(),
	dummy(),
	realifft(), 
	dm_mn(),
	zeropad(0)
	{ 
	  zeropad = MISCMATHS::nextpow2(pxdata.getNumVolumes());
	  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 Volume& epivol, int masksize, const string& epifname, const string& susanpath, int usanthresh, int lag=0);
      void applyConstraints();
      void filter(const ColumnVector& filterFFT);
      void 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 Volume& epivol);

      void getMeanEstimate(ColumnVector& ret);

      VolumeSeries& getEstimates() { return acEst; }
      VolumeSeries& 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);

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

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

      VolumeSeries dminFFTReal;
      VolumeSeries 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