diff --git a/fsl_glm.cc b/fsl_glm.cc
index e7ce0d0df26f1ae32dc33c59655b4f86b9f6f121..0398b2992ce5e9fb7f4b327ca72b86e5f979c866 100644
--- a/fsl_glm.cc
+++ b/fsl_glm.cc
@@ -27,14 +27,14 @@ using namespace std;
 		string(" \n Simple GLM usign ordinary least-squares regression on\n")+
 		string(" time courses and/or 3D/4D imges against time courses \n")+
 		string(" or 3D/4D images\n\n");
-  string examples="fsl_glm -i <input> -d <design> [options]";
+  string examples="fsl_glm -i <input> -d <design> -o <output> [options]";
 
 //Command line Options {
   Option<string> fnin(string("-i,--in"), string(""),
 		string("        input file name (matrix 3D or 4D image)"),
 		true, requires_argument);
   Option<string> fnout(string("-o,--out"), string(""),
-		string("       output file name for GLM parameter estimates"),
+		string("output file name for GLM parameter estimates"),
 		true, requires_argument);
   Option<string> fndesign(string("-d,--design"), string(""),
 		string("file name of the GLM design matrix (time courses or spatial maps)"),
@@ -44,7 +44,7 @@ using namespace std;
 		false, requires_argument);
   Option<string> fncontrasts(string("-c,--contrasts"), string(""),
 		string("matrix of t-statistics contrasts"),
-		false, requires_argument);
+		false, requires_argument, false);
   Option<string> fnftest(string("-f,--ftests"), string(""),
 		string("matrix of F-tests on contrasts"),
 		false, requires_argument);
