FWHM, InterpolationFunction & Solve

Discussion:

(too old to reply)

Mathieu G

2007-08-22 08:42:38 UTC

Hello,
I have a set of discrete data, representing a peak.
I would like to compute the Full Width at Half Maximum (FWHM) of this peak.
For that I would like to know which points correspond to half the peak value.
SO far I was considering using an interpolating function, but this does
not seem to work:

DataFile = Import["FFT.dat"];
CleanDataFile = Part[DataFile, 2 ;; Length[DataFile]];
WorkingData = CleanDataFile[[All, {1, 3}]];

ListLinePlot[WorkingData, PlotRange -> All]

MaximumFFTY = Max[WorkingData[[All, 2]]];
MaximumFFTPosition = Position[WorkingData, MaximumFFTY][[1, 1]];
MaximumFFTX = WorkingData[[MaximumFFTPosition, 1]];

DataInterpolation = Interpolation[WorkingData];

Solve[DataInterpolation[x] == MaximumFFTValue/2, x]

Can you help me please? How would you do that?
I then want to compute the area under the peak:

Integrate[DataInterpolation[x], {x, BegFWHM, EndFWHM}]

Which works fine with the interpolating function.

Thank you for your help!
Mathieu

David Annetts

2007-08-23 05:00:17 UTC

Permalink

Hi Mathieu

< snippage>

Post by Mathieu G
Can you help me please? How would you do that?

When I last needed to do this, Ted Ersek's RootSearch package worked very
well. The package is still available on MathSource & seems to work well
enough

Post by Mathieu G
Integrate[DataInterpolation[x], {x, BegFWHM, EndFWHM}]
Which works fine with the interpolating function.

I'll use (very noisy, just to make it interesting) synthetic data after
loading the package ....

Needs["Enhancements`RootSearch`"]
samp = {#, Sin[#] + RandomReal[{-.1, .1}]} & /@ Range[0, Pi,
Pi/128];
smax = Max[samp[[All, 2]]]
ListPlot[samp, Joined -> True, GridLines -> {None, {smax/2}}, Mesh
-> All]
sint = Interpolation[samp];

We find the boundaries of the peak at half-height using

hmax = t /. RootSearch[sint[t] == smax/2., {t,
***@First@***@sint, ***@First@***@sint}]

And the area using

