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

for( i = 0; i < hn-nx; i++ )

betta[i] = bettacoeff*EMF[hn-2]/((HT[i+1]-HT[i])*(HT[i+1]-HT[i]));

for( i = hn-nx; i < hn; i++ )betta[i] = 0;

fp = previp = 0;

iter = 1;

while( fp < hn-nx && iter <= MaxIter )

{

CholeskySR( hn, A, betta, hn-nx, M );

if( !CholeSolveSR( hn, M, betta, F, DEMF, hn-nx, HEMF ) )return errCannotSolveSLAE;

fp = hn-nx;

for( i = 1; i < hn-nx && fp == hn-nx; i++ )

if( (1.-maxinc)*HEMF[i] > HEMF[i-1] )fp = i-1;

if( fp < previp )

for( i = fp; i < previp; i++ )betta[i] *= mulcoeff;else

for( i = fp; i < hn-nx; i++ )betta[i] *= mulcoeff;

previp = fp;

if( ProgressFunc != NULL )

if( !(*ProgressFunc)( paContinue, fp*100/(hn-nx) ) )return errUserBreak;

iter++;

}

if( iter > MaxIter )return errRegFailed;

return errOk;

}

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

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

{

int i, j, k;

double t,srct,sum;

int oldn;

double oldHTn,oldHEn;

int res;

res = DecomposeImpuls();

if( res != errOk )return res;

if( m < 2 || hn < 2 )return errTooFewPoints;

oldn = emfn = hn - 1;

t = HT[oldn];

while( emfn > 0 && HT[emfn-1]+Coeff[1]-Coeff[2*m-1] > t )emfn--;

if( emfn < 2 )return -1;

oldHTn = HT[oldn]; oldHEn = HEMF[oldn];

HT[oldn] = HT[emfn-1]+Coeff[1]-Coeff[2*m-1];

HEMF[oldn] = HEMF[oldn-1]*pow( HT[oldn-1], 2.5 )*pow( HT[oldn], -2.5 );

if( ProgressFunc != NULL )

(*ProgressFunc)( paSumHvs, 0 );

for( j = 0; j < emfn; j++ )

{

if( ProgressFunc != NULL )

if( !(*ProgressFunc)( paContinue, j*100/(emfn-1) ) )return errUserBreak;

srct = HT[j] + Coeff[1];

sum = Coeff[0]*HEMF[j];

k = j+1;

for( i = 1; i < m; i++ )

{

t = srct - Coeff[2*i+1];

while( HT[k] < t )k++;

k--;

t = (t - HT[k])/(HT[k+1] - HT[k]);

sum += Coeff[2*i]*(1-t)*HEMF[k] + Coeff[2*i]*t*HEMF[k+1];

}

aEMF[j] = sum;

}

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

aET[i] = HT[i] + Coeff[1];

HT[oldn] = oldHTn; HEMF[oldn] = oldHEn;

return errOk;

}

int CEMFTransform::Smooth( double *HT, double *HEMF, int hn, double *SmoothedEMF, double alpha )

{

double *logHT, *logHEMF;

double *SplineNodes;

int i,j,dj;

int monoton, nodenum, prev_node;

int off;

if( !IntMem )return errNoMemAllocated;

logHT = (double *)&IntMem[(off=2*m*sizeof(double))];

logHEMF = (double *)&IntMem[(off+=hn*sizeof(double))];

SplineNodes = (double *)&IntMem[(off+=hn*sizeof(double))];

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

{

if( HT[i] <= 0 || HEMF[i] <= 0 )return errNoLogScale;

logHT[i] = log10(HT[i]);

logHEMF[i] = log10(HEMF[i]);

}

SplineNodes[0] = logHT[0];

nodenum = 1; prev_node = 0;

i = 0;

while( i < hn-1 )

{

j = i + 1;

while( j < hn-1 && fabs(logHEMF[i] - logHEMF[j]) < grouplevel )j++;

while( j < hn-1 && fabs(logHEMF[j] - logHEMF[j+1]) < 0.01 )j++;

dj = j;

while( dj > i+1 && fabs(logHEMF[dj-1] - logHEMF[dj]) < 0.01 )dj--;

i = (j + dj)/2;

if( i > prev_node+2 )

{

SplineNodes[nodenum++] = logHT[i];

prev_node = i;

}

}

SplineNodes[nodenum-1] = logHT[hn-1];

off += nodenum*sizeof(double);

monoton = 0;

while( !monoton )

{

if( nodenum >= 4 )

{

if( !SmoothSpline( logHT, logHEMF, hn, SplineNodes, nodenum,

alpha, SmoothedEMF, &IntMem[off] ) )

return errCannotSolveSLAE;

i = j = 1;

monoton = 1;

while( i < hn && monoton )

{

monoton &= SmoothedEMF[i] < SmoothedEMF[i-1];

if( logHT[i] > SplineNodes[j] )j++;

i++;

}

if( !monoton )

{

nodenum--;

for( i = j; i < nodenum; i++ )SplineNodes[i] = SplineNodes[i+1];

}

}else return errTooFewSplineNodes;

}

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

SmoothedEMF[i] = pow( 10, SmoothedEMF[i] );

return errOk;

}

int CEMFTransform::LoadCoeff( istream& is, int cn )

{

if( !IntMem )return errNoMemAllocated;

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

for( int i = 0; i < cn; i++ )

is >> Coeff[2*i] >> Coeff[2*i+1];

m = cn;

return errOk;

}

int CEMFTransform::SaveCoeff( ostream& os )

{

int res;

res = DecomposeImpuls();

if( res != errOk )return res;