6
$\begingroup$

If I have the following data (which is temperature in the x-axis and heat flow in the y-axis):

Import["https://pastebin.com/raw/SMKZUtbQ", "Package"]

which plotted using ListLinePlot[data, PlotRange -> {{50, 100}, {-0.1, 1}}] gives:

Image

Questions:

1) How can I find the onset value of each peak?. The onset value is defined as the intersection of the tangents of the peak with the extrapolated baseline (which can be taken as 0 in this case). An example of the onset value for two peaks are given in the figure below (done in excel), where the onset would be where both red lines intersect:

onset

PS: Here's a very brief description of onset temperature if you have any doubt: http://www.hydrateweb.org/dsc#:~:text=The%20onset%20and%20the%20offset%20temperatures%20are%20defined%20as%20the,heating%20rate%20or%20sample%20preparation)

EDIT: Here's one clarification about the tangent line to use. The tangent should be the line along the most reasonable length of the peak. The criteria for “reasonable” in this scenario would be the longest length of approximately consistent slope. Here's a picture that shows another example using two tangents lines to the "left" and "right" (I only need the one on the "left" to find the "onset"):

Image

For the baseline, please use 0 as the baseline.

$\endgroup$
23
  • 2
    $\begingroup$ Have you seen FindPeaks? (reference.wolfram.com/language/ref/FindPeaks.html) $\endgroup$
    – C. E.
    Jun 5, 2020 at 21:12
  • $\begingroup$ @ C.E thank you! I think FindPeaks will definitely help me finding question 1. Thank you for the suggestion. Is there any similar function to find the onset value (temperature in this case)?. Also, do you know if I can use FindPeaks only on a given region (let's say from 65 to 94)? $\endgroup$
    – John
    Jun 5, 2020 at 21:14
  • 1
    $\begingroup$ Just take the slice of data that you’re interested in and give that to FindPeaks. As for the onset, there is no built-in function for that. $\endgroup$
    – C. E.
    Jun 5, 2020 at 21:23
  • 2
    $\begingroup$ @John Yes, I get that, but through which point is the tangent drawn? For instance, a tangent line that goes through the tip of the peak will be almost horizontal. THe web site does not define that either. I think they may have implied the tangent through the inflection point in the peak, which is pretty much the only one you can define precisely on the "side" of the peak's shape. $\endgroup$
    – MarcoB
    Jun 5, 2020 at 21:47
  • 3
    $\begingroup$ Multiple subtle choices need to be made in order to estimate 1) the curve's inflection points along with its tangent lines in the presence of noise — its derivatives turned out to be pretty noisy, and 2) the baseline. Too much freedom in the choices, no criteria of correctness but the vague common sense, no Mathematica functions to make all of the choices automagically, and next to no Mathematica users to know the specifics of this particular area. You need to find out — or invent yourself — the precise methodology of how to do it in principle before asking how to do it in Mathematica. $\endgroup$
    – aooiiii
    Jun 5, 2020 at 22:54

1 Answer 1

9
+50
$\begingroup$

John, I hope you will believe me when I tell you that you are severely underestimating the complexity of your problem. I think it laudable that you want to learn the tools of your trade better, and I definitely encourage it!

The process you describe, however, is complex and tedious to program up by hand; it also requires a great deal of technique-specific knowledge to make the "right" decisions. For instance, a piece of software to perform peak detection and integration for Nuclear Magnetic Resonance data will make different decisions than one designed for chromatography, and yet others are used in calorimetry. To try an re-create these nuances in a few lines of code is a bit naïve. Done right, your question is far more complex than you give it credit for.

I also wanted to address your point about "getting better by watching experts do something". Although that is certainly true, programming relies on a lot of trial and error. You try something; it does not work; you painfully, slowly fix your mistakes by trawling through this site and others and by reading the docs; and then whatever you learned will be seared into your brain :-)

But enough chit chat. Here is some code to illustrate a few of my points.

First off, you know that this is a question of zeros in derivatives, so first we need to calculate some derivatives. That might steer you towards interpolation:

int = Interpolation[data];

