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

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

Answer 1

[Star Strider]

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

[steamrice]

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]

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]

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]

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]

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]

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:

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]

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]

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]

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]

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:

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.

[steamrice]

Thanks for the fast response! I will playing around and give you some updates!

[Star Strider]

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]

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]

‘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]

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]

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]

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]

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]

I see! Thanks!

[Star Strider]

As always, my pleasure!

I will help as I can.

[steamrice]

Thanks for your time! Do you have suggestion for the trigger siganls from your pervious experience?

[Star Strider]

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]

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]

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]

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]

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]

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]

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]

That's what I was thinking and planning to do. Is that something closer to what you meant?

credits to: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0036407

[Star Strider]

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]

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]

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]

I see! I am playing things around again and will ask you for more questions! Thank you!

[Star Strider]

As always, my pleasure!

[steamrice]

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]

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]

@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]

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]

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]

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..

[Star Strider]

‘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]

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]

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]

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?)

[steamrice]

What do you think of this method?

https://www.mathworks.com/help/signal/ug/formant-estimation-with-lpc-coefficients.html

[Star Strider]

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]

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]

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]

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

[Star Strider]

Yes!

[steamrice]

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]

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]

Another thing if you don't mind, I was trying the Ipc method:

https://www.mathworks.com/help/signal/ug/formant-estimation-with-lpc-coefficients.html

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);

[steamrice]

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

[Star Strider]

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.

Answer 2

[steamrice]

I attached my data sample if you have time to take a look? (The first two columns are voltage measuremnt whereas the third column is just a signal trigger)