Skip to content
Snippets Groups Projects
Commit 4bb15891 authored by Christian Beckmann's avatar Christian Beckmann
Browse files

new test prog

parent 6855896f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 180 additions and 7 deletions
test.cc
+ 180
− 7
View file @ 4bb15891
...@@ -2,18 +2,191 @@ ...@@ -2,18 +2,191 @@
#include <iomanip> #include <iomanip>
#include "libvis/miscplot.h" #include "libvis/miscplot.h"
#include "libvis/miscpic.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 std;
using namespace Utilities;
using namespace NEWIMAGE;
using namespace MISCPLOT;
using namespace MISCPIC;
template<class T>
int foo(const T& bar){
return 1; 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[]) int main(int argc, char *argv[])
{ {
cerr << foo(int(5));
return 1; 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;
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment