Initial revision

AutoCorrEstimator.bak2

+#include <iostream>
+#include <strstream>
+// Written by Mark Woolrich 6/1/99
+
+#include "AutoCorrEstimator.h"
+#include "miscmaths.h"
+#include "Log.h"
+#include "Volume.h"
+#include "histogram.h"
+
+using namespace NEWMAT;
+using namespace UTILS;
+namespace TACO {
+
+  void AutoCorrEstimator::setDesignMatrix(const Matrix& dm) {
+    Tracer tr("AutoCorrEstimator::setDesignMatrix");
+
+    int sizeTS = xdata.getNumVolumes();
+    int numPars = dm.Ncols();
+    
+    dminFFTReal.ReSize(zeropad, numPars);
+    dminFFTImag.ReSize(zeropad, numPars);
+      
+    ColumnVector dmrow;
+    dmrow.ReSize(zeropad);
+
+    ColumnVector dm_fft_real, dm_fft_imag;
+    ColumnVector dummy(zeropad);
+    ColumnVector realifft(zeropad);
+
+    // FFT design matrix
+    for(int k = 1; k <= numPars; k++)
+      {
+	dummy = 0;
+	dmrow = 0;
+	mn(k) = MISCMATHS::mean(dm.Column(k));
+	dmrow.Rows(1,sizeTS) = dm.Column(k) - mn(k);
+
+	FFT(dmrow, dummy, dm_fft_real, dm_fft_imag);
+	
+	dminFFTImag.Column(k) = dm_fft_imag;
+	dminFFTReal.Column(k) = dm_fft_real;
+      } 
+  }
+
+  void AutoCorrEstimator::preWhiten(ColumnVector& in, ColumnVector& ret, int i, Matrix& dmret) {
+      Tracer tr("AutoCorrEstimator::preWhiten");
+
+      int sizeTS = xdata.getNumVolumes();
+      int numPars = dminFFTReal.getNumSeries();
+
+      ret.ReSize(sizeTS);
+      dmret.ReSize(sizeTS, numPars);
+
+      // FFT auto corr estimate
+      dummy = 0;
+      vrow = 0;
+      vrow.Rows(1,sizeTS/2) = acEst.getSeries(i).Rows(1,sizeTS/2);
+      vrow.Rows(zeropad - sizeTS/2 + 2, zeropad) = acEst.getSeries(i).Rows(2, sizeTS/2).Reverse();
+
+      FFT(vrow, dummy, ac_fft_real, ac_fft_im);      
+
+      // FFT x data
+      dummy = 0;
+      xrow = 0;
+      xrow.Rows(1,sizeTS) = in;
+     
+      FFT(xrow, dummy, x_fft_real, x_fft_im);
+
+      // inverse auto corr to give prewhitening filter:
+      // no DC component so set first value to 0;
+      ac_fft_real(1) = 0.0;	 
+      for(int j = 2; j <= zeropad; j++)
+	{
+	  ac_fft_real(j) = 1.0/ac_fft_real(j);
+	}
+
+      // filter design matrix
+      for(int k = 1; k <= numPars; k++)
+	{
+	  dm_fft_real = dminFFTReal.getSeries(k);
+	  dm_fft_imag = dminFFTImag.getSeries(k);
+	  FFTI(SP(ac_fft_real, dm_fft_real), SP(ac_fft_real, dm_fft_imag), realifft, dummy);
+
+	  // place result into ret:
+	  dmret.Column(k) = realifft.Rows(1,sizeTS);
+	  float std = pow(MISCMATHS::var(dmret.Column(k)),0.5);
+	  dmret.Column(k) = (dmret.Column(k)/std) + mn(k);
+	}
+
+      // Do filtering and inverse FFT:
+      FFTI(SP(ac_fft_real, x_fft_real), SP(ac_fft_real, x_fft_im), realifft, dummy);
+      
+      // place result into ret:
+      ret = realifft.Rows(1,sizeTS);     
+
+    }
+
+  void AutoCorrEstimator::preWhiten(VolumeSeries& in, VolumeSeries& ret)
+    {
+      Tracer tr("AutoCorrEstimator::preWhiten");
+
+      cerr << "Prewhitening... ";
+
+      int sizeTS = xdata.getNumVolumes();
+      int numTS = xdata.getNumSeries();
+
+      ret.ReSize(sizeTS, numTS);
+
+      // make sure p_vrow is cyclic (even function)
+      ColumnVector vrow, xrow;
+      vrow.ReSize(zeropad);
+      xrow.ReSize(zeropad);
+
+      ColumnVector x_fft_real, ac_fft_real;
+      ColumnVector x_fft_im, ac_fft_im;
+      ColumnVector dummy(zeropad);    
+      ColumnVector realifft(zeropad);
+      int co = 1;
+
+      for(int i = 1; i <= numTS; i++)
+	{
+	  // FFT auto corr estimate
+	  dummy = 0;
+	  vrow = 0;
+	  vrow.Rows(1,sizeTS/2) = acEst.getSeries(i).Rows(1,sizeTS/2);
+	  vrow.Rows(zeropad - sizeTS/2 + 2, zeropad) = acEst.getSeries(i).Rows(2, sizeTS/2).Reverse();
+	  FFT(vrow, dummy, ac_fft_real, ac_fft_im);      
+	  
+	  // FFT x data
+	  dummy = 0;
+	  xrow = 0;
+	  xrow.Rows(1,sizeTS/2) = in.getSeries(i).Rows(1,sizeTS/2);
+	  xrow.Rows(zeropad - sizeTS/2 + 2, zeropad) = in.getSeries(i).Rows(2, sizeTS/2).Reverse();
+	  FFT(xrow, dummy, x_fft_real, x_fft_im);
+
+	  // inverse auto corr to give prewhitening filter:
+	  // no DC component so set first value to 0;
+	  ac_fft_real(1) = 0.0;	 
+	  for(int j = 2; j <= zeropad; j++)
+	    {
+	      ac_fft_real(j) = 1.0/ac_fft_real(j);
+	    }
+
+	  // Do filtering and inverse FFT:
+	  FFTI(SP(ac_fft_real, x_fft_real), SP(ac_fft_real, x_fft_im), realifft, dummy);
+	
+	  // place result into ret:
+	  ret.Column(i) = realifft.Rows(1,sizeTS);
+
+	  if(co > 100)
+	    {
+	      co = 1;
+	      cerr << (float)i/(float)numTS << ",";
+	    }
+	  else
+	    co++;
+	}
+
+      cerr << " Completed" << endl;
+    }
+
+  void AutoCorrEstimator::fitAutoRegressiveModel()
+    {
+      Tracer trace("AutoCorrEstimator::fitAutoRegressiveModel");
+      
+      cerr << "Fitting autoregressive model..." << endl;
+      
+      const int maxorder = 10;
+      const int minorder = 1;
+
+      int sizeTS = xdata.getNumVolumes();
+      int numTS = xdata.getNumSeries();
+      
+      // setup temp variables
+      ColumnVector x(sizeTS);
+      ColumnVector order(numTS);
+      ColumnVector betas(maxorder);
+      acEst.ReSize(sizeTS, numTS);
+      acEst = 0;
+      int co = 1;
+       
+      for(int i = 1; i <= numTS; i++)
+	{
+	  x = xdata.getSeries(i).AsColumn();
+	  order(i) = SIGPROC::Pacf(x, minorder, maxorder, betas);
+	  if(order(i) != -1)
+	    {
+	      // Calculate auto corr:
+	      ColumnVector Krow(sizeTS);
+	      Krow = 0;
+	      Krow(sizeTS) = 1;
+	      Krow.Rows(sizeTS-order(i), sizeTS-1) = -betas.Reverse();
+	      
+	      Matrix Kinv(sizeTS, sizeTS);
+	      Kinv = 0;
+	      for(int j = 1; j <= sizeTS; j++)
+		{
+		  Kinv.SubMatrix(j,j,1,j) = Krow.Rows(sizeTS-j+1,sizeTS).t();
+		}
+	      
+	      // Kinv now becomes V:
+	      Kinv = (Kinv.t()*Kinv).i();
+	      acEst.SubMatrix(1,sizeTS/2+1,i,i) = (Kinv.SubMatrix(sizeTS/2, sizeTS/2, sizeTS/2, sizeTS)/Kinv.MaximumAbsoluteValue()).AsColumn();
+	       
+	      if(co > 200)
+		{
+		   co = 1;
+		   cerr << (float)i/(float)numTS << ",";
+		 }
+	      else
+		co++;
+	    }
+	}
+      
+      Log::getInstance().out("order", order);
+      cerr << " Completed" << endl; 
+    }
+  
+  int AutoCorrEstimator::establishUsanThresh(const Volume& epivol)
+    {
+      int usanthresh = 100;      
+      int num = epivol.getVolumeSize();
+
+      Histogram hist(epivol, num/200);
+      hist.generate();
+      float mode = hist.mode();
+      cerr << "mode = " << mode << endl;
+ 
+      float sum = 0.0;
+      int count = 0;
+
+      // Work out standard deviation from mode for values greater than mode:
+      for(int i = 1; i <= num; i++) {
+	if(epivol(i) > mode) {
+	  sum = sum + (epivol(i) - mode)*(epivol(i) - mode);
+	  count++;
+	}
+      }
+
+      int sig = (int)pow(sum/num, 0.5);
+      cerr << "sig = " << sig << endl;
+
+      usanthresh = sig/3;
+
+      return usanthresh;
+    } 
+
+  void AutoCorrEstimator::spatiallySmooth(const string& usanfname, const Volume& epivol, int masksize, const string& epifname, const string& susanpath) {
+    Tracer trace("AutoCorrEstimator::spatiallySmooth");
+    
+    Log& logger = Log::getInstance();
+
+    // Establish epi thresh to use:
+    int usanthresh = establishUsanThresh(epivol);
+     
+    // Setup external call to susan program:
+    char callsusanstr[1000];
+    ostrstream osc3(callsusanstr,1000);
+    string preSmoothVol = "preSmoothVol";
+    string postSmoothVol = "postSmoothVol";
+
+    osc3 << susanpath << " "
+	 << logger.getDir() << "/" << preSmoothVol << " 1 " 
+	 << logger.getDir() << "/" << postSmoothVol << " " 
+	 << masksize << " 3D 0 1 " << usanfname << " " << usanthresh << " " 
+	 << logger.getDir() << "/" << "usanSize" << '\0';
+    
+    // Loop through first third of volumes
+    // assume the rest are zero
+    int factor = 10000;
+
+    // Setup volume for reading and writing volumes:
+    Volume vol(acEst.getNumSeries());
+    vol.setDims(xdata.getDims());
+    vol.setPreThresholdPositions(xdata.getPreThresholdPositions());
+   
+    acEst.setDims(xdata.getDims());
+    acEst.setPreThresholdPositions(xdata.getPreThresholdPositions());
+    int i = 2;
+
+    acEst.unthresholdSeries();
+    
+    cerr << "Spatially smoothing auto corr estimates" << endl;
+
+    cerr << callsusanstr << endl;
+
+    for(; i < 20; i++)
+      {
+	// output unsmoothed estimates:
+	vol = acEst.getVolume(i).AsColumn()*factor;
+	vol.writeAsInt(logger.getDir() + "/" + preSmoothVol);
+
+	// call susan:
+	system(callsusanstr);
+
+	// read in smoothed volume:
+	vol.read(logger.getDir() + "/" + postSmoothVol);
+
+	acEst.getVolume(i) = vol.AsRow()/factor;
+	cerr << ".";
+      }
+    cerr << endl;
+
+    // Clear unwanted written files
+    char rmfilesstr[1000];
+    ostrstream osc(rmfilesstr,1000);
+    osc << "rm -rf " 
+	<< logger.getDir() + "/" + postSmoothVol + "* "
+	<< logger.getDir() + "/" + preSmoothVol + "* "
+	<< logger.getDir() + "/" + epifname + "* "
+	<< logger.getDir() + "/usanSize*" << '\0';
+    
+    cerr << rmfilesstr << endl;
+
+    system(rmfilesstr);
+
+    cerr << "Completed" << endl;
+
+    acEst.thresholdSeries();
+  }
+
+  void AutoCorrEstimator::calcRaw() { 
+    
+    cerr << "Calculating raw AutoCorrs...";      
+
+    Matrix fts;
+    AutoCorr(xdata.t(), acEst, fts, zeropad);
+    acEst = acEst.t();
+
+    cerr << " Completed" << endl;  
+  }
+  
+  void AutoCorrEstimator::filter(const ColumnVector& filterFFT) {
+
+    Tracer tr("AutoCorrEstimator::filter");
+
+    cerr << "Combining temporal filtering effects with AutoCorr estimates... ";
+
+    // This function adjusts the autocorrelations as if the
+    // xdata has been filtered by the passed in filterFFT
+    // DOES NOT filter the xdata itself
+    ColumnVector vrow;
+
+    // make sure p_vrow is cyclic (even function)
+    vrow.ReSize(zeropad);
+
+    ColumnVector fft_real;
+    ColumnVector fft_im;
+    ColumnVector dummy(zeropad);    
+    ColumnVector realifft(zeropad);
+    int sizeTS = xdata.getNumVolumes();
+
+    for(int i = 1; i <= xdata.getNumSeries(); i++)
+      {
+	dummy = 0;
+	vrow = 0;
+	vrow.Rows(1,sizeTS/2) = acEst.getSeries(i).Rows(1,sizeTS/2);
+	vrow.Rows(zeropad - sizeTS/2 + 2, zeropad) = acEst.getSeries(i).Rows(2, sizeTS/2).Reverse();
+      
+	FFT(vrow, dummy, fft_real, fft_im);
+
+	FFTI(SP(fft_real, filterFFT), dummy, realifft, dummy);
+	
+	// place result into acEst:
+	acEst.Column(i) = realifft.Rows(1,sizeTS)/realifft(1);
+      }
+
+    cerr << " Completed" << endl;
+    
+  }
+  
+  void AutoCorrEstimator::pava() {
+    
+    Tracer tr("AutoCorrEstimator::pava");
+    
+    cerr << "Using New PAVA on AutoCorr estimates... ";
+
+    for(int i = 1; i <= xdata.getNumSeries(); i++) {
+	int sizeTS = xdata.getNumVolumes();
+	int stopat = (int)sizeTS/2;
+
+	// 5% point of distribution of autocorr about zero
+	const float th = (-1/sizeTS)+(2/sqrt(sizeTS));
+
+	ColumnVector values = acEst.Column(i);
+	ColumnVector zero(1);
+	zero = 0;
+	values = values.Rows(1,stopat) & zero;
+
+	ColumnVector gm(stopat + 1);
+	for(int j = 1; j <= stopat + 1; gm(j) = j++);
+	  
+	ColumnVector weights(stopat+1);
+	weights = 1;
+
+	bool anyviolators = true;
+
+	while(anyviolators) {
+	  anyviolators = false;
+	  
+	  for(int k = 2; k <= values.Nrows(); k++) {
+	    if(values(k) > values(k-1)) {
+	      anyviolators = true;
+	      values(k-1) = (values(k-1)*weights(k-1) + values(k)*weights(k))/(weights(k-1) + weights(k));
+	      values = values.Rows(1,k-1) & values.Rows(k+1,values.Nrows());
+	      weights(k-1) = weights(k) + weights(k-1);
+	      weights = weights.Rows(1,k-1) & weights.Rows(k+1,weights.Nrows());
+
+	      for(int j = 1; j <= stopat + 1; j++) {
+		if(gm(j) >= k)  
+		  gm(j) = gm(j)-1;
+	      }
+
+	      break;
+	    }
+	  }
+	}
+	
+	acEst.Column(i) = 0.0;
+	int j=1;
+
+	for(; j <= stopat; j++) {
+	  
+	  acEst(j,i) = values(gm(j));
+	  if(acEst(j,i) <= 0.0)
+	    {
+	      acEst(j,i) = 0.0;
+	      break;
+	    }
+	}
+	
+	if(acEst(2,i) < th/2)
+	  {
+	    acEst.SubMatrix(3,stopat,i,i) = 0;
+	  }
+	else if(j > 2)
+	  {
+	    int endst = j;
+	    int stst = j-(int)(1+(j/8.0));
+
+	    const int expwidth = MISCMATHS::Max((endst - stst)/2,1);
+	    const int exppow = 2;
+	
+	    for(j = stst; j <= endst; j++)
+	      {
+		acEst(j,i) = acEst(j,i)*exp(-pow((j-stst)/expwidth,exppow));
+	      }	
+	  }
+    }
+    cerr << " Completed" << endl;
+  }
+  
+  void AutoCorrEstimator::applyConstraints() {
+
+    Tracer tr("AutoCorrEstimator::applyConstraints");
+
+    cerr << "Applying constraints to AutoCorr estimates... ";
+
+    for(int i = 1; i <= xdata.getNumSeries(); i++)
+      {
+	int sizeTS = xdata.getNumVolumes();
+	int j = 3;
+	int stopat = (int)sizeTS/4;
+
+	// found1 is last valid value above threshold
+	int found1 = stopat;
+
+	// 5% point of distribution of autocorr about zero
+	const float thresh = (-1/sizeTS)+(2/sqrt(sizeTS));
+
+	acEst(2,i) = (acEst(2,i)+ acEst(3,i))/2;
+	if(acEst(2,i) < 0)
+	  {
+	    acEst(2,i) = 0;
+	  }
+	else
+	  {
+	    float grad = 0.0;
+
+	    while(j <= stopat && j < found1 + 2)
+	      {    
+		grad = ((acEst(j,i) + acEst(j+1,i))/2 - acEst(j-1,i))/1.5;
+		if(grad < 0)
+		  acEst(j,i) = grad + acEst(j-1,i);
+		else
+		  acEst(j,i) = acEst(j-1,i);
+		
+		// look for threshold
+		if(acEst(j,i) < thresh/3.0 && found1 == stopat)
+		  {
+		    found1 = j;
+		  }
+	   
+		if(acEst(j,i) < 0)
+		  {
+		    acEst(j,i) = 0;
+		  }
+		
+		j++;
+	      }
+	  }
+
+	// set rest to zero:
+	for(; j <= sizeTS; j++)
+	  {
+	      acEst(j,i) = 0;
+	  }
+      }
+    cerr << "Completed" << endl;
+  }
+
+  void AutoCorrEstimator::getMeanEstimate(ColumnVector& ret)
+    {
+      Tracer tr("AutoCorrEstimator::getMeanEstimate");
+
+      ret.ReSize(acEst.getNumVolumes());
+      // Calc global Vrow:
+      for(int i = 1; i <= acEst.getNumVolumes(); i++)
+	{
+	  ret(i) = MISCMATHS::mean(acEst.getVolume(i).AsColumn());
+	}
+    }
+}
+
+
+
+
+
+
+
+
+
+
+
+

