diff --git a/fsl_glm.cc b/fsl_glm.cc
index f6801ccb6b889ffe4265a728e86f21719348af01..3d97261b161b80a76898d59f5613f0c233895c16 100644
--- a/fsl_glm.cc
+++ b/fsl_glm.cc
@@ -48,7 +48,7 @@ using namespace std;
   Option<string> fnftest(string("-f,--ftests"), string(""),
 		string("matrix of F-tests on contrasts"),
 		false, requires_argument,false);
-	Option<int> dofset(string("--dof"),0,
+	Option<int> dofset(string("--dof"),-1,
 		string("        set degrees-of-freedom explicitly"),
 		false, requires_argument);
 	Option<bool> normdes(string("--des_norm"),FALSE,
@@ -103,6 +103,9 @@ using namespace std;
 	Option<string> outvnscales(string("--out_vnscales"),string(""),
 		string("output file name for scaling factors for variance normalisation"),
 		false, requires_argument);
+        Option<vector<string> > voxelwiseConfounds(string("--vxf"), vector<string>(), 
+         string("\tlist of 4D images containing voxelwise EVs (list order corresponds to numbers in vxl option). caution BETA option."), 
+         false, requires_argument);
 		/*
 }
 */
@@ -147,7 +150,7 @@ void saveit(Matrix what, string fname){
 		write_ascii_matrix(what,fname);
 }
 
-int setup(){
+int setup(int &dof){
 	if(fsl_imageexists(fnin.value())){//read data
 		//input is 3D/4D vol
 		volume4D<float> tmpdata;
@@ -199,10 +202,34 @@ int setup(){
 		design = read_ascii_matrix(fndesign.value());
 	}
 
+        dof=(int)ols_dof(design);
+
+	if ( voxelwiseConfounds.set() ) {
+	  vector<Matrix> confounds;
+	  confounds.resize(voxelwiseConfounds.value().size());
+	  cerr << confounds.size() << endl;
+	  volume4D<float> input;
+	  for(unsigned int i=0; i< confounds.size(); i++) {
+	    read_volume4D(input,voxelwiseConfounds.value().at(i));
+	    if ( mask.nvoxels() )	  
+	      confounds.at(i)=input.matrix(mask);
+	    else
+	    confounds.at(i)=input.matrix();
+	  }
+	  for(int voxel=1;voxel<=data.Ncols();voxel++) {
+	    Matrix confound(confounds.at(0).Column(voxel) );
+	    for(unsigned int i=1; i< confounds.size(); i++) 
+	    confound|=confounds.at(i).Column(voxel);
+	    data.Column(voxel)=(IdentityMatrix(confound.Nrows())-confound*pinv(confound))*data.Column(voxel);	    
+	    }
+	  dof-=2;
+	}
+
 	if(perf_demean.value()){
 		if(debug.value())
 			cout << "De-meaning the data matrix" << endl;
 		data = remmean(data,1);
+		dof-=1;
 	}
 	if(normdat.value()){
 		if(debug.value())
@@ -262,12 +289,12 @@ void write_res(){
 }
 
 int do_work(int argc, char* argv[]) {
-  if(setup())
-		exit(1);
-
-	glm.olsfit(data,design,contrasts,dofset.value());
-	write_res();
-	return 0;
+  int dof(-1);
+  if(setup(dof))
+    exit(1);
+  glm.olsfit(data,design,contrasts,dof);
+  write_res();
+  return 0;
 }
 
 ////////////////////////////////////////////////////////////////////////////
@@ -302,6 +329,7 @@ int main(int argc,char *argv[]){
 			options.add(outsigsq);
 			options.add(outdata);
 			options.add(outvnscales);
+			options.add(voxelwiseConfounds);
 	    options.parse_command_line(argc, argv);
 
 	    // line below stops the program if the help was requested or 
diff --git a/melhlprfns.cc b/melhlprfns.cc
index 0ce17c03b5d785ca73f6d68fec05d80e4425f90b..7d5a013f134347379716e127e5bbe71b02b6b461 100644
--- a/melhlprfns.cc
+++ b/melhlprfns.cc
@@ -884,7 +884,7 @@ namespace Melodic{
   }  //Matrix gen_arCorr
 
 	void basicGLM::olsfit(const Matrix& data, const Matrix& design, 
-		const Matrix& contrasts, int DOFadjust)
+		const Matrix& contrasts, int requestedDOF)
 	{
 		beta = zeros(design.Ncols(),1); 
 		residu = zeros(1); sigsq = -1.0*ones(1); varcb = -1.0*ones(1); 
@@ -898,8 +898,9 @@ namespace Melodic{
 			
 			beta = pinvdes * dat;
 			residu = dat - design*beta;
-
-			dof = design.Nrows() - design.Ncols()-1;
+			dof = ols_dof(design);
+			if ( requestedDOF>0)
+			  dof = requestedDOF;
 			sigsq = sum(SP(residu,residu))/dof;
 			
 			float fact = float(dof) / design.Ncols();
diff --git a/melhlprfns.h b/melhlprfns.h
index 7b4dcbef5d24578b9708aa453071ac1593f14873..516445193d0be4d93c541ca1ec65c73378fe453d 100644
--- a/melhlprfns.h
+++ b/melhlprfns.h
@@ -86,7 +86,7 @@ namespace Melodic{
 			~basicGLM(){}
 		
 			void olsfit(const Matrix& data, const Matrix& design, 
-				const Matrix& contrasts, int DOFadjust = 0);
+				const Matrix& contrasts, int DOFadjust = -1);
 
 			inline Matrix& get_t(){return t;}
 			inline Matrix& get_z(){return z;}