Detect valleys of a signal

Question

I have raw values that are coming directly from a sensor (on the fly). These values are coming from a wearable sensor installed on a cows neck. When the cow ruminates (meaning that it chews food that brought back from the 2nd stomach of it), generates some kind of acceleration values that have some kind of specific amplitude and frequency (the chew duration is mostly between 800 msec to 1200 msec). A rumination period is between 30 to 120 seconds every time and the cow ruminates around 10 hours per day.

I do a quick calculation on these values and I end up with a graph like the bellow image.

I would like to be able to detect every time when the signal falls (like the highlighted peaks) as it is a representation of a rumination period. That fall in the graph represents the bolus event of the cow (meaning that the cow just shallow-ed the chewed food and brought back a new one [from the stomach] to chew again).

I tried to do it with a threshold approach (if the signal falls bellow a set threshold), but it fails to operate always correctly. This is because:

1) The signals intensities change slightly during the day and between different sensors eg. On the image, the negative peaks do not fall bellow the black horizontal line which is my threshold. As such I have false detection (not detecting them at all). In case I increase the threshold, I do wrong detections since a lot of other points in the signal might pass it.

2) When the signal is slightly noisy I do wrong detections or no detections at all. As an example, the 9th yellow point of the signal should be somehow detected as valley but it is well above the threshold. Of course it could be mistakenly detected as it has the same amplitude as some nearby points.

Does anybody know how I can detect such points (dips) in time?

Note: These data are coming real-time and they are not stored somewhere. As such I must somehow be able to detect dynamically the fall of the signal without having any knowledge of the future values of it (well I can have some small history of course - lets say last 200 values or something).

Extra info: These data values are actually coming from a 3-axis accelerator sensor. I calculate the magnitudes of them (sqrt(x^2 + y^2 + z^2)) which they deviate around 1g and then I get the Moving Standard Deviation of the last second (62.5 values). That is what is displayed on the graph above. The yellow points are at around 0.01g deviation.

Further: After being asked by some users, I also attach an image of the original signal (only a tiny part of it). With yellow are marked the areas of interest (on the magnitude). The rest 3 graphs are the X,Y,Z (legend exists on the top-right of the graph).

Also, I attach a hyperlink where you can download sample signals in CSV format. The format is as following: 1st, 2nd, 3rd columns are the X, Y and the Z axis accelerometer values. The 4th column is the magnitude of them (sqrt(x^2 + y^2 + z^2)) which you can also manualy calculate if you want so.

https://1drv.ms/f/s!Au5-DmkSvYQzgZlLo2KJ_UBAVeYDWw

There are files with 4 quite good signals. The Good3s.csv contains some small junk at the start of it (some 10000 first values can be discarded). The Good4e.csv has the same at the end. There are also 2 more signal files that are quite noisy and you can parse them if you like hardcore situations.

Answer 1

First, you'll probably have better luck posting this on dsp.stackexchange. That's a more specialized group that does stuff like this all the time.

In terms of your problem, here's a couple of options.

One is a machine learning approach. e.g. create a training set by taking a bunch of data and hand marking the points that are good versus bad (like you've done above). Then take an Artificial Neural Network, feed it the training data and teach it to recognize the pattern you want.

Another approach might be to take a Statistical process control approach. Lots of tools there for detecting when something gets out of line.

Another way might a floating threshhold. Compute the mean of the last X samples, and if the signal drops more than Y below that mean, call it bad. could also compute the standard deviation of the last X samples, and if the signal drops more than Z std devs then call it back.

or, potentially just high-pass filter the signal to remove the slowly drifting mean and then you can use a fixed threshhold.

Answer 2

I did not read carefully through the whole question as I don't have time now, but have you tried some form of robust peak detection? See e.g., https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html

Then you can set parameters such as the minimum (and maximum) distance, the prominence, the minimal height. If you look for minima just flip the sign of your signal.

Answer 3

So your basic task is to detect "quiet" parts where the amplitude of some vibrations drops below a threshold.

Can you make assumptions on the length of the "silent" phase until you want to detect it? After analyzing your data I see that most of the silent phases have a duration of ca. 200 samples and repeat every 3000 samples. Can you make assumptions on repetitions? (For example: after a "detection" there can't be another one for 2000 samples")

Below I've added my current code for GNU Octave (which IMHO is very great for prototyping algorithms, design filters, analyze data and so on)

# https://dsp.stackexchange.com/questions/51269/detect-valleys-of-a-signal
# https://en.wikipedia.org/wiki/Standard_deviation#Rapid_calculation_methods

graphics_toolkit fltk

pkg load image  ## for "nlfilter"
pkg load signal ## for "butter"

FS = 62.5;  # in the question "of the last second (62.5 values)"

fns = glob ("*.csv");

k = 1;

d = dlmread (fns{k}, ";");

## calculate magnitude;
mag = sqrt (sumsq (d(:, 1:3), 2)) - 1;

## remove the slow varying trend
[b, a] = butter (8, 0.01);
mag -= filter (b, a, mag);

## find parts where magitude is below +/- thres
thres = 0.04;
q = abs (mag) < thres;

## ensure, that the first and last sample never is "quiet"
q(1) = q(end) = false;

## starts of quiet part
qs = find(diff(q) == 1);

## ends of quiet part
qe = find(diff(q) == -1);

## length of the quiet parts
len = qe - qs;

## only consider parts longer than 2s
idx = find (len > 2 * FS);

plot (mag)

hold on
# draw line above detection
line ([qs(idx) qe(idx)].', repmat (thres, 2, numel (idx)), "color", "red", "linewidth", 5)
hold off

The red lines are the "detections"

PS: I guess you know https://en.wikipedia.org/wiki/Standard_deviation#Rapid_calculation_methods but just in case not