How do I find all peaks and troughs of tidal data?

I'm working with some ocean tide data that's structured like this:

$data = array('date' => array('time' => array('predicted','observed')));

Here's a sample of real data that I'm using: http://pastebin.com/raw.php?i=bRc2rmpG

And this is my attempt at finding the high/low values: http://pastebin.com/8PS1frc0

Current issues with my code:

• When the readings fluctuate (as seen in the 11/14/2010=>11:30:00 to 11/14/2010=>11:54:00 span in the sample data), it creates a "wobble" in the direction logic. This creates an erroneous Peak and Trough. How can I avoid/correct this?

Note: My method is very "ad-hoc".. I assum开发者_如何转开发ed I wouldn't need any awesome math stuff since I'm not trying to find any averages, approximations, or future estimations. I'd really appreciate a code example of a better method, even if it means throwing away the code I've written so far.

I've had to perform similar tasks on a noisy physiological data. In my opinion, you have a signal conditioning problem. Here is a process that worked for me.

1. Convert your time values to seconds, i.e. (HH*3600)+(MM*60)+(SS), to generate a numeric "X" value.
2. Smooth the resulting X and Y arrays with a sliding window, say 10 points in width. You might also consider filtering data with redundant and/or bogus timestamps in this step.
3. Perform an indication phase detection by comparing the smoothed Y[1] and Y[0]. Similar to the post above, if (Y[1] > Y[0]), you may assume the data are climbing to a peak. If (Y[1] < Y[0]), you may assume the data are descending to a trough.
4. Once you know the initial phase, peak and trough detection may be performed as described above: if Y[i] > Y[i+1] and Y[i] < Y[i-1], you have encountered a peak.
5. You can estimate the peak/trough time by projecting the smoothed X value back to the original X data by considering the sliding window size (in order to compensate for "signal lag induced" by the sliding window). The resulting time value (in seconds) can then be converted back to an HH:MM:SS format for reporting.

You're looking for local minima and maxima, I presume? That's really easy to do:

<?php

$data = array(1, 9, 4, 5, 6, 9, 9, 1);

function minima($data, $radius = 2)
{
  $minima = array();

  for ($i = 0; $i < count($data); $i += $radius)
  {
    $minima[] = min(array_slice($data, $i, $radius));
  }

  return $minima;
}

function maxima($data, $radius = 2)
{
  $maxima = array();

  for ($i = 0; $i < count($data); $i += $radius)
  {
    $maxima[] = max(array_slice($data, $i, $radius));
  }

  return $maxima;
}

print_r(minima($data));
print_r(maxima($data));

?>

You just have to specify a radius of search, and it will give you back an array of local minima and maxima of the data. It works in a simple way: it cuts the array into segments of length $radius and finds the minimum of that segment. This process is repeated for the whole set of data.

Be careful with the radius: usually, you want to select the radius to be the average distance from peak to trough of the data, but you will have to find that manually. It is defaulted to 2, and that will only search for minima/maxima within a radius of 2, which will probably give false positives with your set of data. Select the radius wisely.

You'll have to hack it into your script, but that shouldn't be too hard at all.

I haven't read it in detail, but your approach seems very ad-hoc. A more correct way would probably be to fit it to a function

 f(A,B,w,p;t)=Asin(wt+p)+B 

using a method such as non-linear least squares (which unfortunately has to be solved using an iterative method). Looking at your sample data, it seems like it would be a good fit. When you have calculated w and p, it's easy to locate the peaks and valleys by just taking the time derivative of the function and solving for zero:

t = (pi(1+2n)-2p)/w

But I suppose, that if your code really does what you want, there's no use to complicate things. Stop second-guessing yourself. :)

A problem is I think that the observations are observations and can contain small errors. That at least needs to be accounted for. For example:

• Only change direction if at least the next 2 entries are also in the same direction.

• Don't let decisions be made by data on a too small difference. Throw away insignificant numbers. It will be a lot better probably when you say $error = 0.10; and change your conditions to if $previous - $error > $current etcetera.

How accurate does the peak/valley detection have to be? If you just need to find the exact record where a peak or valley occurs, isn't it enough to check for inflection points?

e.g. considering a record at position 'i', if record[i-1] and record[i+1] are both "higher" than record[i], you've got a valley. and if record[i-1] and record[i+1] are both lower than record[i], you've got a peak. As long as your sampling rate is faster than the tide changes (look up Nyquist frequency), that process should get you your data's peaks/troughs.

If you need to generate a graph from this and try to extrapolate more accurate time points for the peaks/troughs, then you're in for more work.

One way may be to define an absolute or relative deviation past which you classify further peaks/troughs as new ones rather than fluctuations around an existing peak/trough.

Currently, $direction determines whether you are finding a peak or trough, so instead of transiting to the other state (finding the trough or peak) once the derivative changes in sign, you can consider changing the state only when the deviation from the current peak/trough is "large" enough.

Given that you should never see two max or 2 min in less than about 12 hours, a simple solution would be to use a sliding windows of 3-5 hr or so and find the max and min. If it ends up being the in the first or last 30 min, ignore it.

As an example, given the following data:

1 2 3 4 5 6 5 6 7 8 7 6 5 4 3 2 1 2

and a window of size 8, with the first and last 2 ignored and only looking a peeks you would see:

1 2 | 3 4 5 6 | 5 6,  max = 6, ignore = Y
2 3 | 4 5 6 5 | 6 7,  max = 7, ignore = Y
3 4 | 5 6 5 6 | 7 8,  max = 8, ignore = Y
4 5 | 6 5 6 7 | 8 7,  max = 8, ignore = Y
5 6 | 5 6 7 8 | 7 6,  max = 8, ignore = N
6 5 | 6 7 8 7 | 6 5,  max = 8, ignore = N
5 6 | 7 8 7 6 | 5 4,  max = 8, ignore = N
6 7 | 8 7 6 5 | 4 3,  max = 8, ignore = N
7 8 | 7 6 5 4 | 3 2,  max = 8, ignore = Y
8 7 | 6 5 4 3 | 2 1,  max = 8, ignore = Y
7 6 | 5 4 3 2 | 1 2,  max = 7, ignore = Y