deb

+   1.2246468e-16

film_gls.ccbak

+/*  film_gls.cc
+
+    Mark Woolrich, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#include <iostream>
+#include <fstream>
+#include <strstream>
+#define WANT_STREAM
+#define WANT_MATH
+
+#include "newmatap.h"
+#include "newmatio.h"
+#include "VolumeSeries.h"
+#include "Volume.h"
+#include "glim.h"
+#include "sigproc.h"
+#include "miscmaths.h"
+#include "gaussComparer.h"
+#include "Log.h"
+#include "AutoCorrEstimator.h"
+#include "paradigm.h"
+#include "FilmGlsOptions.h"
+#include "glimGls.h"
+
+#include <string>
+
+#ifndef NO_NAMESPACE
+using namespace NEWMAT;
+using namespace SIGPROC;
+using namespace TACO;
+using namespace UTILS;
+#endif
+
+int main(int argc, char *argv[])
+{
+  try{
+    rand();
+    // parse command line to find out directory name for logging:
+    ofstream out2;
+    FilmGlsOptions& globalopts = FilmGlsOptions::getInstance();
+    globalopts.parse_command_line(argc, argv, out2);
+
+    // Setup logging:
+    Log& logger = Log::getInstance();
+    logger.setLogFile("glslogfile");
+    logger.establishDir(globalopts.datadir);
+
+    // parse command line again to output arguments to logfile
+    globalopts.parse_command_line(argc, argv, logger.str());
+
+    // load non-temporally filtered data
+    VolumeSeries x;
+    x.read(globalopts.inputfname);
+
+    // if needed output the 12th volume for use later
+    Volume epivol;
+    if(globalopts.smoothACEst)
+      {
+	epivol = x.getVolume(12).AsColumn();
+	epivol.setDims(x.getDims());
+	
+	epivol.writeAsInt(logger.getDir() + "/" + globalopts.epifname);
+      }
+
+    // This also removes the mean from each of the time series:
+    x.thresholdSeries(globalopts.thresh, true);
+
+    // if needed later also threshold the epi volume
+    if(globalopts.smoothACEst)
+      {
+	epivol.setPreThresholdPositions(x.getPreThresholdPositions());
+	epivol.threshold();
+      }
+
+    int sizeTS = x.getNumVolumes();
+    int numTS = x.getNumSeries();
+    
+    // Load paradigm: 
+    Paradigm parad;
+    parad.load(globalopts.paradigmfname, "", false, sizeTS);
+
+    // Sort out detrending:
+    if(globalopts.detrend)
+      {
+	SIGPROC::Detrend(x, false);
+      }
+
+    if(globalopts.verbose)
+      {
+	logger.out("Gc", parad.getDesignMatrix());
+      }
+
+    // Set up OLS GLM for non-whitened data 
+    Glim glim(x, parad.getDesignMatrix());
+	    
+    cerr << "Computing residuals for non-whitened data... ";
+    VolumeSeries& rnotw = glim.ComputeResids();
+    cerr << "Completed" << endl;
+    
+    // Estimate Autocorrelations:
+    AutoCorrEstimator acEst(rnotw);
+
+    if(globalopts.fitAutoRegressiveModel)
+      {
+	acEst.fitAutoRegressiveModel();
+	if(globalopts.verbose)
+	  {
+	    AutoCorrEstimator acEstForLogging(rnotw);
+	    acEstForLogging.calcRaw();
+	    logger.out("rawac", acEstForLogging.getEstimates());
+	    logger.out("autoregac", acEst.getEstimates());
+	  }
+      }
+    else
+      {
+	acEst.calcRaw();
+	    
+	if(globalopts.verbose)
+	  {
+	    //logger.out("rawac", acEst.getEstimates());
+	  }
+	    
+	// Smooth raw estimates:
+	if(globalopts.smoothACEst)
+	  { 
+	    acEst.spatiallySmooth(logger.getDir() + "/" + globalopts.epifname, epivol, globalopts.ms, globalopts.epifname, globalopts.susanpath, globalopts.epith);
+	    
+	  }
+	
+	// Apply constraints to estimate autocorr:
+	acEst.pava();
+	    
+      }
+    
+    ColumnVector I(sizeTS);
+    I = 0;
+    I(1) = 1;
+    
+    // Setup OLS GLM for temporally filtered data:
+    int numParams = parad.getDesignMatrix().Ncols();
+    GlimGls glimGls(numTS, sizeTS, numParams);
+    
+    ColumnVector xw(sizeTS);
+    ColumnVector xprew(sizeTS);
+
+    acEst.setDesignMatrix(parad.getDesignMatrix());
+    int co = 1;
+    
+    // Loop through voxels:
+    cerr << "Calculating params..." << endl;
+
+    for(int i = 1; i <= numTS; i++)
+      {	
+	// Use autocorr estimate to prewhiten data:
+	xprew = x.getSeries(i);
+	    
+	Matrix designmattw;
+
+	acEst.preWhiten(xprew, xw, i, designmattw);
+
+	if(co > 1000 || i == 5618 || i == 5582)
+	  {
+	    co = 1;
+	    cerr << (float)i/(float)numTS << ",";
+	  }
+	else
+	  co++;
+	    
+	glimGls.setData(xw, designmattw, i);   
+      }
+    cerr << "Completed" << endl;
+    
+
+    if(globalopts.verbose)
+      {	
+	// Save E
+	VolumeSeries& E = acEst.getE();
+	ColumnVector& countLargeE = acEst.getCountLargeE();
+	
+	logger.out("countLargeE", countLargeE);
+	VolumeSeries::Dims dims = x.getDims();
+	dims.v = acEst.getZeroPad();
+
+	E.setDims(dims);
+	E.setPreThresholdPositions(x.getPreThresholdPositions());
+	E.unthresholdSeries();	
+	E.writeAsFloat(logger.getDir() + "/E");
+  
+	VolumeSeries& threshac = acEst.getEstimates();
+	threshac.unthresholdSeries(x.getDims(),x.getPreThresholdPositions());
+	threshac.writeAsFloat(logger.getDir() + "/threshac");
+	  
+	rnotw.unthresholdSeries(x.getDims(),x.getPreThresholdPositions());
+	rnotw.writeAsFloat(logger.getDir() + "/rnotw");
+      }
+
+    // Write out necessary data:
+    cerr << "Saving results... ";
+    glimGls.Save(x.getDims(), x.getPreThresholdPositions());
+    cerr << "Completed" << endl;
+  }  
+  catch(Exception p_excp) 
+    {
+      cerr << p_excp.what() << endl;
+    }
+  catch(...) 
+    {
+      cerr << "Image error" << endl;
+    } 
+  return 0;
+}
+

glim.baccc

+/*  glim.cc
+
+    Mark Woolrich, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#include <strstream>
+
+#include "glim.h"
+#include "miscmaths.h"
+#include "Log.h"
+#include "Volume.h"
+
+#ifndef NO_NAMESPACE
+using namespace MISCMATHS;
+using namespace UTILS;
+using namespace TACO;
+namespace SIGPROC {
+#endif
+
+  Glim::Glim(const VolumeSeries& p_y, const Matrix& p_x):
+    y(p_y),
+    x(p_x),
+    numTS(p_y.Ncols()),
+    sizeTS(p_y.Nrows()),
+    numParams(p_x.Ncols()),
+    r(sizeTS, numTS, p_y.getDims(), p_y.getPreThresholdPositions()),
+    pinv_x(p_x.Ncols(), sizeTS),
+    V(sizeTS,sizeTS),
+    RV(sizeTS,sizeTS),
+    RMat(sizeTS,sizeTS),
+    batch_size(BATCHSIZE),
+    corrections(numTS, numParams*numParams),
+    b(numTS, numParams),
+    sigmaSquareds(numTS),
+    dofs(numTS),
+    traceRVs(numTS)
+    { 
+    }
+
+  void Glim::Save()
+    {
+      // Need to save b, sigmaSquareds, corrections and dof 
+      Log& logger = Log::getInstance();
+     
+      // b:
+      Volume peVol;
+      for(int i = 1; i <= numParams; i++)
+	{
+	  peVol = b.Row(i).AsColumn();
+	  peVol.setDims(y.getDims());
+	  peVol.setPreThresholdPositions(y.getPreThresholdPositions());
+	  peVol.unthreshold();	
+	  
+	  // Add param number to "pe" to create filename:
+	  char strc[4];
+	  ostrstream osc(strc,4);
+	  osc << i << '\0';
+	  
+	  peVol.writeAsFloat(logger.getDir() + "/pe" + strc);
+	}
+
+      // sigmaSquareds:
+      sigmaSquareds.setDims(y.getDims());
+      sigmaSquareds.setPreThresholdPositions(y.getPreThresholdPositions());
+      sigmaSquareds.unthreshold();	
+      sigmaSquareds.writeAsFloat(logger.getDir() + "/sigmasquareds");
+
+      // dof:
+      dofs.setDims(y.getDims());
+      dofs.setPreThresholdPositions(y.getPreThresholdPositions());
+      dofs.unthreshold();	
+      dofs.writeAsFloat(logger.getDir() + "/dofs");
+
+      // traceRVs:
+      traceRVs.setDims(y.getDims());
+      traceRVs.setPreThresholdPositions(y.getPreThresholdPositions());
+      traceRVs.unthreshold();	
+      traceRVs.writeAsFloat(logger.getDir() + "/traceRVs");
+
+      // corrections:
+      logger.out("corrections", corrections);
+    }
+
+  // Called on entire data set:
+  const VolumeSeries& Glim::ComputeResids()
+    {
+      Tracer ts("ComputeResids");
+
+      int batch_pos = 1;
+
+      // pinv(x) = inv(x'x)x'
+      pinv_x = (x.t()*x).i()*x.t();
+
+      // R = I - x*pinv(x)
+      Matrix I(sizeTS, sizeTS);
+      Identity(I);
+
+      RMat = I - x*pinv_x;
+      
+      while(batch_pos <= numTS)
+	{
+	  if(batch_pos+batch_size - 1 > numTS)
+	    r.Columns(batch_pos, numTS) = RMat*y.Columns(batch_pos, numTS);
+	  else
+	    r.Columns(batch_pos, batch_pos+batch_size-1) = RMat*y.Columns(batch_pos, batch_pos+batch_size-1);
+	
+	  batch_pos += batch_size;
+	}
+      
+      return r;
+    }
+
+  // Called on entire data set:
+  void Glim::ComputePes()
+    { 
+      Tracer ts("ComputePe");
+      
+      b = pinv_x*y;
+    }
+
+  void Glim::ConstructV(const ColumnVector& p_vrow)
+    {
+      Tracer ts("ConstructV");
+      V = 0;
+
+      for (int i = 1; i <= sizeTS; i++)
+	{
+	  V.SubMatrix(i,i,i,sizeTS) = p_vrow.Rows(1,sizeTS-i+1).t();
+	  V.SubMatrix(i,i,1,i) = p_vrow.Rows(1,i).Reverse().t();
+	}
+    }
+
+  void Glim::SetVrow(const ColumnVector& p_vrow, const int ind)
+    {
+      Tracer ts("SetVrow");
+      
+      ConstructV(p_vrow);
+
+      // var/e = inv(x'x)x'*V*x*inv(x'x)
+      Matrix corr = pinv_x*V*pinv_x.t();
+      SetCorrection(corr, ind);
+      RV = RMat*V;
+      traceRVs(ind) = Trace(RV);
+
+      dofs(ind) = traceRVs(ind)*traceRVs(ind)/Trace(RV*RV);
+    }
+
+  void Glim::SetGlobalVrow(const ColumnVector& p_vrow)
+    {
+      Tracer ts("Glim::SetGlobalVrow");
+
+      ConstructV(p_vrow);
+      RV = RMat*V;
+
+      traceRVs = Trace(RV);
+ 
+      dofs = Trace(RV)*Trace(RV)/Trace(RV*RV);
+    }
+
+
+  void Glim::UseGlobalVrow()
+    {
+      Tracer ts("Glim::UseGlobalVrow");
+
+      // var/e = inv(x'x)x'*V*x*inv(x'x)
+      Matrix corr = pinv_x*V*x*(x.t()*x).i();
+
+      for(int i = 1; i <= numTS; i++)
+	{
+	  SetCorrection(corr, i);
+	}
+    }
+
+  void Glim::ComputeSigmaSquared(const int ind)
+    {
+      Tracer ts("Glim::ComputeSigmaSquared");
+
+      sigmaSquareds(ind) = (r.Column(ind).t()*r.Column(ind)/traceRVs(ind)).AsScalar();
+    }
+
+  void Glim::SetCorrection(const Matrix& corr, const int ind)
+    {
+      Tracer ts("SetCorrection");
+
+      // puts Matrix corr which is p*p into Matrix correction
+      // as a p*p length row:
+      int p = corr.Nrows();
+      
+      for (int i = 1; i <= p; i++)
+	{
+	  for (int j = 1; j <= p; j++)
+	    {
+	      corrections(ind, (i-1)*p + j) = corr(i,j); 
+	    }
+	}
+    }
+
+#ifndef NO_NAMESPACE
+}
+#endif
+
+
+
+
+
+
+

glim.bach

+/*  glim.h
+
+    Mark Woolrich, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#include <iostream>
+#include <fstream>
+#define WANT_STREAM
+#define WANT_MATH
+
+#include "newmatap.h"
+#include "newmatio.h"
+#include "VolumeSeries.h"
+#include "Volume.h"
+
+#ifndef NO_NAMESPACE
+using namespace NEWMAT;
+using namespace TACO;
+namespace SIGPROC{
+#endif
+
+#if !defined(__glim_h)
+#define __glim_h
+ 
+#define FALSE 0
+#define TRUE 1
+
+#define BATCHSIZE 50
+
+  class Glim
+    {
+    public:
+      Glim(const VolumeSeries& p_y, const Matrix& p_x); 
+
+      void Save();
+
+      const VolumeSeries& ComputeResids();      
+      void ComputePes();
+      void SetVrow(const ColumnVector& p_vrow, const int ind);
+      void SetGlobalVrow(const ColumnVector& p_vrow);
+      void ComputeSigmaSquared(const int ind);
+      void UseGlobalVrow();
+
+    private:
+      Glim();
+      Glim(const Glim&);
+      Glim& operator=(const Glim& p_glim);
+      
+      void SetCorrection(const Matrix& corr, const int ind);
+      void ConstructV(const ColumnVector& p_vrow);
+
+      // y = bx + r
+      const VolumeSeries& y;
+      const Matrix& x;
+
+      int numTS;
+      int sizeTS;
+      int numParams;
+
+      VolumeSeries r;
+      Matrix pinv_x;
+      Matrix V;
+      Matrix RV;
+      Matrix RMat;
+
+      int batch_size;
+
+      // Data to be saved:
+      Matrix corrections;
+      Matrix b;
+      Volume sigmaSquareds;
+      Volume dofs;
+      Volume traceRVs;
+
+    };
+
+
+#ifndef NO_NAMESPACE
+       }
+#endif
+
+#endif
+
+
+
+

ols.ccbak

+/*  ols.cc
+
+    Mark Woolrich, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#include "ols.h"
+#include "miscmaths.h"
+#include "Log.h"
+
+#ifndef NO_NAMESPACE
+using namespace MISCMATHS;
+namespace FILM {
+#endif
+
+  Ols::Ols(const Matrix& p_y, const Matrix& p_x, const Matrix& p_contrasts):
+    y(p_y),
+    x(p_x),
+    contrasts(p_contrasts),
+    numTS(p_y.Ncols()),
+    sizeTS(p_y.Nrows()),
+    r(sizeTS,numTS),
+    pinv_x(p_x.Ncols(), sizeTS),
+    var_on_e(0.0),
+    cb(numTS),
+    b(p_x.Ncols()),
+    var(numTS),
+    dof(sizeTS - p_x.Ncols()),
+    V(sizeTS,sizeTS),
+    RV(sizeTS,sizeTS),
+    RMat(sizeTS,sizeTS),
+    batch_size(BATCHSIZE)
+    {
+      SetContrast(1);
+    }
+  
+//    void Ols::ConstructV(const ColumnVector& p_vrow)
+//      {
+//        Tracer ts("ConstructV");
+//        V = 0;
+
+//        for (int i = 1; i <= sizeTS; i++)
+//         {
+//  	 V.SubMatrix(i,i,i,sizeTS) = p_vrow.Rows(1,sizeTS-i+1).t();
+//  	 V.SubMatrix(i,i,1,i) = p_vrow.Rows(1,i).Reverse().t();
+//         }
+//      }
+
+//    float Ols::SetupWithV(const ColumnVector& p_vrow, bool p_justvarone)
+//      {
+//        Tracer ts("SetupWithV");
+      
+//        ConstructV(p_vrow);
+
+//        // var/e = c'inv(x'x)x'*V*x*inv(x'x)*c
+//        var_on_e = (c.t()*pinv_x*V*x*(x.t()*x).i()*c).AsScalar();
+
+//        if(!p_justvarone)
+//  	{
+//  	  // dof = 2*trace(RV)^2/trace(R*V*R*V);
+//  	  RV = RMat*V;
+//  	  dof = Trace(RV)*Trace(RV)/Trace(RV*RV);
+//  	}
+
+//        return dof;
+//      }
+
+//    float Ols::SetupWithV(const ColumnVector& p_vrow, const ColumnVector& p_kfft, ColumnVector& vrow, bool p_justvarone, const int zeropad)
+//      {
+//        Tracer ts("SetupWithV");
+
+//        // make sure p_vrow is cyclic (even function)
+//        //ColumnVector vrow(zeropad);
+//        vrow.ReSize(zeropad);
+
+//        vrow = 0;
+//        vrow.Rows(1,sizeTS/2) = p_vrow.Rows(1,sizeTS/2);
+//        vrow.Rows(zeropad - sizeTS/2 + 2, zeropad) = p_vrow.Rows(2, sizeTS/2).Reverse();
+
+//        // fft vrow
+//        ColumnVector fft_real;
+//        ColumnVector fft_im;
+//        ColumnVector dummy(zeropad);
+//        dummy = 0;
+      
+//        ColumnVector realifft(zeropad);
+
+//        FFT(vrow, dummy, fft_real, fft_im);
+
+//        FFTI(SP(fft_real, p_kfft), dummy, realifft, dummy);
+      
+//        vrow = realifft.Rows(1,sizeTS);
+
+//        // Normalise vrow:
+//        vrow = vrow/vrow(1);
+
+//        ConstructV(vrow);     
+
+//        // var/e = c'inv(x'x)x'*V*x*inv(x'x)*c
+//        var_on_e = (c.t()*pinv_x*V*x*(x.t()*x).i()*c).AsScalar();
+      
+//        if(!p_justvarone)
+//  	{
+//  	  // dof = 2*trace(RV)^2/trace(R*V*R*V);
+//  	  RV = RMat*V;
+//  	  dof = Trace(RV)*Trace(RV)/Trace(RV*RV);
+//  	}
+
+//        return dof;
+//      }
+
+  const Matrix& Ols::ComputeResids()
+    {
+      Tracer ts("ComputeResids");
+
+      int batch_pos = 1;
+
+      // pinv(x) = inv(x'x)x'
+      pinv_x = (x.t()*x).i()*x.t();
+
+      // R = I - x*pinv(x)
+      Matrix I(sizeTS, sizeTS);
+      Identity(I);
+
+      RMat = I - x*pinv_x;
+      
+      while(batch_pos <= numTS)
+	{
+	  if(batch_pos+batch_size - 1 > numTS)
+	    r.Columns(batch_pos, numTS) = RMat*y.Columns(batch_pos, numTS);
+	  else
+	    r.Columns(batch_pos, batch_pos+batch_size-1) = RMat*y.Columns(batch_pos, batch_pos+batch_size-1);
+	
+	  batch_pos += batch_size;
+	}
+      
+      return r;
+    }
+
+//    const ColumnVector& Ols::Computecb()
+//      { 
+//        Tracer ts("Computecb");
+      
+//        //     cerr << "Computing cbs";
+//        int batch_pos = 1;
+      
+//        while(batch_pos <= numTS)
+//  	{
+//  	  if(batch_pos+batch_size - 1 > numTS)
+//  	    cb.Rows(batch_pos, numTS) = (c.t()*pinv_x*y.Columns(batch_pos, numTS)).t();
+//  	  else
+//  	    cb.Rows(batch_pos, batch_pos+batch_size-1) = (c.t()*pinv_x*y.Columns(batch_pos, batch_pos+batch_size-1)).t();
+//  	  batch_pos += batch_size;
+//  	  //	  cerr << ".";
+//  	}
+//        //cerr << endl;
+//        return cb;
+//      }
+
+//    const float Ols::Computecb(const int ind)
+//      { 
+//        Tracer ts("Computecb");
+
+//        cb(ind) = ((c.t()*pinv_x*y.Column(ind)).t()).AsScalar();
+//        return cb(ind);
+//      }
+
+//    const ColumnVector& Ols::Computeb(const int ind)
+//      { 
+//        Tracer ts("Computeb");
+      
+//        b = pinv_x*y.Column(ind);
+//        return b;
+//      }
+
+//    const ColumnVector& Ols::ComputeVar()
+//      { 
+//        Tracer ts("ComputeVar");
+
+//        //      cerr << "Computing Vars";
+//        int batch_pos = 1;
+
+//        Matrix varmatfull(batch_size, batch_size);
+//        ColumnVector vartempfull(batch_size);
+
+//        while(batch_pos <= numTS)
+//  	{
+//  	  if(batch_pos+batch_size - 1 > numTS)
+//  	    {
+//  	      // var = e*var_on_e
+//  	      // e is the estimate of the variance of the timeseries, sigma^2
+//  	      Matrix varmat = (r.Columns(batch_pos, numTS).t()*r.Columns(batch_pos, numTS)/Trace(RV))*var_on_e;
+//  	      ColumnVector vartemp;
+//  	      //getdiag(vartemp, varmat);  // obsolete fn
+//  	      vartemp = diag(varmat);  // MJ NOTE: new fn
+//  	      var.Rows(batch_pos, numTS) = vartemp;
+//  	    }      
+//  	  else
+//  	    {
+//  	      varmatfull = (r.Columns(batch_pos, batch_pos+batch_size-1).t()*r.Columns(batch_pos, batch_pos+batch_size-1)/Trace(RV))*var_on_e;
+//  	      //getdiag(vartempfull, varmatfull);  // obsolete fn
+//  	      vartempfull = diag(varmatfull);  // MJ NOTE: new fn
+//  	      var.Rows(batch_pos, batch_pos+batch_size-1) = vartempfull;
+//  	    }
+//  	  batch_pos += batch_size;
+//  	  //  cerr << ".";
+//  	}
+      
+//        // cerr << endl;
+//        return var;
+//      }
+
+//    const float Ols::ComputeVar(const int ind)
+//      { 
+//        Tracer ts("ComputeVar");
+   
+//        // var = e*var_on_e
+//        // e is the estimate of the variance of the timeseries, sigma^2
+//        var(ind) = ((r.Column(ind).t()*r.Column(ind)/Trace(RV))*var_on_e).AsScalar();
+//        // var(ind) = (r.Column(ind).t()*r.Column(ind)*var_on_e).AsScalar();
+
+//        return var(ind);
+//      }
+
+#ifndef NO_NAMESPACE
+}
+#endif
+
+
+
+
+
+
+
+
+

ols.hbak

+/*  ols.h
+
+    Mark Woolrich, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#if !defined(__ols_h)
+#define __ols_h
+
+#include <iostream>
+#include <fstream>
+#define WANT_STREAM
+#define WANT_MATH
+
+#include "newmatap.h"
+#include "newmatio.h"
+
+using namespace NEWMAT;
+namespace FILM {
+ 
+#define BATCHSIZE 50
+
+  class Ols
+    {
+    public:
+     
+      Ols(const Matrix& p_y, const Matrix& p_x, const Matrix& p_contrasts); 
+
+      const Matrix& ComputeResids();
+      const ColumnVector& ComputeVar();
+      const ColumnVector& Computecb();
+      const float ComputeVar(const int ind);
+      const float Computecb(const int ind);
+      const ColumnVector& Computeb(const int ind);
+
+      float SetupWithV(const ColumnVector& p_vrow, bool p_justvarone = false);
+      float SetupWithV(const ColumnVector& p_vrow, const ColumnVector& p_krowspec, ColumnVector& p_res, bool p_justvarone = false, const int zeropad = 512);
+
+      const Matrix& Getr() const {return r;}
+     
+      const ColumnVector& Getvar() const {return var;}
+      const ColumnVector& Getcb() const {return cb;}
+      float Getdof() const {return dof;}
+      float Getvarone() const {return var_on_e;}
+
+      void SetContrast(const int p_num) {c = contrasts.Row(p_num).t();}
+      void Setvarone(const float p_varone) {var_on_e = p_varone;}
+      
+    private:
+      Ols();
+      Ols(const Ols&);
+      Ols& operator=(const Ols& p_ols);
+       
+      void ConstructV(const ColumnVector& p_vrow);
+
+      // y = bx + r
+      const Matrix& y;
+      const Matrix& x;
+   
+      // Contrast
+      const Matrix& contrasts;
+      ColumnVector c;
+
+      int numTS;
+      int sizeTS;
+
+      Matrix r;
+
+      Matrix pinv_x;
+      float var_on_e;
+
+      ColumnVector cb;
+      
+      ColumnVector b;
+
+      ColumnVector var;
+      float dof;
+
+      // Covariance matrix is sigma^2*V:
+      Matrix V;
+      Matrix RV;
+
+      Matrix RMat;
+      
+      // batch stuff
+      int batch_size;
+    };
+
+}
+
+#endif
+
+
+

sigprocbak.ccbak

+/*  sigproc.cc
+
+    Mark Woolrich, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#include "ols.h"
+#include "sigproc.h"
+#include "miscmaths.h"
+#include "glm.h"
+
+#ifndef NO_NAMESPACE
+namespace SIGPROC {
+  using namespace MISCMATHS;
+#endif
+
+  void Out(string p_fname, const Matrix& p_mat)
+    {
+      ofstream out;
+      out.open(p_fname.c_str(), ios::out);
+      out << "/Matrix" << endl;
+      out << p_mat;
+      out.close();
+    }
+
+  void Out(string p_fname, const ColumnVector& p_mat)
+    {
+      ofstream out;
+      out.open(p_fname.c_str(), ios::out);
+      out << p_mat;
+      out.close();
+    }
+
+  void Out(string p_fname, const RowVector& p_mat)
+    {
+      ofstream out;
+      out.open(p_fname.c_str(), ios::out);
+      out << p_mat;
+      out.close();
+    }
+
+  void In(string p_fname, ColumnVector& p_mat)
+    {
+      int size = p_mat.Nrows();
+
+      ifstream in;
+      in.open(p_fname.c_str(), ios::in);
+      
+      if(!in)
+	{
+	  string err("Unable to open ");
+	  err =  err +  p_fname;
+	  throw Exception(err.c_str());
+	}
+
+      for(int i = 1; i <= size; i++)
+	{
+	  in >> p_mat(i);
+	}
+
+      in.close();
+    }
+
+  void In(string p_fname, Matrix& p_mat)
+    {
+      ifstream in;
+      in.open(p_fname.c_str(), ios::in);
+      
+      if(!in)
+	{
+	  cerr << "Unable to open " << p_fname << endl;
+	  throw;
+	}
+
+      int numWaves = 0;
+      int numPoints = 0;
+       
+      string str;
+
+      while(true)
+	{
+	  in >> str;
+	  if(str == "/Matrix")
+	    break;
+	  else if(str == "/NumWaves")
+	    {
+	      in >> numWaves;
+	    }
+	  else if(str == "/NumPoints" || str == "/NumContrasts")
+	    {
+	      in >> numPoints;
+	    }
+	}
+
+      if(numWaves != 0)
+	{
+	  p_mat.ReSize(numPoints, numWaves);
+	}
+      else
+	{
+	  numWaves = p_mat.Ncols();
+	  numPoints = p_mat.Nrows();
+	}
+      
+      for(int i = 1; i <= numPoints; i++)
+	{
+	  for(int j = 1; j <= numWaves; j++)    
+	    {
+	      in >> p_mat(i,j);
+	    }
+	  }
+
+      in.close();
+    }
+
+  int EstablishZeroPadding(int sizeTS) {
+    
+    Tracer tr("EstablishZeroPadding");
+    
+    return (int)pow(2,ceil(log(sizeTS)/log(2)));
+  }
+
+  void BandPass(Matrix& p_ts, int lowcut, int highcut)
+  {
+    Tracer tr("BandPass");
+
+      int sizeTS = p_ts.Nrows();
+      int numTS = p_ts.Ncols();
+      
+      int zeropad = EstablishZeroPadding(sizeTS);
+      ColumnVector x(zeropad);
+      x = 0;
+      ColumnVector fft_real;
+      ColumnVector fft_im;
+      ColumnVector dummy(zeropad);
+      ColumnVector dummy2;
+      dummy = 0;
+      ColumnVector realifft(zeropad);
+
+      // convert cutofs to zeropad cutoffs:
+      lowcut = (int)(lowcut*zeropad/(float)sizeTS);
+      highcut = (int)(highcut*zeropad/(float)sizeTS);
+
+      for(int i = 1; i <= numTS; i++)
+	{
+	  x.Rows(1,sizeTS) = p_ts.Column(i);
+
+	  FFT(x, dummy, fft_real, fft_im);
+
+	  // Perform cutoff here:
+	  // everything greater than highcut and lower than lowcut is removed
+	  if(lowcut > 0 && lowcut < sizeTS){
+	    fft_real.Rows(1, lowcut) = 0;
+	    fft_im.Rows(1, lowcut) = 0;
+	    fft_real.Rows(zeropad - lowcut, zeropad) = 0;
+	    fft_im.Rows(zeropad - lowcut, zeropad) = 0;
+	  }
+	  if(highcut < sizeTS && highcut > 0){
+	    fft_real.Rows(highcut, zeropad - highcut) = 0;
+	    fft_im.Rows(highcut, zeropad - highcut) = 0;
+	  }
+	  //////////////////////
+
+	  FFTI(fft_real, fft_im, realifft, dummy2);
+	  p_ts.Column(i) = realifft.Rows(1,sizeTS);
+	} 
+  }
+
+  void AutoCorr(const Matrix& p_ts, Matrix& p_ret, Matrix& fts, int p_zeropad)
+    {
+      Tracer tr("AutoCorr");
+
+      int sizeTS = p_ts.Nrows();
+      int numTS = p_ts.Ncols();
+
+      ColumnVector x(p_zeropad);
+      x = 0;
+      ColumnVector fft_real;
+      ColumnVector fft_im;
+      ColumnVector dummy(p_zeropad);
+      ColumnVector dummy2;
+      dummy = 0;
+      ColumnVector realifft(p_zeropad);
+      p_ret.ReSize(sizeTS,numTS);
+      fts.ReSize(p_zeropad,numTS);
+
+      for(int i = 1; i <= numTS; i++)
+	{
+	  x.Rows(1,sizeTS) = p_ts.Column(i);
+	  FFT(x, dummy, fft_real, fft_im);
+
+	  for(int j = 1; j <= p_zeropad; j++)
+	    {
+	      // (x+iy)(x-iy) = x^2 + y^2
+	      fft_real(j) = fft_real(j)*fft_real(j) + fft_im(j)*fft_im(j);
+	      fft_im(j) = 0;
+	      fts(j,i) = fft_real(j)*fft_real(j) + fft_im(j)*fft_im(j);      
+	    }
+
+	  // normalise fts:
+	  fts.Column(i) = fts.Column(i)/fts.Column(i).Sum();
+	  
+	  FFTI(fft_real, fft_im, realifft, dummy2);
+
+	  float varx = var(ColumnVector(x.Rows(1,sizeTS)));
+	  p_ret.Column(i) = realifft.Rows(1,sizeTS);
+ 
+	  for(int j = 1; j <= sizeTS-1; j++)
+	    {
+	       // Correction to make autocorr unbiased and normalised
+	       p_ret(j,i) = p_ret(j,i)/((sizeTS-j)*varx);
+	    }  
+	}
+
+    }
+
+  int Pacf(const ColumnVector& x, int minorder, int maxorder, ColumnVector& betas)
+    { 
+      Tracer ts("pacf");
+      int order = -1;
+      int sizeTS = x.Nrows();
+
+      // Set c
+      Matrix c(1,1);
+      c(1,1) = 1;
+      
+      Glm glm;
+
+      for(int i = minorder+1; i <= maxorder+1; i++)
+	{
+	  ColumnVector y = x.Rows(i+1,sizeTS);
+
+	  // Setup design matrix
+	  Matrix X(sizeTS-i, i);
+	  X = 0;
+	  
+	  for(int j = 1; j <= i; j++)
+	    {
+	      X.Column(j) = x.Rows(i+1-j,sizeTS-j).AsColumn();	    
+	    }
+	  
+	  glm.setData(y, X, c);
+	  
+	  glm.ComputeResids();
+	  betas = glm.Getb();
+	  
+	  if((abs(betas(i)) < (1/sizeTS) + (2/pow(sizeTS,0.5)) && order == -1) 
+	     || i == maxorder+1)
+          {      
+	    order = i-1;
+	    break;
+	  }
+	}
+      return order; 
+    }
+
+  void Detrend(Matrix& p_ts, bool p_justremovemean)
+    {
+      Tracer trace("SIGPROC::Detrend");
+
+      cerr << "Detrending...";
+      int sizeTS = p_ts.Nrows();
+     
+      // Set c
+      Matrix c(1,1);
+      c(1,1) = 1;
+
+      // p_ts = b*a + p_ret OLS regression
+      Matrix a(sizeTS, 3);
+
+      if(p_justremovemean) a.ReSize(sizeTS,1);
+
+      // Create a
+      for(int i = 1; i <= sizeTS; i++)
+	{
+	  a(i,1) = 1;
+	  if(!p_justremovemean)
+	    {
+	      a(i,2) = (float)i/sizeTS;
+	      a(i,3) = pow((float)i/sizeTS,2);
+	    }
+	}
+          
+      // Do OLS for all TS:
+      Ols ols(p_ts, a, c);
+
+      p_ts = ols.ComputeResids();
+      
+      cerr << " Completed" << endl;
+    }
+
+#ifndef NO_NAMESPACE
+}
+#endif
+
+
+
+
+
+

sigprocbak.hbak

+/*  sigproc.h
+
+    Mark Woolrich, FMRIB Image Analysis Group
+
+    Copyright (C) 1999-2000 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#include <iostream>
+#include <fstream>
+#define WANT_STREAM
+#define WANT_MATH
+#include <string>
+
+#include "newmatap.h"
+#include "newmatio.h"
+  
+#if !defined(__sigproc_h)
+#define __sigproc_h
+  
+using namespace NEWMAT;
+namespace SIGPROC{
+
+  int EstablishZeroPadding(int num);
+  void BandPass(Matrix& p_ts, int lowcut, int highcut);
+  void AutoCorr(const Matrix& p_ts, Matrix& p_ret, Matrix& fts, int p_zeropad);
+  void Detrend(Matrix& p_ts, bool p_justremovemean = false);
+  int Pacf(const ColumnVector& x, int minorder, int maxorder, ColumnVector& betas);
+  void Out(string p_fname, const Matrix& out);
+  void Out(string p_fname, const ColumnVector& out);
+  void Out(string p_fname, const RowVector& out);
+  void In(string p_fname, ColumnVector& p_mat);
+  void In(string p_fname, Matrix& p_mat); 
+
+}
+#endif
+
+

tmp

+x = sin(pi)
+save -ASCII deb x
\ No newline at end of file