I would first use findchangepts (R2016a and later) to segment the signal, then use findpeaks on each segment.
How to use findpeaks to find the maximum for a stimulus that contains multiple events?
9 次查看（过去 30 天）
显示 更早的评论
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
55 个评论
steamrice
2020-2-10
Thanks for the response! That may work! I am just looking at an example
findchangepts(t,'MaxNumChanges',8)
But should I segment a way such that it seperate, says the 8 events, from the inital point to the end point of each event? (would it be possible?)
Star Strider
2020-2-10
My pleasure!
Without your signal to experiment with, it is not possible for me to determine what the best option would be. Consider using other name-value pairs as well, such as 'Statistic' and others so that findchangepts will do what you want more easily. There are similar functions that you may consider experimenting with, such as cusum and ischange (R2017b and later), both with different calling conventions and syntax.
Star Strider
2020-2-10
编辑：Star Strider
2020-2-10
Using:
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',20,'Statistic','rms', 'MinDistance',fix(numel(EMG)/500))
the findchangepts output for the first column of data are:
cpv =
6408
8492
13240
15902
20901
22658
27312
28504
41208
42055
47777
49485
54219
56000
61036
61889
66853
67899
Elapsed time is 555.593266 seconds.
corresponding to:
I will re-run this with 22 change points (for the other 2 spikes).
steamrice
2020-2-10
Yes! That's another point I wanted to point out. findchangepts function took me a while to run and I had to stop the program for a few times.
steamrice
2020-2-10
编辑：steamrice
2020-2-10
Thanks! I am testing this as well. Just a few questions while I am waiting for it to run. How do you determine the number of those changing lines (the 18 of them)? In addiiton, if I want to check the peaks between each section, do I have to do each section one by one? Using findpeaks on each segment and enter the boundary as you indicated I am thinking? Also, is there a reason why it would take about 10 mins to run? (is it because of the number of the data points?)
Star Strider
2020-2-10
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:
and:
%20-%202020%2002%2010.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.)
steamrice
2020-2-11
Thanks for your response, that is really helpful! I am wondering the difference in the x-aixs for the segment 1 and segment 10 graphs. What are they representing?
Star Strider
2020-2-11
As always, my pleasure!
They are representing the absolute times. (I opted for that in order to preserve them in the event that was necessary.) If you want to represent the relative times (so that each starts at 0), subtract the initial time from the rest of the vector.
For example:
plot(segment{10}(:,1)-segment({10}(1,1), segment{10}(:,2))
will plot ‘segment{10}’ beginnining at 0. That applies to all the others. Using relative times might be easier if you are using findpeaks. It would also standardise them if you are interested in exploring morphological similarities.
steamrice
2020-2-11
That helps! I was just testing and learning how to find peaks for the next steps. But I am a bit confused. I tried : (say I wanted to find peak of the first segment)
findpeaks(segment{1) (:,2),'MinPeakHeight',0.1 )
The plots showed the segment 1 and give peak of the amplidue that is higher than 0.1. So for segment 1, I should use segment{1) (:,2) right instead of segment{1) (:,1)? I am a bit confused here.
Also, what did you by "You might be able to do that in the loop that creates the individual ‘segment’ cell array." for the findpeaks?
Currently I have to kinda plot the graph to see its ampliude to set the MinPeakHeight. Is it the standard way to do it? I am wondering would there be a more efficient way to find that amplidue if the signals are differnet
Star Strider
2020-2-11
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:
%20-%202020%2002%2010.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-11
As always, my pleasure!
I wish it could be faster, however findchangepts — for all its wonderful abilities — still takes 4½ minutes for this data set.
steamrice
2020-2-12
Thanks for the response! I think this way would better help me in terms of getting the amplitude because potentially the signal could vary signifcantly. Having manually set the ampltude (like 0.1 for segement would not be efficient!)
I was looking back to the 'segement', for this line:
cpv = findchangepts(EMG, 'MaxNumChanges',20,'Statistic','rms', 'MinDistance',fix(numel(EMG)/500))
Q1: I am not too sure how did you come up with fix(numel(EMG)/500). Why did you divide by 500?
Q2: Could I also think of those segment as 'onset' of the EMG signals (as the stimulus occurs)?
Star Strider
2020-2-12
‘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
Wow! That seems really saving a lot of processing time! Am experiening with it and will see how it goes! Thanks so much!
Star Strider
2020-2-12
As always, my pleasure!
I wish I had thought of it sooner! It is still necessary to use findchangepts, however with the correct number of change points for each record (2*numel(epks)) that should be easier.
steamrice
2020-2-12
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!
Star Strider
2020-2-12
As always, my pleasure!
The window length for envelope is what I determined experimentally to provide a reasonable result. Like most everything in signal processing with physiological signals, the parameter choices are chosen by experimenting and then using the parameter that appears to give the best result. There is nothing inherently ‘magic’ about 2500, so experiment with other values. You may find one that works better.
steamrice
2020-2-12
Thanks for your time! Do you have suggestion for the trigger siganls from your pervious experience?
Star Strider
2020-2-12
In the one multiple-animal study that I did, the stimuli were immediately followed by a muscle contraction, so the stimuli (recorded on a different channel) could be used as the start of the subsequent EMG record, and the end of the previous record. That made segmenting the EMG signals much easier. That does not appear to be the same for your data. Also, I thought ‘trigger’ = ‘stimulus’. It does not with your data.
No worries about that — you designed your study diifferently, and apparently with a different objective, than I did mine.
steamrice
2020-2-12
That's great! I will be using multiple channels for different mucle contration! But I was planning to have a trigger siganl such that it can indicate when events happen. I think I will try things around and ask you for more questions!!
Star Strider
2020-2-12
My experience is a bit limited, and I did my animal studies in the mid-1990s. I will nevertheless offer what assistance I can.
steamrice
2020-2-13
WOW. That's impressive in the 1990s!
That might be a liitle be unrelated to what I orginally asked, but I am getting struck and hope to get some ideas from you. For the EMG measurement, I have an external device connected to the matlab that can send the logic high and logic low as a marker corresponding to the stimulus. But I think it does not work too well. What I was thinking was hoping to get the trigger signal as event marker to tell when the stimulus happens and when the event was being carried out. (event happens as logic high and stimulus being carried out as logic low). This is because the elasped time running for the Matlab code is not matching the stimulus. When you were woring with EMG measruemnt, did you have a way/were you interested in finding/ getting the correct stimulus repsonese's time ?
Star Strider
2020-2-13
I recorded the stimulus on a separate channel, so I had simultaneous stimulus and response (EMG) recotrds and output of the force-displacement transducer (three channels, and a separate simultaneous time channel). If you have that capability, most of your problems in finding the beginning (if not also the end) of a particular stimulus are solved. I certainly recommend that sort of instrumentation setup if you have those options.
steamrice
2020-2-13
If I understand it right, isn't it similar to what I have (mulitple channels for EMG signals and one channel you see from the picture (the square wave) as stimulus response?). I am still in a design stage ish planning on how to get accurate data in which I hope the data analysis steps will be eailer...
Star Strider
2020-2-13
I believe so, however the ‘trigger’ does not seem to be an actual stimulus. That may be the difference.
I was simply describing the setup I used.
steamrice
2020-2-13
That's what I was thinking and planning to do. Is that something closer to what you meant?
Star Strider
2020-2-13
Essentially, yes. I used the stimulus to ‘gate’ the responses in my analyses. I did not use the sort of record depicted in the image, instead using the short square-wave (impulse) stinulus onset for each EMG response. I considered that the previous response ended when the subsequent stimulus began, since some stimulus frequencies were reasonably high. (I was gathering data to create a model reference controller for muscle stimulation, so I used various frequencies and amplitudes.)
steamrice
2020-2-14
Did you implment this method using MATLAB? Or did you do it manually each time? Like I am unsure how would this would work if the stimulus are changed?
Star Strider
2020-2-14
The instrumentation was a dedicated lab instrumentation system hooked up to an Apple machine current in the mid-1990s. (I don't remember the details.) It sampled the signals at 1 kHz, and produced .csv text files that I could easily read with MATLAB. I used MATLAB for all the analyses.
steamrice
2020-2-14
I see! I am playing things around again and will ask you for more questions! Thank you!
steamrice
2020-2-18
Good day!
I was playing around my setup and realized that the problem I have now is more related to the timing setup! I am worried about is the Matlab script that I have for running the experiment. I have a "for loop" that has:
getting the audio file>saving the audio file (by calling a function)>playing the audio (by calling another function). I'd like to have a 'trigger signal' such that the stimulus plays (1s) and record the onset of the person (1s).
But I just checked the timing of the 'for loop' I have by using tic and toc. It says 'Elapsed time' is about 3.3/3.4s for the loop to run. I am thinking that the unmatching timing comes from the MATLAB script... I am unsure of how to fix this at this point since each command seems necessary (I want to keep what I have: opening the audio, saving the audio, playing the audio with defined frequency...)
I have attached a picture trying to show you what my problem is.
Thank you again for your time! 
Star Strider
2020-2-18
I am not certain that I understand what you want to do. I am also not certain that MATLAB is entirely appropriate for this. I have never done something like this with MATLAB, only with dedicated instrumentation, using MATLAB to analyse the data offline later. I have no experience with using MATLAB this way.
steamrice
2020-3-4
@Star Strider. I just noticed that your won't take any email question. Is it okay to ask you my new question here? Please
steamrice
2020-3-4
Can you direct on how to find the formants from a frequency plot (after loading the sound file and using fft)? Is it possbile?
Star Strider
2020-3-4
I am not certain what a ‘formant’ is with enough understanding to detect it. If you are interested in time-frequency plots, the spectrogram function is likely appropriate.
steamrice
2020-3-4
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:
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..
Star Strider
2020-3-4
‘But I am not too sure how to detect the frequency that matchs to what I am looking for.’
Nor am I. This is quite remote from my areas of expertise.
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.
Wavelets could be another way to solve this, however I am not sufficiently familiar with wavelets to craft a specific approach.
steamrice
2020-3-4
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!
Star Strider
2020-3-4
As always, my pleasure!
Doing this transformation:
20*log10(abs(FT_af(Iv,2))*2)
converts ‘FT_af’ to dB.
See the documentation for the firls function to understand how to create the filters. After you design them (using firls or designfilt), use freqz to be sure the filter works as you intend it to work.
The time integration simply uses trapz (or cumtrapz depending on how you formulate the problem) to integrate the outputs of the various filters. The filter with the greatest output is likely the one you want. As I mentioned, this is very far from my areas of expertise, so there may be much more efficient ways of doing this.
If you want to use findpeaks, it will be necessary for you to take the Fourier transforms of specific regions of the time-domain vector that represent the sounds you want to examine. Taking the Fourier transform of the entire record will lllikely give you no useful information.
steamrice
2020-3-4
Thanks for the response! But based on my pervious graph, I am wondering why my peaks are not matching to the plot? (the red peaks are not on the curve?)
Star Strider
2020-3-4
Plot the peaks in dB:
plot(Fv(LocsR), 20*log10(PksR), '^r', 'MarkerFaceColor','r')
They should now be on the curve.
Also, for your purposes, it might be best to use a linear frequency axis rather than a logarithmic one. It may make the information easier to see.
steamrice
2020-3-4
Also, for your purposes, it might be best to use a linear frequency axis rather than a logarithmic one. It may make the information easier to see.
Do you mean fft(PksR) in my case?
Star Strider
2020-3-4
No. The ‘PksR’ vector are the peaks of the absolute value of the fft, so just plot them using the decibel transformation. Nothing else is necessary.
The linear frequency axis is simply a suggestion in order to make the spectrum easier to visualise.
steamrice
2020-3-4
I see what you meant! You meant go to the poperty in the graph and changed it to linear rather than log? It does look better now
steamrice
2020-3-4
编辑：steamrice
2020-3-4
I think I am still trying to find out a way to better do this. If I just want to find the first two higest peaks from the fft plo in a certain freqency range.(does that make much sense to you?) Is it possbile to just use findpeaks() to do so? You taguht me how find peaks worked and I tried to use in this case
I am thinking:
[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);
[PksR,LocsR] = findpeaks(abs(FT_af (:,1)), 'MinPeakHeight',0.006, 'MinPeakDistance',20);
plot(Fv(LocsR), abs(FT_af (:,1)), '^r', 'MarkerFaceColor','r');
But I recived the following error:
Index exceeds the number of array elements (5001).
Error in
plot(Fv(LocsR), abs(FT_af(:,1)), '^r', 'MarkerFaceColor','r');
Star Strider
2020-3-4
Using findpeaks makes sense. I am not certain what peaks you are looking for.
Experiment with 'MinPeakProminence' to detect them. It may be more reliable than 'MinPeakHeight' and 'MinPeakDistance'.
steamrice
2020-3-4
Another thing if you don't mind, I was trying the Ipc method:
load testing
segmentlen = 100;
noverlap = 90;
NFFT = 128;
spectrogram(FT_af (:,1),segmentlen,noverlap,NFFT,Fs,'yaxis')
title('Signal Spectrogram')
dt = 1/Fs;
I0 = round(0.1/dt);
Iend = round(0.25/dt);
x = (I0:Iend);
x1 = x.*hamming(length(x));
preemph = [1 0.63];
x1 = filter(1,preemph,x1);
A = lpc(x1,6);
rts = roots(A);
rts = rts(imag(rts)>=0);
angz = atan2(imag(rts),real(rts));
[frqs,indices] = sort(angz.*(Fs/(2*pi)));
bw = -1/2*(Fs/(2*pi))*log(abs(rts(indices)));
nn = 1;
for kk = 1:length(frqs)
if (frqs(kk) > 90 && bw(kk) <400)
formants(nn) = frqs(kk);
nn = nn+1;
end
end
formants
I am having a trobule with A = lpc(x1,6);
Error using roots (line 23)
Input must be a vector.
Error in testingmethod(line 77)
rts = roots(A);
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Measurements and Feature Extraction 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改，无法完成操作。请重新加载页面以查看其更新后的状态。
Translated by 
您也可以从以下列表中选择网站：
如何获得最佳网站性能
选择中国网站（中文或英文）以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)