#include <iostream.h>
#include <iomanip>
#include "libvis/miscplot.h"
#include "libvis/miscpic.h"
#include <iostream>
#include "newmatap.h"
#include "newmatio.h"
#include "newimage/newimageall.h"
#include "miscmaths/miscmaths.h"
#include "miscmaths/miscprob.h"
#include <string>
#include "utils/log.h"
using namespace std;
using namespace Utilities;
using namespace NEWIMAGE;
using namespace MISCPLOT;
using namespace MISCPIC;
Matrix calcFFT(const Matrix& Mat)
{
Matrix res;
for(int ctr=1; ctr <= Mat.Ncols(); ctr++)
{
ColumnVector tmpCol;
tmpCol=Mat.Column(ctr);
ColumnVector FtmpCol_real;
ColumnVector FtmpCol_imag;
ColumnVector tmpPow;
if(tmpCol.Nrows()%2 != 0){
Matrix empty(1,1); empty=0;
tmpCol &= empty;}
RealFFT(tmpCol,FtmpCol_real,FtmpCol_imag);
tmpPow = pow(FtmpCol_real,2)+pow(FtmpCol_imag,2);
tmpPow = tmpPow.Rows(2,tmpPow.Nrows());
//if(opts.logPower.value()) tmpPow = log(tmpPow);
if(res.Storage()==0){res= tmpPow;}else{res|=tmpPow;}
}
return res;
} //Matrix calc_FFT()
string float2str(float f, int width, int prec, int scientif)
{
ostrstream os;
int redw = int(std::abs(std::log10(std::abs(f))))+1;
if(width>0)
os.width(width);
if(scientif>0)
os.setf(ios::scientific);
os.precision(redw+std::abs(prec));
os.setf(ios::internal, ios::adjustfield);
os << f << '\0';
return os.str();
}
int main(int argc, char *argv[])
{
Matrix data;
cerr<< "Reading mixing matrix ";
data = read_ascii_matrix(string("mixmat"));
if (data.Storage()<=0) {
cerr <<" Please specify the mixing matrix correctly" << endl;
}
cerr << "done " << endl;
cerr << "Matrix size : " << data.Nrows() << " x " << data.Ncols() << endl;
//set up data
if(data.Ncols()==0){
cerr << " create new data " << M_PI <<endl;
ColumnVector X(200);
for(int i=1; i<=X.Nrows(); i++){
X(i) = (float((i-1))/(X.Nrows()-1))*M_PI;
}
RowVector Y(10);
for(int i=1; i<=Y.Ncols(); i++)
Y(i)=3*i;
data = X*Y;
for(int i=1; i<=data.Ncols(); i++)
for(int j=1; j<=data.Nrows(); j++)
data(j,i) = sin(data(j,i));
cerr << "Matrix size : " << data.Nrows() << " x " << data.Ncols() << endl;
Log mmat;
mmat.setDir(string("./res/"),string("mix.txt"));
mmat << data <<endl;
}
Log IChtml;
IChtml.setDir(string("./res/"),string("index.html"));
IChtml << "<HTML> " << endl
<< "<TITLE>MELODIC Component "
<< "</TITLE>" << endl
<< "<BODY BACKGROUND=\"file:" << getenv("FSLDIR")
<< "/doc/images/fsl-bg.jpg\">" << endl
<< "<hr><CENTER><H1>MELODIC Component "
<< "</H1>"<< endl;
float tr = 0.0;
Matrix data2;
data2 = calcFFT(data);
for (int i=1; i<=data.Ncols(); i++)
{
cerr << i << " " ;
{//plot time course
miscplot newplot;
if(tr>0.0)
newplot.timeseries(data.Column(i).t(),
string("./res/t")+num2str(i)+".png",
string("Timecourse (in seconds); TR = ")+
float2str(tr,0,2,0)+" s",
tr);
else
newplot.timeseries(data.Column(i).t(),
string("./res/t")+num2str(i)+".png",
string("Timecourse (in TRs)"));
// write_ascii_matrix(string("./res/t")+num2str(i)+".txt"),
// data.Column(i));
Log tfile;
tfile.setDir(string("./res/"),string("t")+num2str(i)+".txt");
tfile << data.Column(i) << endl;
IChtml << "<A HREF=\"" << string("t")
+num2str(i)+".txt" << "\"> ";
IChtml << "<img BORDER=0 SRC=\""
+string("t")+num2str(i)+".png\"></A><p>" << endl;
}//time series plot
{//plot frequency
miscplot newplot;
int fact = int(std::pow(10.0,int(std::log10(float(data.Nrows())))));
cerr << " FACT : "<<fact;
if(tr>0.0)
newplot.timeseries(data2.Column(i).t(),
string("./res/f")+num2str(i)+".png",
string("FFT of timecourse (in Hz / ") +num2str(fact)+")",
fact/(tr*data2.Nrows()));
else
newplot.timeseries(data2.Column(i).t(),
string("./res/f")+num2str(i)+".png",
string(string("FFT of timecourse (in cycles); ")
+"frequency(Hz)=cycles/("
+num2str(data2.Nrows())
+"* TR); period(s)=("
+num2str(data2.Nrows())
+"* TR)/cycles"));
Log ffile;
ffile.setDir(string("./res/"),string("f")+num2str(i)+".txt");
ffile << data2.Column(i) << endl;
IChtml << "<A HREF=\"" << string("f")
+num2str(i)+".txt" << "\"> ";
IChtml << "<img BORDER=0 SRC=\""
+string("f")+num2str(i)+".png\"></A><p>" << endl;
}//frequency plot
}
IChtml<< "<HR><FONT SIZE=1>This page produced automatically by "
<< "FMRIB Software Library</A>.</FONT>" << endl
<< "</BODY></HTML>" << endl;
cerr << endl;
return 0;
}