How to separate waves from a signal

14 views (last 30 days)
I have the following signal on the image above:
I want to delete all the waves of small amplitudes (under 0.1) and retrace the signal so I can have this figure
  3 Comments
Image Analyst
Image Analyst on 30 Jun 2022
@Jeffrey Clark, your Comment actually looks like an "Answer" to me, not a Comment. Can you therefore move it down to the official "Answers" section? This "Comments" section, where you posted, is supposed to be where people ask the original poster to clarify the wording of the question, attach data, etc. Answers go below in the Answers section. A benefit is that down there (unlike here) you can earn "reputation points" if the original poster Accepts your Answer or if anyone clicks the "Vote" icon. Thanks in advance. 🙂
Abdelhakim Souidi
Abdelhakim Souidi on 1 Jul 2022
% signal
load('signal.mat');
%time array
load('time.mat');

Sign in to comment.

Accepted Answer

Image Analyst
Image Analyst on 1 Jul 2022
If you have the Image Processing Toolbox (used to find the center of the "silent" parts), try this:
% Demo by Image Analyst
% Initialization Steps.
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
clear; % Erase all existing variables. Or clearvars if you want.
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;
s1 = load('signal.mat')
s2 = load('time.mat')
y = s1.signal;
t = s2.time; % Don't use time as the name of a variable since it's a built-in function.
subplot(3, 1, 1);
plot(t, y, 'b-');
grid on;
[upper, lower] = envelope(y, 30, 'peak');
hold on;
plot(t, upper, 'r-', 'LineWidth', 2);
threshold = 0.01;
mask = upper < threshold;
% Throw out any sections narrower than, say, 50 elements.
mask = bwareafilt(mask, [50, inf]);
yline(threshold, 'Color', 'g', 'LineWidth', 2)
y(mask) = 0;
title('Signal plus envelopes and threshold', 'FontSize', fontSize)
subplot(3, 1, 2);
plot(t, mask, 'g-', 'LineWidth', 4)
grid on;
ylim([0, 1.5]);
title('Mask indicating where no signals are', 'FontSize', fontSize)
% Find the centroid (half way point between the pulses and the peak intensities
% because we want to erase the shorter pulses between those half way points.
propsFlat = regionprops(mask, 'Centroid');
xy = vertcat(propsFlat.Centroid);
midPoints = round([1; xy(:, 1); numel(mask)])
% Now scan across each pulse section and if the max value
% in that section is less than 0.5 erase the whole section.
yMasked = y; % Initialize.
for k = 1 : length(midPoints) - 1
index1 = midPoints(k);
index2 = midPoints(k+1);
if max(y(index1 : index2)) <= 0.5
yMasked(index1 : index2) = 0; % Erase this section.
end
end
subplot(3, 1, 3);
plot(t, yMasked, 'b-', 'LineWidth', 1)
grid on;
title('Masked Signals', 'FontSize', fontSize)
  3 Comments
Abdelhakim Souidi
Abdelhakim Souidi on 1 Jul 2022
@Image Analyst Last question if possible, if I want to do the inverse (focusing on the small amplitudes) how to do it ?
Image Analyst
Image Analyst on 1 Jul 2022
Change this line
if max(y(index1 : index2)) <= 0.5
to this line
if max(y(index1 : index2)) >= 0.5

Sign in to comment.

More Answers (2)

Image Analyst
Image Analyst on 30 Jun 2022
Not hard, but you forgot to attach your data. But it would go something like this
[upper, lower] = envelope(y);
mask = upper < 0.1;
y(mask) = 0;
If you have any more questions, then attach your data and code to read it in with the paperclip icon after you read this:
  1 Comment
Abdelhakim Souidi
Abdelhakim Souidi on 1 Jul 2022
thank you so much your answer is very helpful, but some amplitude from the desired wave are distored how to overcome this issue

Sign in to comment.


Jeffrey Clark
Jeffrey Clark on 30 Jun 2022
@Abdelhakim Souidi , I couldn't determine if you had the data or just an image to work with. If you are working with just an image the someone with image processing could provide a more appropriate response. If you are looking to separate the signals in data you have (although they look very separate with the little information you provided) then from the information provided frequency is expected not to be a discriminator and only amplitude at the peaks are. If that is the case please refer to my original comment as one approach to your answer. Please mark as answered if you find this true.
  7 Comments
Abdelhakim Souidi
Abdelhakim Souidi on 1 Jul 2022
Edited: Abdelhakim Souidi on 1 Jul 2022
@Jeffrey Clark I am confused, can you repost all the code my friend this is it
Jeffrey Clark
Jeffrey Clark on 1 Jul 2022
@Abdelhakim Souidi , I assume you are getting the peaks you wanted with your post:
[pks,locs, wdth] = findpeaks(sh1,ty,'MinPeakProminence',0.1,'Annotate','extents','WidthReference','halfheight');
As I said I don't have access to findpeaks but reading the documentation the parameters 'Annotate' and after only apply to when there are no returns requested and a plot is made, so only 'MinPeakProminence',0.1 affects the default behavior. In addition, the data you posted in signal.mat and time.mat doesn't match the figures you posted (compare with my plot from these files in previous post). If the MinPeakProminence you used actually found the large peaks in your plots I don't think the value would work for the data you sent, and I'm not sure what it would need to be to find both the small and large peaks without picking up peaks nearby. Perhaps you should look into using both 'MinPeakHeight' and 'MinPeakDistance' instead.
I added documentation of what the code I provided expectations are:
% locs is the rerurn from findpeaks for the time (ty) relative data for the large amplitude peaks
% locs2 is the rerurn from findpeaks for the times for large and small peaks
% sh1 is the sampled signal amplitude
% ty is the time associated with each sample amplitude
discard = find(~any(locs'==locs2));
for i = discard % clear amplitudes associated with the small peaks
if i == 1
sh1(ty<(locs2(i+1)+locs2(i))/2) = 0;
elseif i==length(locs2)
sh1(ty>=(locs2(i-1)+locs2(i))/2) = 0;
else
sh1(ty>=(locs2(i-1)+locs2(i))/2 & ty<(locs2(i+1)+locs2(i))/2) = 0;
end
end

Sign in to comment.

Products


Release

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!