Commit ebdaa2fb authored by thomas.forbriger's avatar thomas.forbriger Committed by thomas.forbriger
Browse files

proceeding - wow

This is a legacy commit from before 2015-03-01.
It may be incomplete as well as inconsistent.
See COPYING.legacy and README.history for details.


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/trunk
SVN Revision: 1493
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent fcd2c942
......@@ -3,7 +3,7 @@
*
* ----------------------------------------------------------------------------
*
* $Id: sigfit.cc,v 1.9 2004-02-15 20:43:09 tforb Exp $
* $Id: sigfit.cc,v 1.10 2004-02-24 16:02:47 tforb Exp $
* \author Thomas Forbriger
* \date 28/01/2004
*
......@@ -13,13 +13,15 @@
*
* REVISIONS and CHANGES
* - 28/01/2004 V1.0 Thomas Forbriger
* - 24/02/2004 V1.1 introduced search range stabilization
* more output values
*
* ============================================================================
*/
#define SIGFIT_VERSION \
"SIGFIT V1.0 fit signal by trial-signals"
"SIGFIT V1.1 fit signal by trial-signals"
#define SIGFIT_CVSID \
"$Id: sigfit.cc,v 1.9 2004-02-15 20:43:09 tforb Exp $"
"$Id: sigfit.cc,v 1.10 2004-02-24 16:02:47 tforb Exp $"
#include <fstream>
#include <iostream>
......@@ -53,6 +55,8 @@ struct Options {
bool truncate;
std::string residualname;
bool writeresidual;
double searchrange;
bool usesearchrange;
}; // struct Options
// formatted number output
......@@ -71,6 +75,7 @@ int main(int iargc, char* argv[])
{
SIGFIT_VERSION "\n"
"usage: sigfit [-v] [-Tdate v] [-truncate]" "\n"
" [-searchrange[=r]]" "\n"
" signal trial [trial ...]" "\n"
" or: sigfit --help|-h" "\n"
};
......@@ -87,6 +92,11 @@ int main(int iargc, char* argv[])
" v give the tolerance in units of the sampling interval" "\n"
"-truncate truncate all series to the same number of samples" "\n"
"-residual f write waveforms containing residual to file \"f\"" "\n"
"-searchrange[=r] set search range to stabilize the fit in cases" "\n"
" where a trial series does not (or almost not)" "\n"
" contribute to the signal" "\n"
" the search range r is given in units of" "\n"
" rms(signal)/rms(trial signals)" "\n"
"\n"
"signal signal to by explained by a linear combination" "\n"
" of trial signals" "\n"
......@@ -107,6 +117,8 @@ int main(int iargc, char* argv[])
{"truncate",arg_no,"-"},
// 4: output SFF file
{"residual",arg_yes,"-"},
// 5: search range
{"searchrange",arg_opt,"100."},
{NULL}
};
......@@ -159,6 +171,8 @@ int main(int iargc, char* argv[])
opt.truncate=cmdline.optset(3);
opt.writeresidual=cmdline.optset(4);
opt.residualname=cmdline.string_arg(4);
opt.usesearchrange=cmdline.optset(5);
opt.searchrange=cmdline.double_arg(5);
/*----------------------------------------------------------------------*/
......@@ -274,12 +288,13 @@ int main(int iargc, char* argv[])
}
/*----------------------------------------------------------------------*/
double signalrms=ts::rms(signal);
// set up system of linear equations
if (opt.verbose) { cout << "set up system of linear equations..." << endl; }
int N=bundlevec.size();
if (opt.verbose) { cout << "system is of size " << N << "x" << N << endl; }
Tmatrix Matrix(N,N), rhs(N), coeff(N);
Tmatrix Matrix(N,N), rhs(N), coeff(N), normproduct(N), trialrms(N);
for (int i=1; i<=N; ++i)
{
for (int k=i; k<=N; ++k)
......@@ -288,6 +303,29 @@ int main(int iargc, char* argv[])
Matrix(k,i)=Matrix(i,k);
}
rhs(i)=ts::innerproduct(bundlevec[i-1], signal);
trialrms(i)=ts::rms(bundlevec[i-1]);
normproduct(i)=rhs(i)/(signalrms*trialrms(i)*signal.size());
}
// add stabilization
if (opt.usesearchrange)
{
if (opt.verbose)
{
cout << "stabilize fit:" << endl;
cout << " relative search range: " << opt.searchrange << endl;
}
double sreachfac=signal.size()*signalrms;
for (int i=1; i<=N; ++i)
{
double range=opt.searchrange*signalrms/trialrms(i);
if (opt.verbose)
{
cout << " search range for trial series " << i
<< " is " << range << endl;
}
Matrix(i,i)+=pow((range/sreachfac),2);
}
}
// solve
......@@ -337,9 +375,17 @@ int main(int iargc, char* argv[])
}
// calculate rms values
double signalrms=ts::rms(signal);
double residualrms=ts::rms(residual);
double explained=1.-(residualrms/signalrms);
// calculate direction cosines
Tmatrix normcoeff(N);
double length=0.;
for (int i=1; i<=N; ++i)
{ length += coeff(i)*coeff(i); }
length=sqrt(length);
for (int i=1; i<=N; ++i)
{ normcoeff(i) = coeff(i)/length; }
std::string timestring(signal.header.date.timestring());
cout << " time: " << timestring.substr(4,21) << endl;
......@@ -353,10 +399,21 @@ int main(int iargc, char* argv[])
for (int i=1; i<=N; ++i)
{ cout << coeff(i) << " "; }
cout << endl;
cout << " rms coeffic.: " << length << endl;
cout << " normalized coefficients: ";
for (int i=1; i<=N; ++i)
{ cout << normcoeff(i) << " "; }
cout << endl;
cout << " coefficients x rms: ";
for (int i=1; i<=N; ++i)
{
cout << coeff(i)*ts::rms(bundlevec[i-1]) << " ";
cout << coeff(i)*trialrms(i) << " ";
}
cout << endl;
cout << " normalized scalar product: ";
for (int i=1; i<=N; ++i)
{
cout << normproduct(i) << " ";
}
cout << endl;
......@@ -372,9 +429,18 @@ int main(int iargc, char* argv[])
{
cout << formatfloat(coeff(i));
}
cout << formatfloat(length);
for (int i=1; i<=N; ++i)
{
cout << formatfloat(normcoeff(i));
}
for (int i=1; i<=N; ++i)
{
cout << formatfloat(coeff(i)*trialrms(i));
}
for (int i=1; i<=N; ++i)
{
cout << formatfloat(coeff(i)*ts::rms(bundlevec[i-1]));
cout << formatfloat(normproduct(i));
}
cout << endl;
}
......
Markdown is supported
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