How to get maxima that are not on a current datapoint?

Hi all,

I am trying to find the position of major infliction points in a signal. I know there are 2 main infliction points in the signal within the range I am looking at. I am writing a function that pulls these out, which in itself finds them (see below). However, these points are always located on one of the X-values in the dataset. I would like to find them on the actual position if one where to have a continuous line through the data.

I tried using the interpo1 function, on the signals using method pchip, creating 10 times as many points. But when I run that data through the routine, still the calculated maxima are located on the original datapoint location.

What could I be doing wrong?

function [peak] = JF_findpeaks(signatures, offset, startband, endband)

% how many spectra
  numspectra = size(signatures(:,1));

% calculate the second derivative of the signatures
% Set window 5, fit second order spline, return second order derivatives
  secondderivatives = JF_savitzky_derivative(signatures, 2, 5, 2);

% Subset to the range where we know the peaks are
  data = secondderivatives(:,startband-offset:endband-offset);

% create the X-values
  waves = startband:endband;

% now find the position of the two peaks in each signature that occurs in the data
  for n=1:numspectra
   [pksw,locsw] = findpeaks(data(n,:), waves, 'MinPeakProminence',0.00001);
  end

end

[Image Analyst]

Please attach your data and a screenshot of it plotted. Make it easier for us to help you. I have a hard time working blind.

Answer 1

[John D'Errico]

Um, you used pchip? !!!

Yes, I predict that was what you did wrong. Without looking at it more deeply to see if you have some other problem, just use of pchip is an assurance that you got exactly the result you did.

pchip has LOTs of good characteristics. But one "feature" of pchip is that on data that is non-monotonic, it will not introduce extreme between data points. Any extremum in the interpolated function will be exactly at the data points themselves.

x = 0:10;
y = rand(size(x));
xint = linspace(0,10,1000);
pp = pchip(x,y);
yint = ppval(pp,xint);
plot(x,y,'bo',xint,yint,'r-')

Ok. See that the interpolated function ALWAYS has a max o a min at any local max of min of the data. This is a feature of pchip.

Now, try it using a cubic spline instead. (What you get from interp1 with the 'spline' option.)

pps = spline(x,y);
yints = ppval(pps,xint);
plot(x,y,'bo',xint,yints,'r-')

See that this last curve has extrema that need not lie at a data point.

That does not make pchip a bad choice in general. On MANY functional forms, pchip is a far better choice of interpolant, when a classic spline will give complete crap. For example, try using a spline to interpolate a simple heaviside step. The spline will generate ringing (think Gibbs phenomena) like mad, whereas pchip will yield a nicely behaved interpolant, with NO ringing at all.

But you should never use a pchip interpolant when you have non-monotonic data, IF you are looking for an extremum between the data points. And that is exactly what you seem to be asking to find.

[Jelle]

This exactly the problem and its solution!

Thank you

The difference between pchip and spline made the difference. To think I have been trying to solve this for almost 2 hours. |tearing at hair.

In short, this is the functional part of the code, now it is a matter of tweaking the resolution and make it dependent on the input

resolution  = 50 * (endband-startband)
xints       = linspace(startband,endband,resolution);
pps         = spline(waves,data(n,:));
yints       = ppval(pps,xints);

plot(waves,data(n,:),'bo',xints,yints,'r-')

findpeaks(yints, xints, 'MinPeakProminence',0.00001)

[John D'Errico]

All you had to do was ask. OH! You did. ;-) Sorry about the loss of hair.

Seriously, this is a little recognized property of pchip, one that makes it tremendously useful SOME of the time. But sometimes, it just surprises people.

[Jelle]

:)

Answer 2

[Walter Roberson]

The interpolation methods available from interp1(), and also spline(), always go exactly through all points specified. To get a line "through" the point without specifying the point itself, you need to remove the data for the point itself from the interpolation data.

Now if you just delete that one point of data and then ask to interpolate there, then the value will usually be lower than the known data: you have defined the point as being a peak, so the adjacent values are less, so the interpolation will usually be less. An exception could be if you used a spline interpolant, in which case the interpolated arc might happen to go high. But realistically it probably will not go higher -- not unless you have closely-spaced steep data. (The opposite argument applies for the minima: the interpolated points are going to tend to be less than the point.)

Another interpolation method would be to take a portion of the data before and after the inflection and to do a polynomial fit of even order, and use the polynomial to predict the point.

[John D'Errico]

I'd argue the problem is not the use of an interpolant, but specifically pchip in this case. pchip is not designed to allow an extreme value that lies between data points. And that is exactly the failure mode that was reported.