@@ -131,6 +131,8 @@ bool isimage(Matrix what){
 void saveit(Matrix what, string fname){
 	if(isimage(what))
 		save4D(what,fname);
+	else if(fsl_imageexists(fndesign.value()))
+		write_ascii_matrix(what.t(),fname);
 	else
 		write_ascii_matrix(what,fname);
 }
diff --git a/meldata.cc b/meldata.cc
index c6235e09fe1f36261ec3998303c33741e195d74c..6d6fd082b305be2084b8688be173c7cb2f02cf66 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -45,7 +45,7 @@ namespace Melodic{
     tmpData = RawData.matrix(Mask);
     
     //estimate smoothness
-    if(Resels == 0)
+    if((Resels == 0)&&(!opts.filtermode))
       Resels = est_resels(RawData,Mask);
         
     //convert to percent BOLD signal change
@@ -62,7 +62,7 @@ namespace Melodic{
       message(" done" << endl);
     }
 
-    meanC = mean(tmpData,2);
+   // meanC = mean(tmpData,2);
 	//	tmpData = remmean(tmpData,2);
 
     //convert to power spectra
@@ -134,15 +134,24 @@ namespace Melodic{
 
   void MelodicData::set_TSmode()
   {
-    Matrix tmp, tmpT, tmpS, tmpT2, tmpS2;
+    Matrix tmp, tmpT, tmpS, tmpT2, tmpS2, tmpT3;
     tmp = expand_dimred(mixMatrix);
     tmpT = zeros(tmp.Nrows()/numfiles, tmp.Ncols());
     tmpS = zeros(numfiles, tmp.Ncols());
-    krfact(tmp,tmpT,tmpS);
+    explained_var = krfact(tmp,tmpT,tmpS);
+		outMsize("tmp",tmp);
+		outMsize("tmpT",tmpT);
+		outMsize("tmpS",tmpS);
+		
     Tmodes.clear(); Smodes.clear();
     for(int ctr = 1; ctr <= tmp.Ncols(); ctr++){
+			tmpT3 << reshape(tmp.Column(ctr),tmpT.Nrows(),numfiles);
+			outMsize("tmpT3", tmpT3);
       tmpT2 << tmpT.Column(ctr);
       tmpS2 << tmpS.Column(ctr);
+			tmpT3 << SP(tmpT3,pow(ones(tmpT3.Nrows(),1)*tmpS2.t(),-1));
+			if(numfiles>1)
+				tmpT2 |= tmpT3;
 			if(mean(tmpS2,1).AsScalar()<0){
 				tmpT2*=-1.0;
 				tmpS2*=-1.0;
@@ -153,9 +162,9 @@ namespace Melodic{
 
 		//add GLM OLS fit
 		if(Tdes.Storage()){
-			Matrix alltcs = Tmodes.at(0);
+			Matrix alltcs = Tmodes.at(0).Column(1);
 			for(int ctr=1; ctr < (int)Tmodes.size();ctr++)
-				alltcs|=Tmodes.at(ctr);
+				alltcs|=Tmodes.at(ctr).Column(1);
 			if((alltcs.Nrows()==Tdes.Nrows())&&(Tdes.Nrows()>Tdes.Ncols()))
 				glmT.olsfit(alltcs,Tdes,Tcon);
 		}
@@ -170,8 +179,14 @@ namespace Melodic{
 
   void MelodicData::setup()
   { 
-    setup_misc();
     numfiles = (int)opts.inputfname.value().size();
+    setup_misc();
+
+		if(opts.filtermode){ // basic setup for filtering only
+			Data = process_file(opts.inputfname.value().at(0));
+		}
+		else{
+
     if((numfiles > 1) && (opts.approach.value()==string("defl") || opts.approach.value()==string("symm")))
       opts.approach.set_T("tica");
 
@@ -288,15 +303,15 @@ namespace Melodic{
     message(endl << "  Data size : "<<Data.Nrows()<<" x "<<Data.Ncols()<<endl<<endl);
  		outMsize("stdDev",stdDev);
    
-    meanC=mean(Data,2);
+    //meanC=mean(Data,2);
 		if(opts.debug.value())
 			save4D(Data,"concat_data");    
     //save the mean & mask
     save_volume(Mask,logger.appendDir("mask"));
     save_volume(Mean,logger.appendDir("mean"));
-    
+    }
   } // void setup()
-
+	
   void MelodicData::setup_misc()
   {
 
@@ -348,6 +363,13 @@ namespace Melodic{
 		if(opts.fn_SconF.value().length()>0)
 			SconF = read_ascii_matrix(opts.fn_SconF.value());
 		
+		if(numfiles>1 && Sdes.Storage() == 0){
+			Sdes = ones(numfiles,1);
+			if(Scon.Storage() == 0){
+				Scon = ones(1,1);
+				Scon &= -1*Scon;
+			}
+		}
 		Tdes = remmean(Tdes,1);
   }
 
@@ -378,7 +400,6 @@ namespace Melodic{
     if((IC.Storage()>0)&&(opts.output_origIC.value())&&(after_mm==false))
       save4D(IC,opts.outputfname.value() + "_oIC");
       
-
     //Output IC -- adjusted for noise	
       if(IC.Storage()>0){
 				volume4D<float> tempVol;	
@@ -412,7 +433,7 @@ namespace Melodic{
 
     //Output mixMatrix
     if(mixMatrix.Storage()>0){
-      saveascii(mixMatrix, opts.outputfname.value() + "_mix");
+      saveascii(expand_mix(), opts.outputfname.value() + "_mix");
       mixFFT=calc_FFT(expand_mix(), opts.logPower.value());
       saveascii(mixFFT,opts.outputfname.value() + "_FTmix");      
     }
@@ -521,8 +542,7 @@ namespace Melodic{
 		outMsize("meanC",meanC);
 
     newData = Data - noiseMix * noiseIC.t();
-		if(meanC.Storage()>0)
-    	newData = newData + meanC*ones(1,newData.Ncols());
+
 		if(meanR.Storage()>0)
     	newData = newData + ones(newData.Nrows(),1)*meanR;
     
diff --git a/meldata.h b/meldata.h
index 809389e9ffe0a703b128a5c675097b434a27aa67..e20044dd10efb22f0af9fad977395fdb478bf759 100644
--- a/meldata.h
+++ b/meldata.h
@@ -36,7 +36,7 @@ namespace Melodic{
  
       void save();
 
-      Matrix process_file(string fname, int numfiles);
+      Matrix process_file(string fname, int numfiles = 1);
 
       inline void save4D(Matrix what, string fname){
 	 			volume4D<float> tempVol;
@@ -193,7 +193,8 @@ namespace Melodic{
       volumeinfo tempInfo;
       vector<Matrix> DWM, WM;
 			basicGLM glmT, glmS;
-			Matrix Tdes, Tcon, TconF, Sdes, Scon, SconF;			
+			Matrix Tdes, Tcon, TconF, Sdes, Scon, SconF;	
+			RowVector explained_var;		
 
     private:
       MelodicOptions &opts;     
diff --git a/melhlprfns.cc b/melhlprfns.cc
index e353681668569dc8d8b7623b582dbedf32c382e1..91cbf184fcc51ad85471604677181eca2a499650 100644
--- a/melhlprfns.cc
+++ b/melhlprfns.cc
@@ -316,7 +316,7 @@ namespace Melodic{
     evecs = C * Evc;
   }  //void em_pca
 
-  void rankapprox(const Matrix& Mat, Matrix& cols, Matrix& rows, int dim)
+  float rankapprox(const Matrix& Mat, Matrix& cols, Matrix& rows, int dim)
   { 
     Matrix Corr, Evecs, tmpWM, tmpDWM, tmp;
     RowVector Evals;
@@ -325,29 +325,36 @@ namespace Melodic{
     tmp = tmpWM * Mat.t();
     cols = tmp.t();
     rows << tmpDWM;
+		float res;
+		Evals=fliplr(Evals);
+		res = sum(Evals.Columns(1,dim),2).AsScalar()/sum(Evals,2).AsScalar()*100;
+		return res;
   } // rankapprox
 
-  void krfact(const Matrix& Mat, Matrix& cols, Matrix& rows)
+  RowVector krfact(const Matrix& Mat, Matrix& cols, Matrix& rows)
   {
-    Matrix out;
-    for(int ctr1 = 1; ctr1 <= Mat.Ncols(); ctr1++)
-      {
-	Matrix tmpVals(cols.Nrows(),rows.Nrows());
-	for(int ctr2 = 1; ctr2 <= rows.Nrows(); ctr2++)
-	  tmpVals.Column(ctr2) << Mat.SubMatrix(cols.Nrows() * (ctr2 - 1) + 1,cols.Nrows()*ctr2 ,ctr1,ctr1);
+		Matrix out; RowVector res(Mat.Ncols());
+    for(int ctr1 = 1; ctr1 <= Mat.Ncols(); ctr1++){
+			Matrix tmpVals(cols.Nrows(),rows.Nrows());
+			for(int ctr2 = 1; ctr2 <= rows.Nrows(); ctr2++)
+	  		tmpVals.Column(ctr2) << Mat.SubMatrix(cols.Nrows() * 
+				(ctr2 - 1) + 1,cols.Nrows()*ctr2 ,ctr1,ctr1);
 	
-	Matrix tmpcols, tmprows;
-	rankapprox(tmpVals, tmpcols, tmprows);
-	cols.Column(ctr1) = tmpcols;
-	rows.Column(ctr1) = tmprows;
-      }
+			Matrix tmpcols, tmprows;
+	 		res(ctr1) =rankapprox(tmpVals, tmpcols, tmprows);
+			cols.Column(ctr1) = tmpcols;
+			rows.Column(ctr1) = tmprows;
+    }
+		return res;
   } // krfact
 
-  void krfact(const Matrix& Mat, int colnum, Matrix& cols, Matrix& rows)
+  RowVector krfact(const Matrix& Mat, int colnum, Matrix& cols, Matrix& rows)
   {
+		RowVector res;
     cols = zeros(colnum,Mat.Ncols());
     rows = zeros(int(Mat.Nrows() / colnum),Mat.Ncols());
-    krfact(Mat,cols,rows);
+    res = krfact(Mat,cols,rows);
+		return res;
   } // krfact
 
   Matrix krprod(const Matrix& cols, const Matrix& rows)
@@ -558,35 +565,46 @@ namespace Melodic{
     int res = 0;
     ColumnVector PPCA;
  
+		if(which == string("aut"))
+			if(int(estimators(2)) < int(estimators(1)) && int(estimators(2)) > 15){
+				res=int(estimators(2));
+	      PPCA << PPCA2.Column(3);
+			}else{
+				res = int(estimators(1));
+	      PPCA << PPCA2.Column(2);
+			}				
     if(which == string("lap")){
       res = int(estimators(1));
       PPCA << PPCA2.Column(2);
     }
-
     if(which == string("bic")){
       res = int(estimators(2));
-      PPCA << PPCA2.Column(2);
+      PPCA << PPCA2.Column(3);
     }
     if(which == string("mdl")){
       res = int(estimators(3));
       PPCA << PPCA2.Column(4);
     }
+    if(which == string("rrn")){
+      res = int(estimators(4));
+      PPCA << PPCA2.Column(5);
+    }
     if(which == string("aic")){
       res = int(estimators(5));
       PPCA << PPCA2.Column(6);
     }
-    if(res==0){//median estimator
-      PPCA = PPCA2.Column(2);
-
-      for(int ctr=1; ctr<=PPCA2.Nrows(); ctr++){ 
-	RowVector tmp = PPCA2.SubMatrix(ctr,ctr,2,6);
-	PPCA(ctr) = float(tmp.Sum()/5);
-      }
-
-      ctr_i = 1; 
-      while((PPCA(ctr_i) < PPCA(ctr_i+1))&&(ctr_i<maxEV)){
-	res=ctr_i+1;ctr_i++;
-      }
+		if(which == string("median")){
+			RowVector unsorted(estimators);	
+			SortAscending(unsorted);
+			ctr_i=1;
+			res=int(unsorted(3));
+			while(res != int(estimators(ctr_i)))			
+				ctr_i++;
+			PPCA << PPCA2.Column(ctr_i);
+		}
+    if(res==0 || which == string("mean")){//median estimator
+      PPCA = mean(PPCA2.Columns(2,6),2);
+			res=int(mean(estimators,2).AsScalar());
     }
 
     dim = res;
@@ -812,10 +830,10 @@ namespace Melodic{
 			beta = pinvdes * dat;
 			residu = dat - design*beta;
 
-			dof = design.Nrows() - design.Ncols();
+			dof = design.Nrows() - design.Ncols()-1;
 			sigsq = sum(SP(residu,residu))/dof;
 			
-			float fact = float(design.Nrows() - 1 - design.Ncols()) / design.Ncols();
+			float fact = float(dof) / design.Ncols();
 			f_fmf =  SP(sum(SP(design*beta,design*beta)),pow(sum(SP(residu,residu)),-1)) * fact;
 		
 			pf_fmf = f_fmf.Row(1); 
diff --git a/melhlprfns.h b/melhlprfns.h
index 4d95a223a35f8b149e3af491f9458bf70186f253..4d9c7a533150b9679db8936062152630bcaf30c2 100644
--- a/melhlprfns.h
+++ b/melhlprfns.h
@@ -50,9 +50,9 @@ namespace Melodic{
   void em_pca(const Matrix& Mat, Matrix& evecs, RowVector& evals, int num_pc = 1, int iter = 20);
   void em_pca(const Matrix& Mat, Matrix& guess, Matrix& evecs, RowVector& evals, int num_pc = 1, int iter = 20);
 
-  void rankapprox(const Matrix& Mat, Matrix& cols, Matrix& rows, int dim = 1);
-  void krfact(const Matrix& Mat, Matrix& cols, Matrix& rows);
-  void krfact(const Matrix& Mat, int colnum, Matrix& cols, Matrix& rows);
+  float rankapprox(const Matrix& Mat, Matrix& cols, Matrix& rows, int dim = 1);
+  RowVector krfact(const Matrix& Mat, Matrix& cols, Matrix& rows);
+  RowVector krfact(const Matrix& Mat, int colnum, Matrix& cols, Matrix& rows);
   Matrix krprod(const Matrix& cols, const Matrix& rows);
   Matrix krapprox(const Matrix& Mat, int size_col, int dim = 1);
 
diff --git a/melodic.cc b/melodic.cc
index 4aeadd49ec12cc4711d1c99b1ef4bf3fd6619f08..d6d95d8431bb31bb58af5165b25ad2265e422b0d 100644
--- a/melodic.cc
+++ b/melodic.cc
@@ -71,13 +71,12 @@ int main(int argc, char *argv[]){
     if (opts.filtermode || opts.filtermix.value().length()>0){
       if(opts.filtermode){ // just filter out some noise from a previous run
     		melodat.setup();
+				if(opts.debug.value())
+					message("  Denoising data setup completed "<< endl);
 				melodat.remove_components();
-      }
-      else
+      } else
 				mmonly(logger,opts,melodat,report);
-    } 
-    else
-    {  // standard PICA now
+    } else {  // standard PICA now
       int retry = 0;
       bool no_conv;
       bool leaveloop = false;
@@ -104,12 +103,10 @@ int main(int argc, char *argv[]){
 	    			opts.approach.set_T("defl");
 	    			message(endl << "Restarting MELODIC using deflation approach" 
 		    			<< endl << endl);
-	  			}
-	  			else{
+	  			}else{
 	    			leaveloop = true;
 	  			}
-				}
-				else{
+				}else{
 	  			if(no_conv){
 	    			retry++;
 	    			if(opts.pca_dim.value()-retry*opts.retrystep > 
@@ -125,10 +122,10 @@ int main(int argc, char *argv[]){
 	      			}
 	    			}
 	    			if(!leaveloop){
-	      message(endl << "Restarting MELODIC using -d " 
-		      << opts.pca_dim.value() 
-		      << endl << endl);
-	    }
+	      			message(endl << "Restarting MELODIC using -d " 
+		      			<< opts.pca_dim.value() 
+		      			<< endl << endl);
+	          }
 	  			}
 				}
       } while (no_conv && retry<opts.maxRestart.value() && !leaveloop);	
diff --git a/meloptions.cc b/meloptions.cc
index 6d50d52b70524d8bbe542a551a02acb8487914e1..0b92f4d1000d2454b7a8aed3b340e4dc9739d8f4 100644
--- a/meloptions.cc
+++ b/meloptions.cc
@@ -32,7 +32,7 @@ MelodicOptions* MelodicOptions::gopt = NULL;
   		version = p_version;
   		filtermode = false;
   		explicitnums = false;
-  		logfname = string("melodic.log");
+  		logfname = string("log.txt");
 
   		// work out the path to the $FSLDIR/bin directory
   		if(getenv("FSLDIR")!=0){
@@ -112,8 +112,11 @@ MelodicOptions* MelodicOptions::gopt = NULL;
       		print_usage(argc,argv);
       		exit(2);
     		} else {   
+					temporal.set_T(false);
       		filtermode = true;
       		varnorm.set_T(false);
+					pbsc.set_T(false);
+					cerr << "WARNING: melodic denoising is deprecated, please use fsl_regfilt instead!" <<endl;
     		} 
   		}
   		if (threshold.value()<=0){
@@ -182,7 +185,7 @@ MelodicOptions* MelodicOptions::gopt = NULL;
   		inputfname.set_T(tmpfnames);
 
   		//create melodic directory name
-  		if(logdir.value()==string("melodic.log")){
+  		if(logdir.value()==string("log.txt")){
     		logdir.set_T(string(inputfname.value().at(0)+".ica"));
     		logger.makeDir(logdir.value(),logfname);
  			} else{
diff --git a/meloptions.h b/meloptions.h
index 9cbf13c2092738f1550b189624bf31c00ea6d831..15670a270167d379aca14c3365795284b387f230 100644
--- a/meloptions.h
+++ b/meloptions.h
@@ -81,6 +81,7 @@ class MelodicOptions {
 		Option<string> bgimage;
   	Option<float>  tr;
   	Option<bool>   logPower;
+		Option<bool>  addsigchng;
 
 		Option<string> fn_Tdesign;
 		Option<string> fn_Tcon;
@@ -106,8 +107,8 @@ class MelodicOptions {
   	Option<string> guessfname;
   	Option<string> paradigmfname;
 
-  	Option<int>  dummy;
-  	Option<int>  repeats;
+  	Option<int>   dummy;
+  	Option<int>   repeats;
   	Option<float> nlconst1;
   	Option<float> nlconst2;
   	Option<float> smooth_probmap;
@@ -147,7 +148,7 @@ class MelodicOptions {
  }
 
  inline MelodicOptions::MelodicOptions() :
-   logdir(string("-o,--outdir"), string("melodic.log"),
+   logdir(string("-o,--outdir"), string("log.txt"),
 	  string("output directory name\n"), 
 	  false, requires_argument),
    inputfname(string("-i,--in"), std::vector<string>(),
@@ -180,8 +181,8 @@ class MelodicOptions {
    joined_whiten(string("--sep_whiten"), true,
 	   string("switch on separate whitening"), 
 	   false, no_argument),
-	 joined_vn(string("--joined_vn"), false,
-		   string("switch on joined variance nomalisation"), 
+	 joined_vn(string("--sep_vn"), true,
+		   string("switch off joined variance nomalisation"), 
 		   false, no_argument, false),
    numICs(string("-n,--numICs"), -1,
 	   string("numer of IC's to extract (for deflation approach)"), 
@@ -200,7 +201,7 @@ class MelodicOptions {
 	   false, no_argument),
    pspec(string("--pspec"), false,
 	   string("        switch on conversion to powerspectra"), 
-	   false, no_argument),
+	   false, no_argument, false),
    segment(string("--covarweight"), string(""),
 	   string("voxel-wise weights for the covariance matrix (e.g. segmentation information)"),
 	   false, requires_argument),
@@ -249,9 +250,12 @@ class MelodicOptions {
    tr(string("--tr"),  0.0,
 	   string("\tTR in seconds"), 
 	   false, requires_argument),
-   logPower(string("--logPower"),  false,
-	   string("calculate log of power for frequency spectrum\n"), 
-	    false, no_argument),
+	 logPower(string("--logPower"),  false,
+		 string("calculate log of power for frequency spectrum\n"), 
+		 false, no_argument),
+	 addsigchng(string("--sigchng"),  false,
+		 string("add signal change estimates to report pages\n"), 
+		 false, no_argument, false),
 	 fn_Tdesign(string("--Tdes"), string(""),
 	   string("        design matrix across time-domain"),
 		 false, requires_argument),
@@ -382,6 +386,7 @@ class MelodicOptions {
 			options.add(bgimage);
 	    options.add(tr);
 	    options.add(logPower);
+	    options.add(addsigchng);
 			options.add(fn_Tdesign);
 			options.add(fn_Tcon);
 			options.add(fn_TconF);
diff --git a/melreport.cc b/melreport.cc
index 63e4491a0ea75c82e25a4d750129014005bc4e68..33c1eaa90aa6e3235f661b3a1af5053c332bdda6 100644
--- a/melreport.cc
+++ b/melreport.cc
@@ -53,7 +53,7 @@ namespace Melodic{
 	  			IChtml << fixed << setprecision(2) << ICstats(cnum,1) << " % of explained variance";
 	  			if(ICstats.Ncols()>1)
 	    			IChtml << "; &nbsp; &nbsp; " << ICstats(cnum,2) << " % of total variance";
-					if(ICstats.Ncols()>2&&melodat.get_numfiles()==1){
+					if(ICstats.Ncols()>2&&opts.addsigchng.value()){
 	    			IChtml << "<p>" <<endl;
 	    			IChtml << " &nbsp; &nbsp; " << ICstats(cnum,3) << " % signal change (pos peak voxel); &nbsp; &nbsp;" << ICstats(cnum,4) << "% signal change (peak neg. voxel)" << endl ;
 	  			}
@@ -112,16 +112,27 @@ namespace Melodic{
       {//plot time course
     	IChtml << "<H3> Temporal mode </H3><p>" << endl <<endl;
     	miscplot newplot;
-			Matrix tmptc = melodat.get_Tmodes(cnum-1).t();
+			Matrix tmptc = melodat.get_Tmodes(cnum-1).Column(1).t();
+			newplot.col_replace(0,0xFF0000);
 
+			newplot.add_label(string("IC ")+num2str(cnum)+" time course");
 			//add GLM OLS fit
 			if(melodat.Tdes.Storage()>0 &&
 		  melodat.glmT.get_beta().Nrows() == melodat.Tdes.Ncols()){
 				tmptc &= melodat.glmT.get_beta().Column(cnum).t() * melodat.Tdes.t();
-				newplot.add_label(string("IC ")+num2str(cnum)+" time course");
 				newplot.add_label("full model fit");
 			}
 
+			//add deviation around time course
+			if(melodat.get_Tmodes(cnum-1).Ncols()>1){
+				Matrix tmp = stdev(melodat.get_Tmodes(cnum-1).Columns(2,melodat.get_Tmodes(cnum-1).Ncols()).t(),1);
+				tmptc &= melodat.get_Tmodes(cnum-1).Column(1).t()+tmp;
+				tmptc &= melodat.get_Tmodes(cnum-1).Column(1).t()-tmp;
+				newplot.add_label("std error across subjects");
+				newplot.col_replace(tmptc.Nrows()-1,0x808080);
+				newplot.col_replace(tmptc.Nrows()-2,0x808080);				
+			}
+		
 	    if(opts.tr.value()>0.0)
 				newplot.add_xlabel(string("Time (seconds); TR = ")+
 				float2str(opts.tr.value(),0,2,0)+" s");
@@ -137,12 +148,19 @@ namespace Melodic{
 			  string("Timecourse No. ")+num2str(cnum), 
 			  opts.tr.value(),150,12,1,false);
 
+			if(melodat.get_Tmodes(cnum-1).Ncols()>1)
+				tmptc &= melodat.get_Tmodes(cnum-1).Columns(2,melodat.get_Tmodes(cnum-1).Ncols()).t();
 	     write_ascii_matrix(report.appendDir(string("t")
 				+num2str(cnum)+".txt"),tmptc.t());
 	     IChtml << "<A HREF=\"" << string("t")
 	  	  +num2str(cnum)+".txt" << "\"> ";
 			 IChtml << "<img BORDER=0 SRC=\"" 
 	  	  +string("t")+num2str(cnum)+".png\"></A><p>" << endl;	
+	
+			if(melodat.get_numfiles()>1 && melodat.explained_var.Storage()>0 
+				&& melodat.explained_var.Ncols()>=cnum)
+				IChtml << "Rank-1 approximation of individual time courses explains " 
+				<< melodat.explained_var(cnum) << "% of variance.<p>" << endl;
 			}//time series plot
 
 	 		if(!opts.pspec.value())
@@ -156,7 +174,7 @@ namespace Melodic{
 				else
 					newplot.add_ylabel(string("Power"));
 		
-				Matrix fmixtc = calc_FFT(melodat.get_Tmodes(cnum-1), opts.logPower.value());
+				Matrix fmixtc = calc_FFT(melodat.get_Tmodes(cnum-1).Column(1), opts.logPower.value());
 			  
 				newplot.set_Ylabel_fmt("%.0f");
 				newplot.set_yrange(0.0,1.02*fmixtc.Maximum());
@@ -239,11 +257,13 @@ namespace Melodic{
   
       if(cnum <= (int)melodat.get_Smodes().size())
 	    {//plot subject mode 
+			  
 	  		Matrix smode;
 	  		smode = melodat.get_Smodes(cnum-1);
 	
 	  		if(smode.Nrows() > 1){
-	      	miscplot newplot;
+	  	  	IChtml << "<hr><H3> Sessions/Subjects mode </H3><p>" << endl <<endl;
+		    	miscplot newplot;
 
 					//add GLM OLS fit
 					//newplot.setscatter(smode,2);
@@ -260,7 +280,7 @@ namespace Melodic{
 //					newplot.set_xysize(smode.Nrows()*80,150);
 	      	newplot.timeseries(smode.t(), 
 			    	report.appendDir(string("s")+num2str(cnum)+".png"),
-			      string("Subject/Session mode"));
+			      string("Subject/Session mode No. ") + num2str(cnum));
 					newplot.clear_xlabel();
 					newplot.clear_labels();
 					newplot.set_xysize(120,200);
@@ -722,7 +742,7 @@ namespace Melodic{
   void MelodicReport::PPCA_rep(){
     
     {//plot time course
-      report << "<p><hr><b>PCA estimates </b> <br>" << endl;
+      report << "<p><hr><b>PCA estimates </b> <p>" << endl;
  
       Matrix what;
       miscplot newplot;
@@ -753,12 +773,18 @@ namespace Melodic{
   }
 	
 	void MelodicReport::Smode_rep(){
-		report << "<b>TICA Subject/Session modes </b> <p>" << endl;
+		report << "<p><hr><b>TICA Subject/Session modes </b> <br>" << endl;
 		miscplot newplot;
+		report << "Boxplots show the relative response amplitudes across the "
+			<< "session/subject domain (" << melodat.get_numfiles() 
+			<< " input files). Components have been sorted in decreasing order of "
+			<< " the median response per component. <br><br>";
 		Matrix allmodes = melodat.get_Smodes().at(0);
 		for(int ctr = 1; ctr < (int)melodat.get_Smodes().size();++ctr)
 			allmodes |= melodat.get_Smodes().at(ctr);
 		
+		newplot.add_xlabel("Component No.");
+		newplot.add_ylabel("");
 		newplot.set_xysize(100+30*allmodes.Ncols(),300);
 		newplot.boxplot(allmodes,report.appendDir(string("bp_Smodes.png")),
     string("Subject/Session modes"));
diff --git a/melreport.h b/melreport.h
index e6ce996d57add175323886f03c4934341a5a80e4..8fbe5f18ed9cf40707810bddb1538a85704239f1 100644
--- a/melreport.h
+++ b/melreport.h
@@ -82,12 +82,12 @@ namespace Melodic{
 						}
 						report << "<IFRAME  height=80px width=100% src=nav.html frameborder=0></IFRAME><p>"<< endl;		
 						loghtml << "<IFRAME  height=100px width=100% src=nav.html frameborder=0></IFRAME><p>"
-							<<"<pre>../melodic.log</pre>" <<endl;		
+							<<"<IFRAME width=100% height=100% src=\"../log.txt\"></IFRAME>" <<endl;		
 						navigator <<"<CENTER><TABLE BORDER=0><TR>" << endl
 							<<"<TD ALIGN=CENTER WIDTH=100%><FONT SIZE=-1>"<<endl
 							<<"<A HREF=\"00index.html\" target=\"_top\">Main</A>&nbsp;-&nbsp;";
-			//			if(opts.guireport.value()=="")
-			//				navigator << "<A HREF=\"log.html\" target=\"_top\">Log</A>&nbsp;-&nbsp;";
+						if(opts.guireport.value()=="")
+							navigator << "<A HREF=\"log.html\" target=\"_top\">Log</A>&nbsp;-&nbsp;";
 						navigator	<<"Components: ";
 						navigator.flush();
 	  			}
@@ -113,9 +113,23 @@ namespace Melodic{
 				if( bool(opts.genreport.value()) ){
 
 	  			report << "<b>Analysis methods</b> <br>"<<endl;
-	  			report << "Analysis was carried out using MELODIC (Multivariate Exploratory Linear Decomposition into Independent Components) Version "<< version <<", part of FSL (FMRIB's Software Library, <A HREF=\"http://www.fmrib.ox.ac.uk/fsl/\">www.fmrib.ox.ac.uk/fsl</A>), an implementation for the estimation of a Probabilistic Independent Component Analysis model [Beckmann 2004]."<<endl;
-	  
-	  			report << "The following melodic pre-processing was applied to the input data file: "<< endl;
+					report << "Analysis was carried out using ";
+	
+					if(opts.approach.value() != string("tica"))
+						report << "Probabilistic Independent Component Analysis"
+						  <<" [Beckmann 2004] as implemented in "<<endl;
+					else	
+						report << "Tensorial Independent Component Analysis "
+					  	<<"[Beckmann 2005] as implemented in "<< endl;
+	
+					report << " MELODIC (Multivariate"
+						<<" Exploratory Linear Decomposition into Independent Components)"
+						<<" Version "<< version <<", part of FSL (FMRIB's Software"
+						<<" Library, <A HREF=\"http://www.fmrib.ox.ac.uk/fsl/\">"
+						<<"www.fmrib.ox.ac.uk/fsl</A>).<br>";
+			
+	  			report << "The following data pre-processing was applied to"
+						<<" the input data: "<< endl;
 
 	  			if(opts.use_mask.value())
 	    			report << " masking of non-brain voxels;";
@@ -124,35 +138,52 @@ namespace Melodic{
 	  
 	  			if(opts.varnorm.value())
 	    			report << " normalisation of the voxel-wise variance; ";
-	  
+	  	    if(opts.pbsc.value())
+						report << " conversion to %BOLD signal change; ";
 	  			report << "<br>"<<endl;
 	  
-	  			report << " Pre-processed data was whitened and projected into a " 
+	  			report << " Pre-processed data were whitened and projected into a " 
 		 				<< melodat.get_mix().Ncols()<< "-dimensional subspace using ";
 	  			if(melodat.get_PPCA().Storage()>0){	    
-	    			report << "probabilistic Principal Component Analysis where the number of dimensions was estimated using ";
+	    			report << "probabilistic Principal Component Analysis where the"
+							<<" number of dimensions was estimated using ";
 	    			if(opts.pca_est.value() == string("lap"))
-	      			report << "the Laplace approximation to the Bayesian evidence of the model order [Minka 2000, Beckmann 2004]. " << endl;
+	      			report << "the Laplace approximation to the Bayesian"
+							<<" evidence of the model order [Minka 2000, Beckmann 2004]. " << endl;
 	    			else
 	      			if(opts.pca_est.value() == string("bic"))
-							report << "the <em> Bayesian Information Criterion</em> (BIC) [Kass 1993]. " << endl;
+							  report << "the <em> Bayesian Information Criterion</em>"
+								  <<" (BIC) [Kass 1993]. " << endl;
 	      		else
 							if(opts.pca_est.value() == string("mdl"))
-		  				report << "<em> Minimum Description Length</em> (MDL) [Rissanen 1978]. " << endl;
+		  				  report << "<em> Minimum Description Length</em> (MDL)"
+						    	<<" [Rissanen 1978]. " << endl;
 						else
 		  				if(opts.pca_est.value() == string("aic"))
-		    			report << "the <em> Akaike Information Criterion</em> (AIC) [Akaike 1969]. " << endl;
+		    			report << "the <em> Akaike Information Criterion</em>"
+								<<" (AIC) [Akaike 1969]. " << endl;
 		  			else
-		    			report << "a committee of approximations to Bayesian the model order [Beckmann 2004]. " << endl;
+		    			report << " approximations to Bayesian the"
+								<<" model order [Beckmann 2004]. " << endl;
 	  			}	  
 	  			else
 	    			report << "Principal Component Analysis. ";
 	  
-	  			report << " The whitened observations were decomposed into a set of time-courses and spatial maps by optimising for non-Gaussian spatial source distributions using a fixed-point iteration technique [Hyv&auml;rinen 1999]. " << endl;	  
-	  			report << "Estimated Component maps were divided by the standard deviation of the residual noise";
+	  			report << " <BR>The whitened observations were decomposed into "
+						<<" sets of vectors which describe signal variation across"
+						<<" the temporal domain (time-courses)";
+					if(opts.approach.value() == string("tica") ||
+						opts.approach.value() == string("concat"))
+						report << ", the session/subject domain ";
+					report <<"  and across the spatial domain (maps) by optimising for"
+						<<" non-Gaussian spatial source distributions using a"
+						<<" fixed-point iteration technique [Hyv&auml;rinen 1999]. " << endl;	  
+	  			report << "Estimated Component maps were divided by the standard"
+						<<" deviation of the residual noise";
 	  
 	  			if(opts.perf_mm.value())
-	    			report << " and thresholded by fitting a mixture model to the histogram of intensity values [Beckmann 2004]. <p>" << endl;
+	    			report << " and thresholded by fitting a mixture model "
+	 						<<"to the histogram of intensity values [Beckmann 2004]. <p>" << endl;
 	  			else
 	    			report <<".<p>" << endl;
 	 
@@ -164,26 +195,56 @@ namespace Melodic{
 				if( bool(opts.genreport.value()) ){
 	  			report << "<b>References</b> <br>"<<endl;
 	  
-	  			report << "[Hyv&auml;rinen 1999] A. Hyv&auml;rinen. Fast and Robust Fixed-Point Algorithms for Independent Component Analysis. IEEE Transactions on Neural Networks 10(3):626-634, 1999.<br> " << endl;
- 					report << "[Beckmann 2004] C.F. Beckmann and S.M. Smith. Probabilistic Independent Component Analysis for Functional Magnetic Resonance Imaging. IEEE Transactions on Medical Imaging 23(2):137-152 2004. <br>" << endl;
+	  			report << "[Hyv&auml;rinen 1999] A. Hyv&auml;rinen. Fast and"
+						<<" Robust Fixed-Point Algorithms for Independent Component"
+						<<" Analysis. IEEE Transactions on Neural Networks 10(3):"
+						<<"626-634, 1999.<br> " << endl;
+ 						report << "[Beckmann 2004] C.F. Beckmann and S.M. Smith."
+							<<" Probabilistic Independent Component Analysis for Functional"
+							<<" Magnetic Resonance Imaging. IEEE Transactions on Medical"
+							<<" Imaging 23(2):137-152 2004. <br>" << endl;
+					if(opts.approach.value() == string("tica") || 
+						opts.approach.value() == string("concat") )
+						report << "[Beckmann 2005] C.F. Beckmann and S.M. Smith."
+							<<" Tensorial extensions of independent component analysis"
+							<< " for multisubject FMRI analysis. Neuroimage "
+							<< " 25(1):294-311 2005. <br>";
 
 	  			if(melodat.get_PPCA().Storage()>0){	    
-	    			report << "[Everson 2000] R. Everson and S. Roberts. Inferring the eigenvalues of covariance matrices from limited, noisy data. IEEE Trans Signal Processing, 48(7):2083-2091, 2000<br>"<<endl;
-	    			report << "[Tipping 1999] M.E. Tipping and C.M.Bishop. Probabilistic Principal component analysis. J Royal Statistical Society B, 61(3), 1999. <br>" << endl;
-	    		/*  report << "[Beckmann 2001] C.F. Beckmann, J.A. Noble and S.M. Smith. Investigating the intrinsic dimensionality of FMRI data for ICA. In Seventh Int. Conf. on Functional Mapping of the Human Brain, 2001. <br>" << endl;*/
+	    			report << "[Everson 2000] R. Everson and S. Roberts."
+							<<" Inferring the eigenvalues of covariance matrices from"
+							<<" limited, noisy data. IEEE Trans Signal Processing,"
+							<<" 48(7):2083-2091, 2000<br>"<<endl;
+	    			report << "[Tipping 1999] M.E. Tipping and C.M.Bishop."
+							<<" Probabilistic Principal component analysis. J Royal"
+							<<" Statistical Society B, 61(3), 1999. <br>" << endl;
+	    		  report << "[Beckmann 2001] C.F. Beckmann, J.A. Noble and"
+							<<" S.M. Smith. Investigating the intrinsic dimensionality"
+							<<" of FMRI data for ICA. In Seventh Int. Conf. on Functional"
+							<<" Mapping of the Human Brain, 2001. <br>" << endl;
 	    			if(opts.pca_est.value() == string("lap"))
-	      			report << "[Minka 2000] T. Minka. Automatic choice of dimensionality for PCA. Technical Report 514, MIT Media Lab Vision and Modeling Group, 2000. <BR>"<< endl;
+	      			report << "[Minka 2000] T. Minka. Automatic choice of"
+							  <<" dimensionality for PCA. Technical Report 514, MIT"
+								<<" Media Lab Vision and Modeling Group, 2000. <BR>"<< endl;
 	    			else
 	      			if(opts.pca_est.value() == string("bic"))
-							report << "[Kass 1995] R.E. Kass and A. E. Raftery. Bayes factors. Journal of the American Statistical Association, 90:733-795, 1995 <br>" << endl;
+							report << "[Kass 1995] R.E. Kass and A. E. Raftery. Bayes"
+								<<" factors. Journal of the American Statistical"
+								<<" Association, 90:733-795, 1995 <br>" << endl;
 	      		else
 							if(opts.pca_est.value() == string("mdl"))
-		  				report << "[Rissanen 1978]. J. Rissanen. Modelling by shortest data description. Automatica, 14:465-471, 1978. <br>" << endl;
+		  				report << "[Rissanen 1978]. J. Rissanen. Modelling by"
+							  <<" shortest data description. Automatica,"
+								<<" 14:465-471, 1978. <br>" << endl;
 						else
 		  				if(opts.pca_est.value() == string("aic"))
-		    			report << "[Akaike 1974]. H. Akaike. A new look at statistical model identification. IEEE Transactions on Automatic Control, 19:716-723, 1974. <br>" << endl;
+		    			report << "[Akaike 1974]. H. Akaike. A new look at"
+							  <<" statistical model identification. IEEE Transactions"
+								<<" on Automatic Control, 19:716-723, 1974. <br>" << endl;
 		  			else
-		    			report << "[Minka 2000]. T. Minka. Automatic choice of dimensionality for PCA. Technical Report 514, MIT Media Lab Vision and Modeling Group, 2000. <BR>" << endl;
+		    			report << "[Minka 2000]. T. Minka. Automatic choice of"
+								<<" dimensionality for PCA. Technical Report 514, MIT"
+								<<" Media Lab Vision and Modeling Group, 2000. <BR>" << endl;
 
 					}
 				}