slimit = Flatten[#] & /@ (Join[{t, #}] & /@ Partition[hmax, 2, 1])
areas = Plus @@ (Integrate[sint[t], #] & /@ slimit)

Which compares favourably to the correct value

areaa = Integrate[Sin[t], Join[{t}, ***@hmax]]

Limits above appear complicated, but really only account for noisy or
pathological data.

May I ask about your application? For some, anomaly width at 2/3 max might
be a more robust quantity to use than anomaly width at 1/2 height.

Regards,

Dave.

Sseziwa Mukasa

2007-08-23 05:07:25 UTC

Permalink

Post by Mathieu G
Hello,
I have a set of discrete data, representing a peak.
I would like to compute the Full Width at Half Maximum (FWHM) of
this peak.
For that I would like to know which points correspond to half the
peak value.
SO far I was considering using an interpolating function, but this
does
DataFile = Import["FFT.dat"];
CleanDataFile = Part[DataFile, 2 ;; Length[DataFile]];
WorkingData = CleanDataFile[[All, {1, 3}]];
ListLinePlot[WorkingData, PlotRange -> All]
MaximumFFTY = Max[WorkingData[[All, 2]]];
MaximumFFTPosition = Position[WorkingData, MaximumFFTY][[1, 1]];
MaximumFFTX = WorkingData[[MaximumFFTPosition, 1]];
DataInterpolation = Interpolation[WorkingData];
Solve[DataInterpolation[x] == MaximumFFTValue/2, x]
Can you help me please? How would you do that?

Use FindRoot not Solve:

FindRoot[DataInterpolation[x]==MaximumFFTValue/2,{x,MaximumFFTX}]

But if you have a model for your signal you should fit that with
FindFit rather than doing interpolation. At any rate this will only
work if there is one extremum in your data.

Regards,

Ssezi

Mathieu G

2007-08-28 06:18:23 UTC

Permalink

Post by Sseziwa Mukasa

FindRoot[DataInterpolation[x]==MaximumFFTValue/2,{x,MaximumFFTX}]
But if you have a model for your signal you should fit that with
FindFit rather than doing interpolation. At any rate this will only
work if there is one extremum in your data.
Regards,
Ssezi

Hi!
Even though I read the documentation I still do not understand the
difference between FindRoot and Solve, can you help?
Regards,
Mathieu

Sseziwa Mukasa

2007-08-29 08:20:24 UTC

Permalink

Post by Mathieu G

Post by Sseziwa Mukasa

Hi!
Even though I read the documentation I still do not understand the
difference between FindRoot and Solve, can you help?

Solve is essentially for solving linear systems, it also attempts to
solve systems involving trigonometric, exponential and expressions.
It tries to find all possible solutions to an equation which in the
case of your interpolation function may be quite a large set.
FindRoot just searches for the first solution found, which is not
necessarily the closest, given the starting point, if you have only
on peak in your data FindRoot should be able to bracket and locate
the FWHM pretty quickly.

Regards,

Ssezi

Andrzej Kozlowski

2007-08-31 03:53:36 UTC

Permalink

Post by Sseziwa Mukasa
Solve is essentially for solving linear systems, it also attempts to
solve systems involving trigonometric, exponential and expressions.

Linear??? So Groebner basis is just a fancy name for Gassian
elimination?

Andrzej Kozlowski

2007-08-23 05:15:36 UTC

Permalink

Hi Mathieu,

if your data has noise you may have some troubles with Max and Solve. In

this case a fit will be more adequate. If the peak has gaussian form I

would fit a Gaussian (take log of data to make the fit linear!). If you

have the Gaussian, you can get FWHM analytically.

hope this helps, Daniel

Post by Mathieu G
Hello,
I have a set of discrete data, representing a peak.
I would like to compute the Full Width at Half Maximum (FWHM) of this peak.
For that I would like to know which points correspond to half the peak value.
SO far I was considering using an interpolating function, but this does
DataFile = Import["FFT.dat"];
CleanDataFile = Part[DataFile, 2 ;; Length[DataFile]];
WorkingData = CleanDataFile[[All, {1, 3}]];
ListLinePlot[WorkingData, PlotRange -> All]
MaximumFFTY = Max[WorkingData[[All, 2]]];
MaximumFFTPosition = Position[WorkingData, MaximumFFTY][[1, 1]];
MaximumFFTX = WorkingData[[MaximumFFTPosition, 1]];
DataInterpolation = Interpolation[WorkingData];
Solve[DataInterpolation[x] == MaximumFFTValue/2, x]
Can you help me please? How would you do that?
Integrate[DataInterpolation[x], {x, BegFWHM, EndFWHM}]
Which works fine with the interpolating function.
Thank you for your help!
Mathieu

Kevin J. McCann

2007-08-24 06:14:09 UTC

Permalink

If your data are noisy, I would recommend the use of a Savitzky-Golay
filter. The gist of this is that a polynomial of some degree d is fit
(LSQ) to 2N+1 data points, both d and N are user determined. Then the
middle data point of the 2N+1 is replaced with the curve fit value. The
process proceeds by sliding one point to the right. All of these
operations happen on the original data.

This is a smoothing filter that removes some of the "higher frequency"
stuff on the assumption that it is due to "noise".

I have a NB with some details. If you would like, I could send it to you.

Kevin

Post by Mathieu G
Hello,
I have a set of discrete data, representing a peak.
I would like to compute the Full Width at Half Maximum (FWHM) of this peak.
For that I would like to know which points correspond to half the peak value.
SO far I was considering using an interpolating function, but this does
DataFile = Import["FFT.dat"];
CleanDataFile = Part[DataFile, 2 ;; Length[DataFile]];
WorkingData = CleanDataFile[[All, {1, 3}]];
ListLinePlot[WorkingData, PlotRange -> All]
MaximumFFTY = Max[WorkingData[[All, 2]]];
MaximumFFTPosition = Position[WorkingData, MaximumFFTY][[1, 1]];
MaximumFFTX = WorkingData[[MaximumFFTPosition, 1]];
DataInterpolation = Interpolation[WorkingData];
Solve[DataInterpolation[x] == MaximumFFTValue/2, x]
Can you help me please? How would you do that?
Integrate[DataInterpolation[x], {x, BegFWHM, EndFWHM}]
Which works fine with the interpolating function.
Thank you for your help!
Mathieu