Разложение импульса.Разложение ЭДС.Использование регуляризации, страница 13

// Данная процедура допускает совпадение адресов ET/HT и Eds/HEds,

// но не частичное перекрывание этих областей.

//

//  Вход:

//   Coeff - коэффициенты разложения (см. DecomposeImpuls)

//       m - кол-во коэффициентов

//      HT - времена кривой ЭДС хевисайда

//    HEMF - значения ЭДС хевисайда, в соответсвующих временах

//      hn - количество точек кривой ЭДС хевисайда

//

//  Выходные:

//     aET - времена для кривой ЭДС, соответсвующей импульсу

//    aEMF - значения ЭДС

//    emfn - количество точек новой кривой ЭДС,

//           обычно emfn = hn - 1

//

int MakeImpulsEMF( double *HT, double *HEMF, int hn,

double *aET, double *aEMF, int &emfn );

int Smooth( double *HT, double *HEMF, int hn, double *SmoothedEMF, double alpha = 0 );

double GetEpsI() { return epsi; }

double GetCThick() { return cthick; }

void SetSelfAlloc( int value )

{

if( SelfAlloc && IntMem )delete[] IntMem;

SelfAlloc = value;

SetIntMem( NULL );

}

void SetDecomposeParams( double aepsi, double acthick )

{

epsi = aepsi;

cthick = acthick;

Coeff = NULL; m = 0;

SetIntMem( NULL );

}

void SetImpuls( double *aIT, double *aII, int ain )

{

IT = aIT;

II = aII;

in = ain;

Coeff = NULL; m = 0;

SetIntMem( NULL );

}

void SetEMFCurve( double *aET, double *aEMF, int aen )

{

ET = aET;

EMF = aEMF;

n = aen;

Coeff = NULL; m = 0;

SetIntMem( NULL );

}

void SetIntMem( char *mem );

int LoadCoeff( istream& is, int cn );

int SaveCoeff( ostream& os );

};

char* ErrorMessage( int err );

#endif

Файл hvs.cpp

#include <mem.h>

#include <math.h>

#include <iostream.h>

#include <fstream.h>

#include "hvs.h"

#include "misc.h"

// !!! Описание всех процедур находятся в файле hvs.h !!!

#define min(a,b) (a < b ? a : b)

#define max(a,b) (a > b ? a : b)

#define sign(x)  (x > 0 ? 1. : (x < 0 ? -1. : 0.))

struct TDPoint{

double x;

double y;

};

// Class CEMFTransform

void CEMFTransform::SetIntMem( char *mem )

{

if( SelfAlloc && IntMem != NULL )delete[] ( IntMem );

if( SelfAlloc && mem == NULL )

IntMem = new char[CountInternalMem()*sizeof(char)];else

{

SelfAlloc = 0;

IntMem = mem;

}

}

int CEMFTransform::CountCoefficients()

{

double eps, max = 0, cn = 0;

int i;

if( !II || !IT )return 0;

for( i = 0; i < in; i++ )

if( fabs(II[i]) > max )max = fabs(II[i]);

eps = max*epsi;

if( cthick == 1 )

{

for( i = 0; i < in-1; i++ )

cn += fabs( II[i] - II[i+1] )/eps;

}else

if( cthick > 1 )

{

for( i = 0; i < in-1; i++ )

{

cn += log(fabs(II[i]-II[i+1])*(cthick-1)/eps+1)/log(cthick);

eps = max*epsi*pow( cthick, floor(cn) );

}

}else

{

for( i = 0; i < in-1; i++ )

{

cn += log(fabs(II[i]-II[i+1])*(cthick-1)/eps+1)/log(cthick);

eps = max*epsi*pow( cthick, ceil(cn) );

}

}

return (int)(cn+4);

}

int CEMFTransform::CountInternalMem()

{

// коэффициенты + 2 треугольника (исходная матрица и разложение)

// + 2 вектора ( F & betta ) + профиль

int cn = CountCoefficients();

return ( 2*cn*sizeof(double)+

2*(n+1)*(n+2)/2*sizeof(double)+

2*(n+1)*sizeof(double)+

(n+2)*sizeof(int) );

}

int CEMFTransform::CountInternalMem( int cn )

{

return ( 2*cn*sizeof(double)+

2*(n+1)*(n+2)/2*sizeof(double)+

2*(n+1)*sizeof(double)+

(n+2)*sizeof(int) );

}

int CEMFTransform::DecomposeImpuls()

{

if( !IntMem )return errNoMemAllocated;

if( Coeff != NULL )return errOk;

if( !IT || !II )return errNoImpuls;

if( in < 2 )return errTooFewPoints;

TDPoint t1, t2, t3, cur;

double max = 0, s = 0, sum = 0;

double x, epst;

int i;

Coeff = (double *)&IntMem[0];

for( i = 0; i < in; i++ )

if( fabs(II[i]) > max )max = fabs(II[i]);

epsi = max*epsi;

epst = IT[in-1]*1e-15;

t1.x = IT[in-1]; t1.y = II[in-1];

t2.x = IT[in-2]; t2.y = II[in-2];

cur = t1;

i = in - 3;

m = 0;

if( ProgressFunc != NULL )

(*ProgressFunc)( paDecompose, 0 );

while(1)

{

if( fabs( cur.y - t2.y ) >= epsi )

{

t3.y = cur.y + sign(t2.y - cur.y)*epsi;

t3.x = (t2.x - t1.x)/(t2.y - t1.y)*(t3.y - t1.y) + t1.x;

if( fabs(cur.x-t3.x) > epst )