#include <DampingResults.h>
Classes | |
| struct | Result |
Public Slots | |
| void | save (QString fileName=QString::null) |
Public Member Functions | |
| bool | compute (int ig, const Signal *sig, const Parameters ¶m) |
| void | createObjects (SubSignalPool *subPool) |
| DampingResults (QWidget *parent=0) | |
| void | setWindowTitle (QString title) |
| ~DampingResults () | |
Public Attributes | |
| QMenu * | menuTools |
Static Protected Member Functions | |
| static double | misfitFunk (double x[]) |
Protected Attributes | |
| LineLayer ** | _averageCurveLayers |
| TextEdit ** | _comments |
| LineLayer ** | _fitCurveLayers |
| QString | _log |
| QVector< Result > | _results |
| LineLayer ** | _stddevCurveLayers |
Static Protected Attributes | |
| static Curve< Point2D > | _aveCurve |
| static double | _fitW = 0 |
| static int | _nSampFit = 0 |
| DampingResults::DampingResults | ( | QWidget * | parent = 0 | ) |
References SciFigs::GraphicSheetMenu::addMenu(), menuTools, save(), QGpCoreTools::tr(), and TRACE.
:
GraphicSheetMenu(parent)
{
TRACE;
Settings::getSize(this, "DampingResults" );
QAction * a;
menuTools=addMenu(tr( "&Tools" ));
a=new QAction(tr( "&Save results" ), this);
a->setStatusTip(tr( "Save results of spectral analysis in a text file with process information (.log file)" ));
connect(a, SIGNAL(triggered()), this, SLOT(save()) );
menuTools->addAction(a);
}
References _averageCurveLayers, _comments, _fitCurveLayers, _stddevCurveLayers, and TRACE.
{
TRACE;
Settings::setSize(this, "DampingResults" );
delete [] _comments;
delete [] _averageCurveLayers;
delete [] _stddevCurveLayers;
delete [] _fitCurveLayers;
}
| bool DampingResults::compute | ( | int | ig, |
| const Signal * | sig, | ||
| const Parameters & | param | ||
| ) |
References _aveCurve, _averageCurveLayers, _comments, _fitCurveLayers, _fitW, _log, _nSampFit, _results, _stddevCurveLayers, GeopsyCore::TimeRangeParameters::absoluteRange(), SciFigs::LineLayer::addLine(), QGpCoreTools::amoeba(), QGpCoreTools::Curve< pointType >::at(), SciFigs::GraphContent::boundingRect(), SciFigs::LineLayer::clear(), QGpCoreTools::Curve< pointType >::clear(), CONST_LOCK_SAMPLES, GeopsyCore::DoubleSignal::copySamplesFrom(), SciFigs::AxisWindow::deepUpdate(), GeopsyCore::DoubleSignal::deltaT(), GeopsyCore::TimeRange::end(), QGpCoreTools::endl(), QGpCoreTools::exp(), Parameters::filter(), GeopsyCore::DoubleSignal::filter(), Parameters::fitLength(), SciFigs::GraphContent::graph(), SciFigs::GraphContentLayer::graphContent(), GeopsyCore::SignalTemplate< sampleType >::initValues(), Parameters::isFilter(), GeopsyCore::TimeRange::lengthSeconds(), misfitFunk(), MSG_ID, GeopsyCore::Signal::nameComponent(), QGpCoreTools::Curve< pointType >::resize(), GeopsyCore::DoubleSignal::restoreType(), GeopsyCore::DoubleSignal::saveType(), GeopsyCore::DoubleSignal::setDeltaT(), SciFigs::Axis::setRange(), SciFigs::TextEdit::setText(), QGpCoreTools::sin(), QGpCoreTools::sqrt(), GeopsyCore::TimeRange::start(), GeopsyCore::DoubleSignal::subtractValue(), GeopsyCore::Signal::t0(), Parameters::timeRange(), Parameters::toString(), QGpCoreTools::tr(), TRACE, UNLOCK_SAMPLES, SciFigs::TextEdit::update(), w, Parameters::windowLength(), QGpCoreTools::Point2D::x(), QGpCoreTools::Rect::x1(), QGpCoreTools::Rect::x2(), SciFigs::AxisWindow::xAxis(), QGpCoreTools::Point2D::y(), QGpCoreTools::Rect::y1(), QGpCoreTools::Rect::y2(), and SciFigs::AxisWindow::yAxis().
{
TRACE;
_log=param.toString();
_averageCurveLayers[ig]->clear();
_stddevCurveLayers[ig]->clear();
_fitCurveLayers[ig]->clear();
// Cut the time window to process
TimeRange tw=param.timeRange().absoluteRange(sig);
double dt=sig->deltaT();
int nSampSig=(int)round((tw.end()-tw.start())/dt);
int nSampWin=(int)ceil(param.windowLength()/dt);
_nSampFit=(int)ceil(param.fitLength()/dt);
if(nSampWin>nSampSig) nSampWin=nSampSig;
if(_nSampFit>nSampWin) {
Message::warning(MSG_ID, tr("Damping toolbox"),
tr("The fitting length cannot be greater than the window length. "
"Continuing with a fitting length reduced to window length"), Message::ok(), true);
_nSampFit=nSampWin;
}
// Get result vectors ready
_aveCurve.clear();
Curve<Point2D> minCurve, maxCurve, fitCurve;
_aveCurve.resize(nSampWin);
minCurve.resize(nSampWin);
maxCurve.resize(nSampWin);
fitCurve.resize(_nSampFit);
// Set x coordinates for curves, and init y
double t=0;
int iw;
for(iw=0;iw < _nSampFit;iw++, t += dt) {
_aveCurve[iw].setX(t);
_aveCurve[iw].setY(0);
minCurve[iw].setX(t);
minCurve[iw].setY(0);
maxCurve[iw].setX(t);
maxCurve[iw].setY(0);
fitCurve[iw].setX(t);
fitCurve[iw].setY(0);
}
for( ;iw < nSampWin;iw++, t += dt) {
_aveCurve[iw].setX(t);
_aveCurve[iw].setY(0);
minCurve[iw].setX(t);
minCurve[iw].setY(0);
maxCurve[iw].setX(t);
maxCurve[iw].setY(0);
}
int nWin=0;
double tau, val; // time shift relative to signal sampling
const DoubleSignal * sigfft=sig->saveType(DoubleSignal::Waveform);
DoubleSignal * tSig=new DoubleSignal(nSampSig);
tSig->initValues(0.0, 0, nSampSig);
tSig->setDeltaT(dt);
tSig->copySamplesFrom(sigfft, tw.start()-sig->t0(), 0.0, tw.lengthSeconds());
tSig->subtractValue();
sig->restoreType(sigfft);
if(param.isFilter()) {
tSig->filter(param.filter());
}
CONST_LOCK_SAMPLES(double, samp, tSig)
// Stack windows (sum for ave, sum2 for min)
int ns=nSampSig - nSampWin - 1;
for(int is=0;is < ns;is++ ) {
if(samp[ is ] <= 0 && samp[ is + 1 ] >= 0) {
nWin++;
if(samp[ is ]==0)
tau=0;
else
tau=samp[ is ]/(samp[ is ] - samp[ is + 1 ] );
for(iw=1;iw < nSampWin;iw++ ) {
val=samp[ is + iw ] + tau * (samp[ is + iw + 1 ] - samp[ is + iw ] );
_aveCurve[iw].setY(_aveCurve.at(iw).y() + val);
minCurve[iw].setY(minCurve.at(iw).y() + val * val);
}
}
}
UNLOCK_SAMPLES(tSig);
delete tSig;
// Compute average and std deviation
for(iw=1;iw < nSampWin;iw++ ) {
val=_aveCurve.at(iw).y()/(double) nWin;
_aveCurve[iw].setY(val);
tau=sqrt(minCurve[iw].y()/(double) nWin - val * val);
minCurve[iw].setY(val - tau);
maxCurve[iw].setY(val + tau);
}
// Fit function A/w*e^(-Zwt)*sin(wt), A, w, Z are the unknowns
// Use the downhill simplex algorithm to find the minimum
// Init the algorithm with good estimates. Estimate the amplitude of signal
// average at the first absolute maximum
double z=0.01, a=0;
int iw0=0, nw0=0;
int nw=_nSampFit - 1;
for(iw=1;iw < nw;iw++ ) {
if(_aveCurve[iw].y() < 0 && _aveCurve[iw + 1].y() > 0) {
nw0++;
iw0=iw;
}
}
if(iw0==0) {
Message::warning(MSG_ID, tr("Damping"), tr("For signal %1, the first maximum cannot be found. This is "
"likely that the specified time range contains no signal at all. Check the time "
"limits.").arg(sig->nameComponent()), Message::cancel());
return false;
}
_fitW=2 * nw0 * M_PI/(iw0 * dt);
// a=aveCurve->at(1).y/dt;
for(iw=1;iw < _nSampFit;iw++ ) {
if(fabs( _aveCurve.at(iw).y()) > a) {
a=fabs(_aveCurve.at(iw).y());
iw0=iw;
}
}
a *= _fitW/exp( -z * _fitW * iw0 * dt);
// Init an initial tetrahedron
double * p[ 3 ], pv[ 6 ], y[ 3 ];
p[ 0 ]=pv;
p[ 1 ]=pv + 2;
p[ 2 ]=pv + 4;
// p[0][0]=a;
// p[0][1]=w;
// p[0][2]=z;
// y[0]=misfitFunk(p[0]);
// p[1][0]=a*1.01;
// p[1][1]=w*1.1;
// p[1][2]=z*1.5;
// y[1]=misfitFunk(p[1]);
// p[2][0]=a*0.99;
// p[2][1]=w*0.9;
// p[2][2]=z*0.5;
// y[2]=misfitFunk(p[2]);
// p[3][0]=a*1.005;
// p[3][1]=w*1.05;
// p[3][2]=z*1.2;
// y[3]=misfitFunk(p[3]);
// p[0][0]=a*0.9;
// p[0][1]=w;
// p[0][2]=z*0.5;
// y[0]=misfitFunk(p[0]);
// p[1][0]=a*1.1;
// p[1][1]=p[0][1];
// p[1][2]=p[0][2];
// y[1]=misfitFunk(p[1]);
// p[2][0]=a;
// p[2][1]=w;
// p[2][2]=p[0][2];
// y[2]=misfitFunk(p[2]);
// p[3][0]=a;
// p[3][1]=w;
// p[3][2]=z*1.5;
// y[3]=misfitFunk(p[3]);
p[ 0 ][ 0 ]=a * 0.9;
p[ 0 ][ 1 ]=z * 0.9;
y[ 0 ]=misfitFunk(p[ 0 ] );
p[ 1 ][ 0 ]=a * 0.9;
p[ 1 ][ 1 ]=z * 1.1;
y[ 1 ]=misfitFunk(p[ 1 ] );
p[ 2 ][ 0 ]=a * 1.1;
p[ 2 ][ 1 ]=z;
y[ 2 ]=misfitFunk(p[ 2 ] );
int nfunk;
amoeba(p, y, 2, 2.5e-8, misfitFunk, &nfunk);
App::stream() << tr( "%1 function evaluations were necessary to converge").arg(nfunk) << endl;
// Get the mean a,w and z
a=0.25 * (p[ 0 ][ 0 ] + p[ 1 ][ 0 ] + p[ 2 ][ 0 ] );
z=0.25 * (p[ 0 ][ 1 ] + p[ 1 ][ 1 ] + p[ 2 ][ 1 ] );
// get the fitCurve
for(iw=1;iw < _nSampFit;iw++ ) {
t=fitCurve.at(iw).x();
fitCurve[iw].setY(a/_fitW * exp( -z * _fitW * t) * sin(_fitW * t) );
}
// Set comments
_comments[ ig ] ->setText(tr("f=%1 Hz\nz=%2%\nN win=%3")
.arg(_fitW/(2 * M_PI), 0, 'f', 2)
.arg(z * 100.0, 0, 'f', 4)
.arg(nWin));
_comments[ ig ] ->update();
// Add curves to graph
static_cast<PlotLine2D *>(_averageCurveLayers[ig]->addLine())->setCurve(_aveCurve);
static_cast<PlotLine2D *>(_stddevCurveLayers[ig]->addLine())->setCurve(minCurve);
static_cast<PlotLine2D *>(_stddevCurveLayers[ig]->addLine())->setCurve(maxCurve);
static_cast<PlotLine2D *>(_fitCurveLayers[ig]->addLine())->setCurve(fitCurve);
// Set limits correctly
GraphContent * gc=_averageCurveLayers[ig]->graphContent();
AxisWindow * w=gc->graph();
Rect r=gc->boundingRect();
w->xAxis()->setRange(r.x1(), r.x2());
double maxY=-r.y1();
if(r.y2()>maxY) maxY=r.y2();
maxY*=1.05;
w->yAxis()->setRange(-maxY, maxY);
w->deepUpdate();
// Keep results
_results[ig].frequency=_fitW/(2*M_PI);
_results[ig].damping=z*100.0;
_results[ig].nWindows=nWin;
return true;
}
| void DampingResults::createObjects | ( | SubSignalPool * | subPool | ) |
References _averageCurveLayers, _comments, _fitCurveLayers, _results, _stddevCurveLayers, SciFigs::GraphicSheetMenu::addGraph(), SciFigs::GraphicSheetMenu::addText(), GeopsyCore::SubSignalPool::at(), GeopsyCore::SubSignalPool::count(), geopsyGui, GeopsyCore::Signal::nameComponent(), SciFigs::GraphicSheetMenu::setGraphGeometry(), SciFigs::Axis::setNumberPrecision(), SciFigs::Axis::setNumberPrecisionInversedScale(), SciFigs::Axis::setNumberType(), SciFigs::AbstractLine::setPen(), SciFigs::LineLayer::setReferenceLine(), SciFigs::GraphicSheet::setStatusBar(), SciFigs::Axis::setTitle(), SciFigs::Axis::setTitleInversedScale(), SciFigs::GraphicSheetMenu::sheet(), TRACE, w, SciFigs::AxisWindow::xAxis(), and SciFigs::AxisWindow::yAxis().
{
TRACE;
sheet()->setStatusBar(geopsyGui->statusBar());
double y=0.5;
// Resize containers to fit the number of signals
int n=subPool->count();
_comments=new TextEdit* [n];
_averageCurveLayers=new LineLayer* [n];
_stddevCurveLayers=new LineLayer* [n];
_fitCurveLayers=new LineLayer* [n];
_results.resize(n);
for(int ig=0; ig<n; ig++) {
AxisWindow * w=addGraph();
setGraphGeometry(w, 0.5, 15.0, y, 6.0);
Signal * sig=subPool->at(ig);
w->xAxis()->setTitle("Time (s)");
w->xAxis()->setNumberPrecision(2);
w->xAxis()->setTitleInversedScale("1/Time (1/s)");
w->xAxis()->setNumberPrecisionInversedScale(2);
w->yAxis()->setTitle(sig->nameComponent());
w->yAxis()->setNumberType(Number::Scientific);
w->yAxis()->setNumberPrecision(1);
w->yAxis()->setTitleInversedScale(sig->nameComponent());
w->yAxis()->setNumberPrecisionInversedScale(1);
PlotLine2D * line;
LineLayer * layer;
layer=new LineLayer(w);
line=new PlotLine2D;
line->setPen(Pen(Qt::black, 0.6, Qt::SolidLine));
layer->setReferenceLine(line);
layer->setObjectName("average");
_averageCurveLayers[ig]=layer;
layer=new LineLayer(w);
line=new PlotLine2D;
line->setPen(Pen(Qt::black, 0.6, Qt::DotLine));
layer->setReferenceLine(line);
layer->setObjectName("stddev");
_stddevCurveLayers[ig]=layer;
layer=new LineLayer(w);
line=new PlotLine2D;
line->setPen(Pen(Qt::red, 0.3, Qt::SolidLine));
layer->setReferenceLine(line);
layer->setObjectName("fit");
_fitCurveLayers[ig]=layer;
_comments[ig]=addText(16, y, 4, 5.0);
_results[ig].name=sig->nameComponent();
y+=6.0;
}
}
| double DampingResults::misfitFunk | ( | double | x[] | ) | [static, protected] |
Damping sinusoide to fit to calculated data, returns a misfit Used with call back in amoeba
References _aveCurve, _nSampFit, A, QGpCoreTools::exp(), misfit(), QGpCoreTools::sin(), TRACE, w, and Z.
Referenced by compute().
| void DampingResults::save | ( | QString | fileName = QString::null | ) | [slot] |
References _log, _results, QGpCoreTools::endl(), MSG_ID, QGpCoreTools::tr(), and TRACE.
Referenced by DampingResults().
{
TRACE;
if(fileName.isEmpty()) {
fileName=Message::getSaveFileName(tr("Save results"), tr("Damping results (*.damping)"));
}
if(fileName.isEmpty()) {
return;
}
QString message;
QFileInfo fi(fileName);
if(fi.exists()) {
message += tr( "File %1 already exists.\n" ).arg(fi.filePath());
}
QDir d(fi.absolutePath());
QFileInfo filog(d.absoluteFilePath(fi.baseName() + ".log" ));
if(filog.exists()) {
message += tr( "File %1 already exists.\n" ).arg(filog.filePath());
}
if(!message.isEmpty()) {
if(Message::warning(MSG_ID, tr( "Saving result files" ),
tr( "%1\nDo you want to replace it(them)?" ).arg(message),
Message::no(), Message::yes(), true)==Message::Answer0) return;
}
// Save a log file with parameters
QFile flog(filog.filePath());
if( !flog.open(QIODevice::WriteOnly) ) {
Message::warning(MSG_ID, tr( "Saving log file" ),
tr( "Unable to open file %1 for writing" ).arg(filog.filePath()), Message::cancel(), true);
return;
}
QTextStream slog(&flog);
slog << "### Parameters ###\n"
<< _log
<< "### End Parameters ###" << endl;
// Save a result file with frequency, damping, number of windows
QFile fres(fileName);
if( !fres.open(QIODevice::WriteOnly) ) {
Message::warning(MSG_ID, tr( "Saving result file" ),
tr( "Unable to open file %1 for writing" ).arg(fileName), Message::cancel(), true);
return;
}
QTextStream sres(&fres);
sres << "# Name | Frequency | Damping | Number of windows" << endl;
for(QVector<Result>::iterator it=_results.begin();it!=_results.end(); it++) {
sres << it->name << "\t"
<< QString::number(it->frequency) << "\t"
<< QString::number(it->damping) << "\t"
<< QString::number(it->nWindows) << endl;
}
}
| void DampingResults::setWindowTitle | ( | QString | title | ) | [inline] |
{GraphicSheetMenu::setWindowTitle(title);}
Curve< Point2D > DampingResults::_aveCurve [static, protected] |
Referenced by compute(), and misfitFunk().
LineLayer** DampingResults::_averageCurveLayers [protected] |
Referenced by compute(), createObjects(), and ~DampingResults().
TextEdit** DampingResults::_comments [protected] |
Referenced by compute(), createObjects(), and ~DampingResults().
LineLayer** DampingResults::_fitCurveLayers [protected] |
Referenced by compute(), createObjects(), and ~DampingResults().
double DampingResults::_fitW = 0 [static, protected] |
Referenced by compute().
QString DampingResults::_log [protected] |
int DampingResults::_nSampFit = 0 [static, protected] |
Referenced by compute(), and misfitFunk().
QVector<Result> DampingResults::_results [protected] |
Referenced by compute(), createObjects(), and save().
LineLayer** DampingResults::_stddevCurveLayers [protected] |
Referenced by compute(), createObjects(), and ~DampingResults().
| QMenu* DampingResults::menuTools |
Referenced by DampingResults().
