Thank you. The progress is very exciting.
To complete this I'd like to share some code on Savitzky-Golay smoothing and differentiation. Unfortunately, we're running out of sync with the TASources unit, and so I am posting only the essential procedure for calculation of the filter coefficients.
By chance I googled that old paper by Madden (reference in the code below) which contains pre-calculated expressions for the most important filter coefficients for any filter length.
- AWindowWidth is the length of the filtering array which is moving over the data. The central point of that array corresponds to the point to be filtered. AWindowWidth must be an odd number. It is related to your AccumulationRange, but I did not yet test the new CalculatedChartSource.
- ADeriv is the order of the derivative (0 = smoothing). The paper goes up to ADeriv=5, but I stopped typing at ADeriv=2 for the second derivative -- the others are less important. I extended the TAccumulationMethod by camDerivative2 and camSmooth to tell which one is needed.
- APolyOrder is the order of the fitting polynomial, must be less than AWindowWidth - 1. TCalculatedChartSource should get a corresponding property for that.
- The coefficients go into a double array FCoeffs which you can plug directly into the CalcDerivative method. As already noted, the point to be fitted corresponds to the center of the FCoeffs array, so the WeightedSum must run upwards and downwards. The AWindowWidth div 2 points at the edges of the data range are dropped for simplicity (set to NaN). Of course, it should be possible to use a window that gets narrower the closer the filtered index is at the edge.
BTW: The code is not quite efficient since I am filling the entire coefficient array. Due to symmetry, it would be better to use only the (AWindowSize div 2 + 1) elements, like you do with the finite difference arrays.
The CalcCoeffs procedure below is called only once after changing of AccumulationRange, PolyOrder, or AccumulationMethod. I am keeping the filter coefficients in a separate class which I FreeAndNil at these circumstances to signal an invalid filter. But a simple boolean flag will do as well. In the CalcDerivative method I check if the filter coefficients are valid, if not, CalcCoeffs gets called.
See the enclosed screenshots as an example on the impressive noise behavior of the Savitzky-Golay differentiation.
procedure TChartSourceSavGolFilter.CalcCoeffs(AWindowWidth,APolyOrder,ADeriv:integer);
{ coefficients from H.H.Madden, Anal.Chem. vol 50 (1978) p.1383
http://nsm1.nsm.iup.edu/jford/courses/CHEM421/Resources/CommentsOnTheSavitzkyGolayConvolutionMethodForLeastSquaresFitSmoothingAndDifferentiationOfDigitalData.pdf
}
function denom1(n,m1,m2:integer) : double;
var
m : integer;
begin
result := 1.0;
m := m1;
while m >= m2 do begin
result := result * (n + m);
dec(m);
end;
end;
function denom2(n,m1,m2:integer) : double;
var
m : integer;
begin
result := 1.0;
m := m1;
while (m >= m2) do begin
result := result * (2*n + m);
dec(m, 2);
end;
end;
var
j : integer;
nf:integer;
nf2,nf3,nf4,nf5,nf6,nf7,nf8 : double;
begin
nf := AWindowWidth div 2;
nf2 := nf*nf;
SetLength(FCoeffs, AWindowWidth);
for j:=0 to AWindowWidth-1 do
FCoeffs[j] := 0;
case ADeriv of
0 : case APolyOrder of
2, 3 : // ok
for j:=0 to nf do begin
FCoeffs[nf-j] := 3*(3*nf2+3*nf-1-5*j*j)/denom2(nf, 3,-1);
FCoeffs[nf+j] := FCoeffs[nf-j];
end;
4, 5 : // ok
begin
nf3 := nf2*nf;
nf4 := nf3*nf;
for j:= 0 to nf do begin
FCoeffs[nf-j] := 15/4 *
( (15*nf4+30*nf3-35*nf2-50*nf+12) - 35*(2*nf2+2*nf-3)*sqr(j) + 63*power(j,4) ) /
( denom2(nf, 5, -3) );
FCoeffs[nf+j] := FCoeffs[nf-j];
end;
end else
raise Exception.Create('Smoothing by Savitzky-Golay only implemented for polyonomial order < 6');
end;
1 : case APolyOrder of
2 : // ok
for j:=1 to nf do begin
FCoeffs[nf-j] := 3*j/((2*nf+1)*(nf+1)*nf);
FCoeffs[nf+j] := -FCoeffs[nf-j];
end;
3,4 : //ok
begin
nf3 := nf2*nf;
nf4 := nf3*nf;
for j:=1 to nf do begin
FCoeffs[nf-j] := 5*
( 5*(3*nf4+6*nf3-3*nf+1)*j - 7*(3*nf2+3*nf-1)*power(j,3) ) /
( denom2(nf, 3,-1) * denom1(nf, 2,-1) );
FCoeffs[nf+j] := -FCoeffs[nf-j];
end;
end;
5, 6 : // ok
begin
nf3 := nf2*nf;
nf4 := nf3*nf;
nf5 := nf4*nf;
nf6 := nf5*nf;
nf7 := nf6*nf;
nf8 := nf7*nf;
for j:=1 to nf do begin
FCoeffs[nf-j] := 21/4 *
( 7*(25*nf8+100*nf7-50*nf6-500*nf5-95*nf4+760*nf3+180*nf2-300*nf+72)*j
-105*(6*nf6+18*nf5-15*nf4-60*nf3+17*nf2+50*nf-12)*power(j,3)
+33*(15*nf4+30*nf3-35*nf2-50*nf+12)*power(j,5) ) /
( denom2(nf, 5, -3) * denom1(nf, 3, -2) );
FCoeffs[nf+j] := -FCoeffs[nf-j];
end;
end;
else
raise Exception.Create('Derivative Savitzky-Golay only implemented for polynomial orders < 7');
end;
2 : case APolyOrder of
2, 3 :
begin
for j := 0 to nf do begin
FCoeffs[nf-j] := 30*(3*j*j - nf*(nf+1)) / (denom2(nf, 3,-1) * denom1(nf, 1,0));
FCoeffs[nf+j] := FCoeffs[nf-j];
end;
end;
4, 5 :
begin
nf3 := nf2*nf;
nf4 := nf3*nf;
nf5 := nf4*nf;
for j:=0 to nf do begin
FCoeffs[nf-j] := -105/2 *
( 15*(6*nf2+6*nf-5)*power(j,4) - 21*(4*nf4+8*nf3-4*nf2-8*nf+5)*j*j
+5*nf*(2*nf5+6*nf4-nf3-12*nf2-nf+6) ) /
( denom2(nf, 5,-3)*denom1(nf, 2,-1) );
FCoeffs[nf+j] := FCoeffs[nf-j];
end;
end;
else
raise Exception.Create('2nd derivative only implemented for polynomial orders < 6');
end;
else
raise Exception.Create('2nd order derivatives or higher not implemented.');
end;
end;
Recent
