How can I detect shoulder peaks in a signal?

19 次查看（过去 30 天）
Tilemahos
Tilemahos 2018-2-25
评论： Blair 2024-2-13
Hello,
I am trying to detect subtle, "shoulder" peaks in a signal, which do not fit the usual peak definition and are not found by findpeaks. They definitely can be seen in the gradient, but so can other wider peaks that are not so prominent in the original signal. I'm guessing there is a specific threshold that makes these shoulder peaks stand out in the gradient, but I can't think of how to find it. Are there any ideas on how to detect those peaks?
Thanks in advance
P.S. Below is a signal with 2 shoulder peaks (left) and its gradient (right)

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2018-2-25
The gradient is definitely the way to go. If you have the Signal Processing Toolbox, use the findpeaks (link) function on the gradient vector. This will find the peaks. In order to find the ‘valleys’, negate your original signal and find the peaks on the negated (inverted) signal. Use the second output (the indices of the peaks/valleys), since they will be correct without transformation for both the original and negated signals. The findpeaks function has a number of name-value pair arguments that can make this much easier.
Example —
ds = gradient(signal)
[pks,pkidx] = findpeaks(ds); % Identify Peaks
[vls,vlidx] = findpeaks(-ds); % Identify Valleys
  10 个评论
显示 8更早的评论
Blair
Blair 2024-2-13
So, absolutely massive amount of time later, but I feel like there's something missing here. The location of the peaks in the gradient of the signal don't correspond to the peak locations in the signal itself.
Peaks in the signal correspond to points where the gradient is equal to zero (or at least changes sign, since it's unlikely your real signal equals exactly zero). In the case of shoulders with no prominence, like the original example, they approximately correspond to the closest (local) approach to zero without actually changing sign.
Here's a little MWE that replicates the data given in the original question:
% generate data
x = 0:1:300; % x coordinates
width = 20; % width for all gaussians
g1 = (300*exp(-((x-150).^2/(2*width^2)))); % generate 5 gaussians
g2 = (150*exp(-((x-200).^2/(2*width^2))));
g3 = (150*exp(-((x-100).^2/(2*width^2))));
g4 = (50*exp(-((x-250).^2/(2*width^2))));
g5 = (50*exp(-((x-50).^2/(2*width^2))));
y = g1 + g2 + g3 + g4 + g5; % sum them together
% display data and gradient
figure(1)
plot(x, y)
dataGradient = gradient(y);
figure(2)
plot(x, dataGradient)
% run peak finding on data
[pks, locs] = findpeaks(y,x);
disp("Peak locations from data:")
disp(locs)
% run peak finding on data gradient
[gradpks, gradlocs] = findpeaks(dataGradient,x);
disp("Peak locations from data gradient:")
disp(gradlocs)
% run peak finding on inverse data gradient
[igradpks, igradlocs] = findpeaks(-dataGradient,x);
disp("Peak locations from data gradient inverse:")
disp(igradlocs)
It generates peaks at 5 positions: 50, 100, 150, 200, 250. If you run the code, you'll notice that when findpeaks is run on the raw data, it only detects a peak at 150, since the shoulders have no prominence.
When you run findpeaks on the gradient of the data, it detects peaks at 31, 84, 134, 192 and 246.
When you run findpeaks on the inverse (negative) of the gradient of the data, it detects peaks at 54, 108, 166, 216 and 269.
The first thing you might notice is that across the two gradient methods, the peak that is least accurately determined is actually the large central peak at 150 (predicted as 134 from the normal gradient, and 166 from the inverse gradient). That's because the peak actually occurs at the point of inflection in the middle of the plot of the gradient, where it passes through zero, rather than any of the minima or maxima (peaks or valleys).
The next thing you might notice is that shoulders on the left side of the main peak are best picked up by the inverse gradient method (54 and 108 vs. 31 and 84), and shoulders on the right side of the main peak are best picked up by the normal gradient method (192 and 246 vs. 216 and 269). This relates to what we said earlier that shoulders with zero prominence approximately correspond to points of closest approach to zero in the gradient. On the left hand side of the main peak, the gradient is positive, so a "closest approach to zero" is a valley in the gradient (and therefore a peak in the inverse of the gradient). On the right hand side of the main peak, the gradient is negative, so a "closest approach to zero" is a peak in the gradient. This graphic might help a bit. The main central peak is marked in red, the shoulders on the left side of the main peak marked in green and the shoulders on the right side of the main peak marked in magenta:
It's important not to forget that the shoulder positions are only approximately equal to these points of "closest approach to zero". The true shoulder positions will be slightly one side or the other of these points. The reason being, the gradient at any one point is the sum of the gradients of each individual peak. So at the true shoulder peak position, the shoulder peak contribution to the gradient is zero, while the contribution from its neighbours is slightly positive/negative. When we move slightly off the center of the shoulder, such that the contribution from the shoulder is opposite in sign to its neighbour, that's where these points of "closest approach to zero" occur.
You might think that you could average the pairs of values for the gradient and inverse gradient peaks (i.e. average 31 and 54 to 42.5, 84 and 108 to 96, etc. such that your final list of peaks is 42.5, 96, 150, 204, 257.5. And you definitely could, the derivative of a gaussian is a pair of peaks (one peak, one valley) equally spaced either side of the peak center, so the average precisely gets you the peak position (for an isolated peak). Which is why taking the average gets you close to the real peak positions in this example. But I think the best way would be to probe the gradient OF the gradient, i.e. the curvature (second derivative). If I get time later I might see if I can put together an answer based on the curvature.
Anyway, this got much longer than first anticipated but I mostly just started writing stream of consciousness to try and reason through it. Sorry for reviving a 6 year old post!

请先登录，再进行评论。

更多回答（1 个）

Image Analyst
Image Analyst 2018-2-25
How about if the gradient is between 0.2 and 0.4?
shoulderLocations = gradientSignal > 0.2 & gradientSignal < 0.4;
This will give the logical index vector where there is one of those shoulder locations. You can then find the centroid of that with regionprops() if you want. You can then inspect the signal between all shoulder centroids to see if the signal is high or low.
  2 个评论
Tilemahos
Tilemahos 2018-2-25
Thanks for the answer. It seems like a good idea. As I'm working on a lot of different signals, is there a more general way of automatically estimating these 0.2 / 0.4 thresholds properly?
Image Analyst
Image Analyst 2018-2-25
You just have to look at a bunch of signals and try to determine what gradient most of the shoulders are at.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Descriptive Statistics 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by