How to use findpeaks to find the maximum for a stimulus that contains multiple events?

2 次查看（过去 30 天）

steamrice 2020-2-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/504694-how-to-use-findpeaks-to-find-the-maximum-for-a-stimulus-that-contains-multiple-events

评论： Star Strider 2020-3-4

Hi, I am wondering if I have a signal that contains mulitple events. Is it possbile to locate the the max peak for each event (see the attached picture) by uaing findpeaks? Thanks in advance

采纳的回答

Star Strider 2020-2-10

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/504694-how-to-use-findpeaks-to-find-the-maximum-for-a-stimulus-that-contains-multiple-events#answer_414832

I would first use findchangepts (R2016a and later) to segment the signal, then use findpeaks on each segment.

55 个评论
显示 53更早的评论隐藏 53更早的评论

Star Strider 2020-2-10

在 MATLAB Online 中打开

As always, my pleasure!

My apologies for the delay — I had an errand to run.

The complete code (now that I have it working):

D = load('pilotexample.mat');
data = D.data;
EMG = data(:,1);
Ts = 0.5;
t = linspace(0, size(data,1), size(data,1)).'*Ts;
cpv = findchangepts(EMG, 'MaxNumChanges',24,'Statistic','rms', 'MinDistance',fix(numel(EMG)/500))
tic
EMG = data(:,1);
Ts = 0.5;
t = linspace(0, size(data,1), size(data,1)).'*Ts;
cpv = findchangepts(EMG, 'MaxNumChanges',24,'Statistic','rms', 'MinDistance',fix(numel(EMG)/500))
toc
figure
plot(t, data(:,1))
hold on
yl = ylim;
plot(t([cpv cpv].'), (ones(size(cpv))*ylim).', 'r')
hold off
grid
for k = 1:fix(numel(cpv)/2)
    segment{k} = [t(cpv(k):cpv(k+1)), EMG(cpv(k):cpv(k+1))];                  % Save Segments In Cell Array
end
figure
subplot(2,1,1)
plot(segment{1}(:,1), segment{1}(:,2))
grid
title('Segment 1')
subplot(2,1,2)
plot(segment{10}(:,1), segment{10}(:,2))
grid
title('Segment 10')

with the segmented data plot:

How to use findpeaks to find the maximum for a stimulus that contains multiple events - 2020 02 10.png

and:

How to use findpeaks to find the maximum for a stimulus that contains multiple events (2) - 2020 02 10.png

This gives you a way of saving the segments and then addressing them as I do in figure(2). The code creating figure(1) plots the change points as vertical red lines.

Yes, you would have to use findpeaks on each ‘segment’ (as I call them here) to get the maximia. You might be able to do that in the loop that creates the individual ‘segment’ cell array. (I suggest saving the ‘segment’ array to a .mat file if you want to use it again. It would save the time of re-segmenting it.)

It takes half as long with 24 change points as with 20, likely because findchangepts does not have to ‘decide’ which 20 to report in its output. Using 24 also gets all of them. (Unfortunately, it may not be possible to completely automate a way to correctly determining the number of change points.)

Star Strider 2020-2-11

在 MATLAB Online 中打开

I was thinking of something like this:

for k = 1:fix(numel(cpv)/2)
    segment{k} = [t(cpv(k):cpv(k+1)), EMG(cpv(k):cpv(k+1))];                  % Save Segments In Cell Array
    [pks{k},locs{k}] = findpeaks(EMG(cpv(k):cpv(k+1)), 'MinPeakHeight',sqrt(mean(EMG(cpv(k):cpv(k+1)).^2)));
end
figure
subplot(2,1,1)
plot(segment{1}(:,1), segment{1}(:,2))
hold on
plot(segment{1}(locs{1}), pks{1}, '^r')
hold off
grid
title('Segment 1')
subplot(2,1,2)
plot(segment{10}(:,1), segment{10}(:,2))
hold on
plot(segment{10}(locs{10}), pks{10}, '^r')
hold off
grid
title('Segment 10')

doing each findpeaks call as the ‘segment’ arrays are created, then producing this plot:

How to use findpeaks to find the maximum for a stimulus that contains multiple events (3) - 2020 02 10.png

One of the peaks in ‘segment{1}’ definitely appears to be higher than 0.1. Is that a problem? (Note that if there is a strict limit on the amplitudes that you expect, the apparent ‘overshoot’ could be due to a sampling artifact. In my own studies, a value that occurs at the same time as a sampling time in a rapidly-changing signal can be distorted by the sampling process. I have no idea why that occurs, I only know that it does.)

You can extend the subplot figure to include all of the ‘segment’ vectors. Since there are 12 in this record, having 6 rows of subplots and 2 columns might work.

I used the RMS value of each EMG ‘segment’ vector for 'MinPeakHeight', since the amplitudes vary significantly in the EMG between the ‘segment’ vectors, and that appeared to be a relatively reasonable and adaptive approach. Use multiples of the RMS value to get different results, or your own value for 'MinPeakHeight' and other name-value pairs that findpeaks offers. The amplitudes appear to having been returned correctly here.

Star Strider 2020-2-12

在 MATLAB Online 中打开

‘I am not too sure how did you come up with fix(numel(EMG)/500). Why did you divide by 500?’

I adapted it to the length of the record. Use any value you want.

‘Could I also think of those segment as 'onset' of the EMG signals (as the stimulus occurs)?’

Yes. The first change point in each segment would work for that purpose. (I initially thought the trigger channel (third column) might be useful. However the third column of your .mat file only appeared to have two triggering events, so I did not use it. Also, they did not correspond to the other data.)

In terms of speed, I would use 24 change points for this record. That takes half the time that 20 change points does.

I just now thought of a way to automate determining the number of segments (half the number of change points):

EMG = data(:,1);
Ts = 0.5;
t = linspace(0, size(data,1), size(data,1)).'*Ts;
env = envelope(EMG, 2500, 'rms');
[epks,elocs] = findpeaks(env, 'MinPeakProminence', 0.0005);
plot(t, EMG)
hold on
plot(t, env)
plot(t(elocs), epks, '^g')
hold off
grid
<REST OF THE CODE>

If you have the Signal Processing Toolbox (R2015b and later), you have the envelope function.

For this record, it returns 12 peaks, so it does count the correct number of EMG events. Experiment with that to see if it works correctly with your other records. It might be necessary to change the 'MinPeakProminence' value, however the 0.0005 value works with this record.

steamrice 2020-2-12

在 MATLAB Online 中打开

You have been so helpful!! I am just looking into more detail regarding the segment related to the number of the events.

(I initially thought the trigger channel (third column) might be useful. However the third column of your .mat file only appeared to have two triggering events, so I did not use it. Also, they did not correspond to the other data.)

The trigger signal that I sent was supposed to correlate when a stimuls happens. It's not so useful at the moment as I pogrammed it only 'logic high' when the stimuls happens 'during that stimlus period' using Matlab. So if you see the first picture I posted you would see there is first 4 events and pause and the last 4 events (Total of 8 events). The trigger signal will be high for the first 4 and low during the pause and high again for the last 4. This is because I put the trigger siganl within the events in Matlab. I am hoping to change it in such a way that the trigger signal only high when there is an event from the EMG. (Like a series of logic high and logic low for each event instead of each 'period'). I do not know whether it will be possbile....(This may affect the way I analyze the data).

For

env = envelope(EMG, 2500, 'rms');

The sliding window of length wl samples (2500) in this case. How should we determine this number? Is it dependent on the frequency of the signal? (I think it was 2000Hz)

Thanks again!

steamrice 2020-3-4

在 MATLAB Online 中打开

Thanks for your repsosne as always!

Baesd on the definiiton: "In classifying speech sounds, the lowest strong frequency band is termed the first formant. The next highest is termed the second formant, and so on."

I have thought about using spectrogram as well to have a plot in terms of frequency vs time. But I am not too sure how to detect the frequency that matchs to what I am looking for.

what I am thinking and having is to load the audio, then plot it in terms of dB vs freuency. Then I am just wondering how to select the appropriate "formant". This is what I am trying to acheive:

Creadits to:

https://home.cc.umanitoba.ca/~krussll/phonetics/acoustic/formants.html

I have read a post and thats what I have:

[audio_file, Fs] = audioread("testing.wav");
L = 10000;
Ts = 1/Fs;                                                      
Fn = Fs/2;                                                    
FT_af = fft(audio_file)/L;                                      
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                             
Iv = 1:numel(Fv);                                               
[PksL,LocsL] = findpeaks(abs(FT_af(Iv,1))*2, 'MinPeakHeight',0.006, 'MinPeakDistance',20);
[PksR,LocsR] = findpeaks(abs(FT_af(Iv,2))*2, 'MinPeakHeight',0.006, 'MinPeakDistance',20);
figure(1);
semilogx(Fv, 20*log10(abs(FT_af(Iv,2))*2)), xlabel('Frequency, Hz'),ylabel('magnitude');
hold on
plot(Fv(LocsR), PksR, '^r', 'MarkerFaceColor','r')
hold off
xlabel('Frequency (Hz)')
ylabel('Amplitude')

However, the scale is not appropriate yet and I haven't got the the formant properly..

steamrice 2020-3-4

在 MATLAB Online 中打开

Thanks for your resposne!!

I am wondering that is what is the y-axis represents for when we plot the audio in fft? We do have to use

semilogx(Fv, 20*log10(abs(FT_af(Iv,2))*2))

to convert it to dB right?

I would solve it by creating different filters (the firls function could be appropriate here), and then take the time-integrated output from all the filters that filtered the same signal and look for the greatest result for a specific time range of a signal. That result would match the filter for that formant.

Could you direct me on this? Creating different filiters (like low pass and high pass ) and select the frequency range that we are interested in using firls? I am unsure the what do you mean by take the time-integrated output ...

Since I just want to compute the first two formants, will findpeak do the job by looking for the first two lowest strong band in the dB vs Hz plot?

Sorry for asking many questions. Your suggestions have really pushed me to exlore different options and helped me to formule possible solutions!

steamrice 2020-3-4

My email comes up with a delay! But thanks for helping it throughout!

Star Strider 2020-3-4

I have not used that function or explored that example. The roots function returns the roots of a polynomial, so the argument must be a vector of polynomial coefficients.

请先登录，再进行评论。