diff --git a/AutoCorrEstimator.bak2 b/AutoCorrEstimator.bak2
new file mode 100755
index 0000000000000000000000000000000000000000..ea876a31398704b39370c2ad0f748aae1e583699
--- /dev/null
+++ b/AutoCorrEstimator.bak2
@@ -0,0 +1,536 @@
+#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());
+ }
+ }
+}
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/deb b/deb
new file mode 100755
index 0000000000000000000000000000000000000000..48ffaa2e39a2f4819b679cc4126f71fb24050398
--- /dev/null
+++ b/deb
@@ -0,0 +1 @@
+ 1.2246468e-16
diff --git a/film_gls.ccbak b/film_gls.ccbak
new file mode 100755
index 0000000000000000000000000000000000000000..7cb5e6d56b45b1b860887e7a846bbb9ed0dfa3b0
--- /dev/null
+++ b/film_gls.ccbak
@@ -0,0 +1,216 @@
+/* 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;
+}
+
diff --git a/glim.baccc b/glim.baccc
new file mode 100755
index 0000000000000000000000000000000000000000..e8e87d201cfbceb91faadce83169956f053177c5
--- /dev/null
+++ b/glim.baccc
@@ -0,0 +1,210 @@
+/* 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
+
+
+
+
+
+
+
diff --git a/glim.bach b/glim.bach
new file mode 100755
index 0000000000000000000000000000000000000000..41cfe361a7002a73964743d2380c7347204518fc
--- /dev/null
+++ b/glim.bach
@@ -0,0 +1,89 @@
+/* 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
+
+
+
+
diff --git a/logfile b/logfile
new file mode 100755
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/ols.ccbak b/ols.ccbak
new file mode 100755
index 0000000000000000000000000000000000000000..0cfe9730faba4812d91eac5dd217c3dd6a3cf463
--- /dev/null
+++ b/ols.ccbak
@@ -0,0 +1,238 @@
+/* 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
+
+
+
+
+
+
+
+
+
diff --git a/ols.hbak b/ols.hbak
new file mode 100755
index 0000000000000000000000000000000000000000..6755f729ec3dc2e706323e143c6756696c1cf881
--- /dev/null
+++ b/ols.hbak
@@ -0,0 +1,96 @@
+/* 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
+
+
+
diff --git a/sigprocbak.ccbak b/sigprocbak.ccbak
new file mode 100755
index 0000000000000000000000000000000000000000..fed675c4ee315dd83dedb1bf874e2a17f25a0fa0
--- /dev/null
+++ b/sigprocbak.ccbak
@@ -0,0 +1,302 @@
+/* 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
+
+
+
+
+
+
diff --git a/sigprocbak.hbak b/sigprocbak.hbak
new file mode 100755
index 0000000000000000000000000000000000000000..68a13540c7c9bb5e4e5386c931ed85c4740392f7
--- /dev/null
+++ b/sigprocbak.hbak
@@ -0,0 +1,38 @@
+/* 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
+
+
diff --git a/tmp b/tmp
new file mode 100755
index 0000000000000000000000000000000000000000..5f2f042bdd3ae5e3785fae363fa1f89517ac623b
--- /dev/null
+++ b/tmp
@@ -0,0 +1,2 @@
+x = sin(pi)
+save -ASCII deb x
\ No newline at end of file