How to divide a signal into groups of peaks
39 views (last 30 days)
Show older comments
Hi, i have an audio sample , which im trying to divide into groups of signals , for energy and time analysis . here is the signal :
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/152135/image.png)
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
4 Comments
Luis Gomes Perini
on 13 Jul 2015
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
Accepted Answer
Star Strider
on 13 Jul 2015
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/177685/image.png)
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.
2 Comments
Star Strider
on 13 Jul 2015
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.
As always, my pleasure. Biomedical signal processing is of particular interest to me, so I learn something from every new problem.
More Answers (1)
Image Analyst
on 13 Jul 2015
Edited: Image Analyst
on 26 Dec 2020
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;
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/177687/image.jpeg)
8 Comments
Image Analyst
on 13 Feb 2019
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 Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!