{Data and signal processing utilities
author: Nikolai V. Shokhirev
created: 2002.02.02
last modified: 2004.12.26
ŠNikolai V. Shokhirev, 2002-2004 }
unit uDynDataUtils;
interface
uses
uMatTypes, uDynArrays, uDynData;
type
TPoint2 = 1..2;
TPoint3 = 1..3;
TPoint4 = 1..4;
TPoint5 = 1..5;
{ (dn/d Base)**n (Base**Exponent) }
function dnPower(const Base: TFloat; const Exponent, n: Integer): TFloat;
{ df/dx at xn ; *** does not depend on n *** }
function FD2d1(x1,x2,y1,y2: TFloat; n: TPoint2): TFloat;
{ df/dx at xn }
function FD3d1(x1,x2,x3,y1,y2,y3: TFloat; n: TPoint3): TFloat;
{ d2f/dx2 at xn }
function FD3d2(x1,x2,x3,y1,y2,y3: TFloat; n: TPoint3): TFloat;
{ df/dx at xn ; *** does not depend on n *** }
function FD4d1(x1,x2,x3,x4,y1,y2,y3,y4: TFloat; n: TPoint4): TFloat;
{ d2f/dx2 at xn }
function FD4d2(x1,x2,x3,x4,y1,y2,y3,y4: TFloat; n: TPoint4): TFloat;
{ d3f/dx3 at xn ; *** does not depend on n *** }
function FD4d3(x1,x2,x3,x4,y1,y2,y3,y4: TFloat; n: TPoint4): TFloat;
{ 5-point 2nd degree average.
INPUT:
y1, y2, y3, y4, y5 - function values at equi-distant points
OUTPUT:
averaged (smoothed) value at the point n
c1 c2 c3 c4 c5 n
( 31, 9, -3, -5, 3)/35 1
^
( 9, 13, 12, 6, -5)/35 2
^
( -3, 12, 17, 12, -3)/35 3
^
( -5, 6, 12, 13, 9)/35 4
^
( 3, -5, -3, 9, 31)/35 5
^ }
function Avg5(const y1, y2, y3, y4, y5: TFloat; n: TPoint5): TFloat;
{ 5-point 2nd degree smoothed derivative.
INPUT:
y1, y2, y3, y4, y5 - function values at equi-distant points
h - step
OUTPUT:
1st derivative value at the point n
c1 c2 c3 c4 c5 n
(-54, 13, 40, 27,-26)/70 1
^
(-34, 3, 20, 17, -6)/70 2
^
( -2, -1, 0, 1, 2)/10 3
^
( 6,-17,-20, -3, 34)/70 4
^
( 26,-27,-40,-13, 54)/70 5
^ }
function Avg5d1(const y1, y2, y3, y4, y5, h: TFloat; n: TPoint5): TFloat;
{ 5-point 2nd degree smoothed 2nd derivative.
INPUT:
y1, y2, y3, y4, y5 - function values at equi-distant points
h - step
OUTPUT:
2nd derivative value (does not depend on n)
c1 c2 c3 c4 c5
( 2, -1, -2, -1, 2)/7 }
function Avg5d2(const y1, y2, y3, y4, y5, h: TFloat; n: TPoint5): TFloat;
{ 2-point forward 1st derivative: (Ain[i+1]-Ain[i])/h. Ain can be equal to D1Aout. }
procedure deriv1f(const Ain: IFSignal1D; const Aout: IFArr1D);
{ 2-point backward 1st derivative: (Ain[i]-Ain[i-1])/h. Ain can be equal to D1Aout. }
procedure deriv1b(const Ain: IFSignal1D; const Aout: IFArr1D);
{ 3-point centered 1st derivative: (Ain[i+1]-Ain[i-1])/(2*h). Ain can be equal to D1Aout. }
procedure deriv1c(const Ain: IFSignal1D; const Aout: IFArr1D);
{ 3-point 2nd derivative: {Ain[i+1]-2*Ain[i]+Ain[i-1])/(h*h). Ain can be equal to D1Aout. }
procedure deriv2c(const Ain: IFSignal1D; const Aout: IFArr1D);
{ 5-point average. Ain can be equal to Aout
( 31, 9, -3, -5, 3)/35 1
^
( 9, 13, 12, 6, -5)/35 2
^
( -3, 12, 17, 12, -3)/35 3...N-2
^
( -5, 6, 12, 13, 9)/35 N-1
^
( 3, -5, -3, 9, 31)/35 N
^ }
procedure Avg5Arr(const Ain, Aout: IFArr1D);
{ 5-point average derivative. Ain can be equal to Aout
(-54, 13, 40, 27,-26)/70 1
^
(-34, 3, 20, 17, -6)/70 2
^
( -2, -1, 0, 1, 2)/10 3...N-2
^
( 6,-17,-20, -3, 34)/70 N-1
^
( 26,-27,-40,-13, 54)/70 N
^ }
procedure Avg5Arrd1(const Ain: IFSignal1D; const Aout: IFArr1D);
{ 5-point average the second derivative. Ain can be equal to Aout
( 2, -1, -2, -1, 2)/7 }
procedure Avg5Arrd2(const Ain: IFSignal1D; const Aout: IFArr1D);
{ replaces a smooth function with a step-wise one }
procedure Sharpen(const Ain, Aout: IFArr1D);
{3-point derivative }
procedure DerivData1D(const Ain, Aout: IFData1D);
{3-point 2-nd derivative }
procedure Deriv2Data1D(const Ain, Aout: IFData1D);
{
PURPOSE:
Generation of the coefficients for smoothing over M = 2*n+1 points using
the least squares cubic (degree = 3) or quadratic (degreet = 2) fit
INPUT:
n - smoothing half-interval (M = 2*n+1)
order = 0 - function (Savitzky-Golay smoothing filter)
= 1 - 1st derivative, Avg should be divided by h, space step
= 2 - 2nd derivative, Avg should be divided by sqr(h)
OUTPUT:
c - the matrix of smoothing coefficients
c[0,m] - symmetric formula, used for the center of the smoothing interval
-m, ... ,-2,-1, 0, 1, 2, ... ,m
^
c[-k,m] - offset k to the left, used for the first n points
-m, ... ,-2,-1, 0, 1, 2, ... ,m
^
c[+k,m] - offset k to the right, used for the last n points
-m, ... ,-2,-1, 0, 1, 2, ... ,m
^ }
function GenSmoothCoeff(const n, order: TInt; degree: TInt = 3):IFArr2D;
{
PURPOSE:
Calculation of smoothed values
NOTE:
This function calculates the smoothed values of a function (Savitzky-Golay
smoothing) or its derivatives depending on the generated coefficients
INPUT:
idx - the center of the smoothing interval
offset < 0 - offset to the left
> 0 - offset to the right
c - coefficients generated in GenSmoothCoeff
OUTPUT:
<f[(idx-offset)*h]> if order = 0 (Savitzky-Golay smoothing filter)
<h*f'[(idx-offset)*h]> if order = 1
<h*h*f"[(idx-offset)*h]> if order = 2 }
function Avg(const y: IFArr1D; const idx, offset: TInt; const c:IFArr2D): TFloat;
implementation
uses
math;
end.
Back
Generated by Lore's Source to HTML Converter(http://www.newty.de/lsc/index.html)