Derivative Chart Source lazarus / 55-TAChart ---------------------------------------------------------------------- wp · 2011-09-04 18:05 In one of my Lazarus projects I needed to plot the derivative of a data curve. Since such a procedure might be nice for reuse in other projects I decided to write a DerivativeChartSource component for TAChart which, along with a demo project, is enclosed in the attachment. Differentiation is done by a simple difference formula as described in http://en.wikipedia.org/wiki/Numerical_differentiation. There are better methods such as the Savitzky-Golay filtering method (http://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_smoothing_filter) , but I don't know yet whether the needed fitting routines are available in fpc. I am donating the TDerivativeChartSource to the TAChart development. I am aware that it may be too specific for line series, and it neglects the YList of the ChartSource. One particular issue exists when the derivative cannot be calculated, e.g. when the x difference is zero (see option "Data from Table" in the attached demo, at x=5). In this case, it would be desirable to skip that particular invalid data point in the chart. In Delphi's TChart, this can be done by assigning the color clNone to a data item. But this does not seem to work with TAChart. Feel free to modify the component to better stream-line it with the TAChart project. ---------------------------------------------------------------------- Ask · 2011-09-04 18:39 I have taken a quick look at your code, will look better tomorrow. Two quick points: 1) Calculating a derivative from the tabulated data seems suboptimal. Is your data really just a list of points, or is it actually a function? In the second case, maybe consider using TFuncSeries instead and expanding it to support derivatives. 2) As a small hint, "Length(Array)-1" can be spelled shorter as "High(Array)" ---------------------------------------------------------------------- wp · 2011-09-04 20:23 QuoteIs your data really just a list of points, or is it actually a function? Yes, the data are from a particular measurement of a voltage as a function of time. The data lie on a straight line which gently changes its slope towards the end of the measurement. Looking at the derivative, that change is much more easily seen. Therefore, my main intention is to calculate the derivative for tabulated data. BTW, maybe the derivative would fit into the scope of the TCalculatedChartSource. However, I got confused by the HistoryBuffer and the percentage calculations, and decided to create the DerivativeChartSource on its own. ---------------------------------------------------------------------- Ask · 2011-09-05 11:33 Quote from: wp on September 04, 2011, 10:23:38 pmmaybe the derivative would fit into the scope of the TCalculatedChartSource. That is correct, I have implemented it in r32184 (set AccumulationMethod=camDerivative and AccumulationRange=2 to 7). For an example see "Statistics" page of the line demo since r32185. Note that I have implemented backwards finite difference, while your code calculated central one. Correct calculation of the central difference needs more complex code, so maybe backwards is enough? Also note that AccumulationRange controls the order of accuracy. Quote from: wp on September 04, 2011, 10:23:38 pmHowever, I got confused by the HistoryBuffer and the percentage calculations, and decided to create the DerivativeChartSource on its own. The calculated source must access previous chart items, which may be inefficient, as described here: http://wiki.lazarus.freepascal.org/TAChart_document ation#Source_optimization_notes. So the history buffer caches last AccumulationRange items, and tries to update the cache intelligently if the calculated source itself is accessed sequentially. ---------------------------------------------------------------------- marcov · 2011-09-06 07:59 There is basic function derivation code available in unit symbolic ---------------------------------------------------------------------- Ask · 2011-09-06 08:16 I am aware of that, but I can not see how it may be useful in this case. What do you propose? ---------------------------------------------------------------------- wp · 2011-09-06 11:57 The CalculatedChartSource with the Derivative option works nicely and will replace my home-made component. Thank you. A few comments, though. I am attaching a modified version of my demo program, you can switch the series datasource between my DerivativeDataSource and the new CalculatedChartSource: - The allowed AccumulationRange 2..7 conflicts with the default value of 1. The derivative series is not displayed with this value. - Higher AccumulationRanges create more noisy data, see the option "noisy Gaussian" in the demo. This seems to be due to the calculation method. A better way would be to fit a polynomial to the data within the AccumulationRange. This leads to the Savitzky-Golay filtering. I found that the fpc numlib does contain some procedures on solving matrices, maybe I can adapt the algorithm in "Numerical recipes", or maybe I simply find the filter coefficients somwhere to provide a table like you do. - The calculation method of the derivative has issues at the beginning of the data set maybe when there are not enough data to fill the accumulation range. Select the sine or exp example, and toggle the chart source: you will see a discontinuity at the left of the series whenever CalculatedChartSource is active. - It seems that your implementation assumes equidistant x values. This may be tolerable for display purposes, but may lead to confusion otherwise. I am aware that even the algorithm in "Numerical Recipes" assumes equidistant x values. But, at least, you should add a note on that possible issue to the documentation. ---------------------------------------------------------------------- Ask · 2011-09-07 10:23 QuoteThe allowed AccumulationRange 2..7 conflicts with the default value of 1Default changed in r32201. QuoteHigher AccumulationRanges create more noisy dataThat is more or less expected, since "higher accuracy" in this case translates to "more sensitivity to the noise". QuoteA better way would be to fit a polynomial to the data within the AccumulationRange. This leads to the Savitzky-Golay filtering. I found that the fpc numlib does contain some procedures on solving matrices, maybe I can adapt the algorithm in "Numerical recipes", or maybe I simply find the filter coefficients somwhere to provide a table like you do. Yes, I think the table would be the best option from both performance and code complexity points of view. If you do calculate the coefficients, please use backward difference (aka "casual filter" in Numerical Recipes). Savitzky-Golay can also be used for smoothing without derivatives. QuoteThe calculation method of the derivative has issues at the beginning of the data set maybe when there are not enough data to fill the accumulation range. Yes, that is a problem. I imagine several solutions: 1) Calculate the single-sided derivative -- this is what you have done,   and I implemented similar change in r32204. 2) Less conveniently, but perhaps mathematically more correctly,   would be to skip first data point. As you have already noted,   TAChart does not yet provide point skipping mechanism,   so I'll have to implement it first. 3) Implement central (and forward) differences -- somewhat hard,   but will also be useful for summation and averaging.   Calculation may then gracefully downgrade accuracy near the edges. Quoteyou should add a note on that possible issue to the documentationDone, see http://wiki.lazarus.freepascal.org/TAChart_docu mentation#Calculated_source ---------------------------------------------------------------------- wp · 2011-09-07 10:46 QuoteDefault changed in r32201. Doesn't this interfere with the other usages of the TDerivedChartSource? Maybe better to leave the default at 1, but to force the value into the range allowed. ---------------------------------------------------------------------- Ask · 2011-09-07 12:15 The historical reason for default value of 1 is that the calculated source performs several related, but different functions -- reordering, accumulation and percentage calculation. Theoretically the proper OO design would be to create a separate component for each function. However, the number of components in the TAChart package is already 16, and growing, so I tried to compromise. AccumulationRange=1 makes sense for the cases where the accumulation is not needed. Fortunately, those cases may be distinguished by AccumulationMethod=camNone, so (I think) the change is safe. ---------------------------------------------------------------------- Ask · 2011-09-13 17:48 Quote from: wp on September 04, 2011, 08:05:36 pmOne particular issue exists when the derivative cannot be calculated, e.g. when the x difference is zero (see option "Data from Table" in the attached demo, at x=5). In this case, it would be desirable to skip that particular invalid data point in the chart. In Delphi's TChart, this can be done by assigning the color clNone to a data item. But this does not seem to work with TAChart. Since r32323, you can skip points by setting X or Y value to NaN (see http://wiki.lazarus.freepascal.org/TAChart_documentation#Skipping _source_items and "nan" demo). Since r32324, the leftmost point and the points where dx=0 are skipped using this mechanism. Please test. Note that currently not all series support point skipping -- I have actually tested only those in the demo. ---------------------------------------------------------------------- wp · 2011-09-13 22:38 That's very nice, now I can continue with the derivatives... http://www.holoborodko.com/pavel/numerical-methods/numerical- derivative/smooth-low-noise-differentiators/ contains a series of coefficients for noise-optimized derivatives. The one-sided differentiators, however, lead to significant distortion of my test function. I've made some progress also with the Savitzky-Golay method. Applying the "ipfpol" procedure from the fpc numlib directory it is easy to fit any polynomial to the filter-generating array. The elements of this array are all zeros, except for a one at the point to be processed. Usually, the one is at the center of the window. Using the fitted polynomial coefficients, one can calculate the function values or the 1st, 2nd etc derivatives of the fitted polynomial at the array points. These function values are used like the coefficient table in your CalculatedChartSource. Of course, this fitting needs to be done only once when characteristic parameters, such as the degree of the fitting polynomial, or the window width, are changed. Unfortunately, differentiation by Savitzky-Golay does not seem to work when the one is the first array element which would correspond to a backward filter. My question: Do you have an idea how to implement random data access in the CalculatedChartSource? Or would it better to put these filtering calculations into a separate chart source which has the data in a sorted buffer? ---------------------------------------------------------------------- Ask · 2011-09-13 23:59 QuoteDo you have an idea how to implement random data access in the CalculatedChartSource? Or would it better to put these filtering calculations into a separate chart source which has the data in a sorted buffer? I do not quite understand you. CalculatedChartSource does implement random access, it is just slower than sequential one [O(N^2) vs O(N)]. If you want to cache the result for faster access, you can use ListChartSource.CopyFrom procedure. Or do you need something else? ---------------------------------------------------------------------- wp · 2011-09-14 10:04 So this was a misunderstanding on my side of one of your previous postings. ---------------------------------------------------------------------- Ask · 2011-09-18 15:48 Implemented AccumulationDirection. Settings AccumulationDirection=cadCenter, AccumulationRange=2 are mostly equivalent to your original code. See chartsource demo since r32408. ---------------------------------------------------------------------- wp · 2011-09-18 19:38 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. Code: [Select]procedure TChartSourceSavGolFilter.CalcCoeffs(AWindowWid th,APolyOrder,ADeriv:integer); { coefficients from H.H.Madden, Anal.Chem. vol 50 (1978) p.1383   http://nsm1.nsm.iup.edu/jford/courses/CHEM421/Resources/CommentsOnTh eSavitzkyGolayConvolutionMethodForLeastSquaresFitSmoothingAndDifferent iationOfDigitalData.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; ---------------------------------------------------------------------- Ask · 2011-09-18 21:37 Thanks, I shall look at your code. Meanwhile, I implemented smoothed derivative based on Pavel Holobrodko's site you have suggested -- see r32412 ---------------------------------------------------------------------- wp · 2011-09-24 10:29 Thank you for the smoothed derivative. It works perfectly - even the backward differences.