How to divide a signal into groups of peaks

Hi, i have an audio sample , which im trying to divide into groups of signals , for energy and time analysis . here is the signal :

what i'm trying to do is, to get the start and the end indexes of this 4 chunks. right now i dont have a clue how to get this done . Thank you

[Azzi Abdelmalek]

What are your criterion? What do you mean by start and end?

[cob]

and the same for the next two

[Luis Gomes Perini]

Hi,

your question is how to develop an algorithm that captures such peaks, or how to determine the start/end of peaks ?

Cheers, Luis Gomes Perini

[cob]

the main purpose is to count how many chunks/peaks, with similar energy and duration between them , in a specific length , for now i want to determine the start index and the end index of such peak/chunk , which will help me in the future analysis . Of course its not a place to ask for algorithm development , its on my own.

Thank you in advance.

Answer 1

[Star Strider]

This is likely the most difficult signal you’ve presented. Probably the best way to produce a derived signal that would allow you to isolate the parts of your signal of interest is to use a bandpass filter and then a Savitzky-Golay filter.

The code:

cs = load('cob signal.mat');
S1  = cs.coarse_d(:,1);
S2  = cs.coarse_d(:,2);
Fs = cs.fs_d;
Fn = Fs/2;
Ts = 1/Fs;
T  = [0:length(S1)-1]*Ts;

Wp = [950 1000]/Fn;
Ws = [0.5 1/0.5].*Wp;
Rp =  1;
Rs = 30;
[n,Wn] = buttord(Wp, Ws, Rp, Rs);
[b, a] = butter(n, Wn);
[sos, g] = tf2sos(b, a);
S1F = filtfilt(sos, g, S1);

SGS1 = sgolayfilt(abs(S1F), 1, 901);
figure(1)
plot(T, S1)
hold on
% plot(T, S1F)
plot(T, SGS1, 'LineWidth',1.3)
hold off
grid
% axis([0.5  1.5    ylim])

The first ‘burst’ of the original signal (with the Savitzky-Golay filter of the bandpass filtered signal superimposed) is:

You could simply threshold the sgolayfilt output, or use findpeaks to determine the locations of each ‘burst’. I didn’t comment-document the code, because it’s similar to my previous filter designs (documented), the only novelty being the addition of sgolayfilt.

[cob]

WOW , thank you very much , i've asked only one little thing , and you just solved the whole problem , again thank you . Can you please explain why do you use the commands ts2sos , and filtfilt and how you determined Wp valuse ?

again thank you very much

[Star Strider]

My pleasure!

One item I left out was the fft code, which is how I arrived at the filter design Wp values. I include it here:

FS1 = fft(S1)/length(S1);
Fv  = linspace(0, 1, fix(length(FS1)/2)+1)*Fn;
Ix  = 1:length(Fv);

figure(2)
plot(Fv, abs(FS1(Ix)))
grid
axis([0  5E3    ylim])

I chose the frequency band that seemed to have the most energy, in order to make it easier for sgolayfilt. It took a bit of experimentation with the bandpass and Savitzky-Golay filters to get a result I liked well enough to post.

Transfer-function designs can be unstable and produce anomalous results, so I use the second-order section representation of a filter (therefore tf2sos) to provide stability. I usually use freqz to assess filter stability of the transfer-function and second-order-section implementations, just to be sure they work. (MATLAB suggests this in its documentation, and one of the filter design tools does it automatically.)

The filtfilt function (as opposed to filter) is phase-neutral (the sort of behaviour that favours the Bessel design in hardware implementations). The filtfilt function renders all discrete filters phase-neutral.

My detailed explanation of filter design is here.

As always, my pleasure. Biomedical signal processing is of particular interest to me, so I learn something from every new problem.

Answer 2

[Image Analyst]

Well, being an image analyst, I came up with a different approach. You can get a logical array defining what's a signal and what's the quiet parts with just two lines of code, one to threshold and one to do a morphological filter. Here's a full blown demo:

% Setup - read in the signal and plot it.
s=load('signal.mat')
signal = s.coarse_d;
subplot(3,1,1);
plot(signal, '-', 'Color', [0,.5,0]);
grid on;

% MAIN CODE IS THE FOLLOWING TWO LINES OF CODE.
% Threshold the signal
lowSignal = abs(signal) < 0.1;
% Remove stretches less than 10,000 elements long.
% And invert it to get the high signal areas instead of the quiet areas.
lowSignal = ~bwareaopen(lowSignal, 10000);

% Now we're done.  Plot the results.
subplot(3,1,2);
plot(lowSignal, 'b-', 'LineWidth', 2);
grid on;

% Plot both on the same graph now.
subplot(3,1,3);
plot(signal, '-', 'Color', [0,.5,0]);
hold on;
plot(lowSignal, 'b-', 'LineWidth', 2);
grid on;

[cob]

Thank you , its another approach , it seems that it also doing the right thing . can you please explain what does the ~bwareaopen command do ? again thank you.

[Image Analyst]

If you have a binary signal (true or false), then bwareaopen() lets you say how big are the "connected components" that you want to keep. For example in the "signal" part, the signal goes up and down very rapidly so there are hundreds of small chunks that are more than the threshold with hundreds of small chunks in between them where the signal value is low. We don't want all those low amplitude parts that occur in the main pulses of the signal. We only want to keep a quiet part if it's more than, say, 10,000 elements long. If it's only like 20 or 50 long, we want to throw those away - we don't want them because they're in the oscillating main signal pulses, not in the long stretch of quiet in between the main pulses. So calling that tells it to throw away short quiet chunks and only keep the really long quiet chunks. Does that explain it?

[Star Strider]

That seems analogous to findpeaks with the appropriate 'MinPeakDistance' value set to ignore peaks within a set range of the designated x-variable (index or defined).

I have the Image Processing Toolbox, but haven’t used it to process 2D data. I may be overlooking significant functionality.

[cob]

Yes, its explain it very well , just one more question , how you decided the value of 10,000 elements long ?

[Image Analyst]

I just looked at the length of the gaps and picked something that I thought would be way longer than any short segment in the pulse, but not so long as to get rid of the long quiet stretches if your pulses happened to be closer together than what you showed (so it does not fail in some other data set).

[MA-Winlab]

after identifying the non-silent parts, how can you put their data points in a colomn of a matrix. i.e., can you create a matrix that has the data points of each identified portion as a colomn. For example if there are 10 splikes (non silent portions), then how to store them (as data points) in 10 colomns in a matrix.

Thank you

[Image Analyst]

You'd have to extract each segment using indexing then put them into a cell array because each segment won't have the same number of elements and thus cannot be put into a double array (unless the ends were padded with zeros or nans or something.

See the FAQ: https://matlab.wikia.com/wiki/FAQ#What_is_a_cell_array.3F

[Surabhi KS]

You can make a class object.