MapThread[
 Plot[
   D[int[x], {x, #1}] /. x -> t, {t, 45, 110},
   PlotRange -> All, Axes -> {True, False}, Frame -> True,
   ImageSize -> Medium, PlotLabel -> #3,
   PlotStyle -> #2] &,
 {{0, 1, 2}, 
  {Black, Red, Blue}, 
  {"interpolated data", "first derivative", "second derivative"} }
]

interpolation results

The data has a wandering baseline, but the first derivative still looks pretty good, although a bit noisy. In fact, it is a good example of why data is often presented in "derivative form" when the baseline matters little and the position of the peaks is more important (of course, the peak positions correspond to the zero-crossings in the first derivative).

The second derivative looks very noisy though. We need to find the zeroes of $f''$ because those are the positions of the inflection points of those peaks. This is too noisy though; it would be challenging to work with this. You would want to smooth it.

In fact, Savitzky - Golay smoothing would be a common choice in this case; convolution with an appropriate Savitzky - Golay kernel can give you smoother data, but also first and second derivatives directly (see Savitzky - Golay filter on Wikipedia, (124928), (37380), (190857), and SavitzkyGolayMatrix.

The thing is, your data is time-stamped and of course you would want to apply smoothing only to the ordinate, not to the times. You would also have to keep track of the points you "lose" by convolution with the filter kernel, etc. etc. Best not done by hand; fortunately, the TimeSeries machinery in Mathematica is perfect for this kind of stuff. All operations will be carried out on the intensities, and the time stamps will be correctly and automatically carried along. Creating a TimeSeries object from your data is simple: TimeSeries[data].

With that in hand, we can apply appropriate Savitzky - Golay filters to smooth the data, and to obtain its smoothed first and second derivatives:

{smoothed, firstderivative, secondderivative} = 
 ListConvolve[SavitzkyGolayMatrix[{10}, 3, #], TimeSeries[data]] & /@ 
  Range[0, 2] 

This applies a smoothing kernel of radius 10 (hand-wavingly, considering runs of ten points in your data), performs polynomial regression of degree 3 (pretty standard choice), and produces the $n^{th}$ derivative. With $n=0$ you get smoothed data, with $n=(1,2)$ you get the smoothed first and second derivatives, respectively:

We can then use DateListPlot to show the results. We can select a particular time window to plot using TimeSeriesWindow, to focus on the region at 50 to 100 seconds (or minutes, or whatever your time unit is, which you did not specify): that is where your peaks are

Here are smoothed data and first derivative:

DateListPlot[
  TimeSeriesWindow[#, {52, 105}] & /@ {smoothed, 5 firstderivative},
  PlotStyle -> {Black, Red}, PlotRange -> All,
  GridLines -> {None, {0}}, GridLinesStyle -> Darker@Gray,
  DateTicksFormat -> {"Minute", ":", "Second"},
  PlotLegends -> {"smoothed data", "first derivative"}
]

data and first derivative

... and here are smoothed data and second derivative:

DateListPlot[
  TimeSeriesWindow[#, {52, 105}] & /@ {smoothed, 30 secondderivative},
  PlotStyle -> {Black, Blue}, PlotRange -> All,
  GridLines -> {None, {0}}, GridLinesStyle -> Darker@Gray,
  DateTicksFormat -> {"Minute", ":", "Second"},
  PlotLegends -> {"smoothed data", "second derivative"}
]

data and second derivative

Much better, no?

Alright, that is a good starting point. We can work with this. So now we "only" have to:

  • find the zeroes of the first derivative (the positions of the peaks, for reference);
  • find the relevant zeroes of the second derivative (the positions of the inflection points), two for each peak ("left" and "right");
  • calculate the values of the first derivative at the inflection point, derive the equation of the tangent line through that point with that slope.
  • estimate a LOCAL baseline for each peak (there obviously is no global baseline here, since the drift is significant); perhaps subtract it from the peak?
  • calculate the intersection between baseline and that tangent.
  • repeat for the other side;
  • repeat for all peaks.

I hope to convey that this is a very complicated task. I am not going to attempt the rest, as it is laborious and time-consuming. But I do strongly encourage you to do so, if you still want to! You will learn A LOT if you do.

$\endgroup$
10
  • $\begingroup$ Very impressive! May I please ask you to comment on the problem of very sharp peaks, resembling Exp[-Abs[x]], with their two inflection points almost coinciding. I thought that symmetric smoothing filters would blur the peak, resulting in either underestimation of the derivative's absolute values, or too much noise, depending on the filter size. My idea was to convolve the signal with a causal filter strongly skewed to zero, to estimate one side of the peak, and its time reversed version for another, so that the slopes wouldn't cancel each other out. What do you think? $\endgroup$
    – aooiiii
    Jun 6, 2020 at 3:29
  • $\begingroup$ @aooiiii Smoothing will certainly blur the peak, but should not typically modify the position of the maximum; the magnitude of the derivative is rarely important in peak detection, its zeroes more so. I'd be interested in seeing your approach though, but I was wondering how often one would encounter such sharp peaks as Exp[Abs[x]]in physical processes. $\endgroup$
    – MarcoB
    Jun 6, 2020 at 4:16
  • $\begingroup$ Could it be that I misunderstood something? According to my approach, the slopes of the tangent lines are defined to be equal to the maximum and minimum of the two derivative estimates. Thus underestimating the derivative's magnitude would have resulted in flatter tangent lines. Regarding Exp[-Abs[x]], I meant that the curve, although smooth, has peaks with tiny rounded corners and inflection points a few sampling periods apart. Precisely this region, smaller than the filter size, defines the slopes of the tangent lines, as the derivative there is the highest in magnitude. I thought. $\endgroup$
    – aooiiii
    Jun 6, 2020 at 5:54
  • $\begingroup$ MarcoB this is a very impressive code. Yes, you are right. Perhaps I am not aware of how difficult (and interesting) this problem is. I will take your code as basis and if I solve it I will post it here so that other people in the community can benefit as well. Thanks $\endgroup$
    – John
    Jun 6, 2020 at 15:57
  • $\begingroup$ MarcoB do you think it is better and easier to create a code that allows the user to select the tangent (out of the many possibles) that the user thinks is correct to get the onset base on the baseline of 0 and once this is done the code will give the onset value. This would be similar to what some old thermogram softwares do as well?, Would this be more approachable to do than the complications that arise from doing it similar to how you are doing it? $\endgroup$
    – John
    Jun 6, 2020 at 17:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.