How to draw a rectangle around the area that shows an energy variation

13 views (last 30 days)
Hello all,
I have several wav files and I want to draw a rectangle around the area here a specific energy variation happend. Let's say, since the beggining til the end (in seconds)
In fact is something similar than the image attached
Thanks in advance for your help.
  3 Comments
Dyuman Joshi
Dyuman Joshi on 10 Oct 2023
Edited: Dyuman Joshi on 10 Oct 2023
Please attach the code and the data you are working with. Use the paperclip button to attach.
Also, specify the region around which you need to draw a rectangle.

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 10 Oct 2023
It would help to have your data.
An adaptive approach using contour3
Fs = 500;
L = 300;
t = linspace(0, L*Fs, L*Fs+1)/Fs;
signal = sum(sin((100:50:200).'*2*pi*t) .* exp(-(t-150).^2*1E-3));
figure
plot(t, signal)
grid
xlabel('Time (s)')
ylabel('Amplitude')
[s,f,t,p] = spectrogram(signal,[],[],[],Fs); % Return Extra Outputs
figure
surf(f,t,p.')
colormap(turbo)
[h,c] = contour3(f,t,p.', 20); % Get Contours
Lvls = c.LevelList
Lvls = 1×20
0.5617 1.1235 1.6852 2.2470 2.8087 3.3705 3.9322 4.4940 5.0557 5.6175 6.1792 6.7409 7.3027 7.8644 8.4262 8.9879 9.5497 10.1114 10.6732 11.2349
ChooseLevel = 1;
idx = find(h(1,:) == Lvls(ChooseLevel));
for k = 1:numel(idx) % Get Contour (x,y) Coordinates
clen = h(2,idx(k));
x{k} = h(1,idx(k)+1:idx(k)+clen);
y{k} = h(2,idx(k)+1:idx(k)+clen);
end
[xmax,xmin] = bounds(cat(2,x{:})) % Extreme Values
xmax = 99.9763
xmin = 200.0246
[ymax,ymin] = bounds(cat(2,y{:})) % Extreme Values
ymax = 100.3145
ymin = 199.6891
colormap(turbo)
title('Contours')
figure
% spectrogram(signal,[],[],[],Fs)
surf(f,t,mag2db(p.'), 'EdgeColor','interp')
hold on
plot3([1 1]*xmin, [ymin ymax], [1 1]*Lvls(ChooseLevel), '-k', 'LineWidth',2)
plot3([xmin xmax], [1 1]*ymin, [1 1]*Lvls(ChooseLevel), '-k', 'LineWidth',2)
plot3([xmin xmax], [1 1]*ymax, [1 1]*Lvls(ChooseLevel), '-k', 'LineWidth',2)
plot3([1 1]*xmax, [ymin ymax], [1 1]*Lvls(ChooseLevel), '-k', 'LineWidth',2)
hold off
colormap(turbo)
view(0,90)
.
  4 Comments

Sign in to comment.

More Answers (2)

Mathieu NOE
Mathieu NOE on 10 Oct 2023
hello
try this
IMO, the rectangle position is to be computed from your spectrogram results (define thresholds)
Fs = 1e3;
t=0:1/Fs:5; % 5 secs @ 1kHz sample rate
y=chirp(t,10,3,100,'q'); % Start @ 10Hz, cross 100Hz at t=3sec
% amplitude is increasing with time
y = y.*t./t(end);
[S,F,T] = spectrogram(y,hanning(128),64,128,Fs);
imagesc(T,F,20*log10(abs(S)));
caxis([-50 30]);
set(gca,'Ydir','normal');
colormap('jet');
colorbar('vert');
title('Quadratic Chirp: start at 100Hz and cross 200Hz at t=1sec');
xlabel('Time (s)')
ylabel('Frequency (Hz)')
% rectangle
fmin = 100;
fmax = 200;
tmin = T(end)/2;
tmax = T(end)*0.99;
rectangle('Position',[tmin fmin (tmax-tmin) (fmax-fmin)],'Curvature',0.25*[1 1],'Linewidth',3)
  3 Comments
Mathieu NOE
Mathieu NOE on 10 Oct 2023
hello Ricardo
yes of course I suspected drawing the rectangle is not the issue here :)
here my new suggestion : assuming we can find the peak value of S , you can define a threshold (here 30 dB below the peak value (in dB)) so the rectangle is automatically drawn
you can change this threshold value in this line :
threshold = S_dB_max - 30; % define here threshold you want to draw the rectangle : here 30 dB below the peak amplitude (= S_dB_max)
full code
Fs = 1e3;
t=0:1/Fs:5; % 5 secs @ 1kHz sample rate
y=chirp(t,10,3,100,'q'); % Start @ 10Hz, cross 100Hz at t=3sec
% amplitude is changing with time
window = exp(-(t-3.5).^2/0.75); % gaussian window
y = y.*window;
[S,F,T] = spectrogram(y,hanning(128),64,128,Fs);
S_dB = 20*log10(abs(S));
S_dB_max = max(S_dB,[],'all');
threshold = S_dB_max - 30; % define here threshold you want to draw the rectangle : here 30 dB below the peak amplitude (= S_dB_max)
ind = S_dB>=threshold;
imagesc(T,F,S_dB);
caxis([-60 30]);
set(gca,'Ydir','normal');
colormap('jet');
colorbar('vert');
title('Quadratic Chirp');
xlabel('Time (s)')
ylabel('Frequency (Hz)')
% rectangle
indF = any(ind,2); % indices for F vector
indT = any(ind,1); % indices for T vector
fmin = min(F(indF),[],'all');
fmax = max(F(indF),[],'all');
tmin = min(T(indT),[],'all');
tmax = max(T(indT),[],'all');
rectangle('Position',[tmin fmin (tmax-tmin) (fmax-fmin)],'Curvature',0.25*[1 1],'Linewidth',3)

Sign in to comment.


Florian Bidaud
Florian Bidaud on 10 Oct 2023
Edited: Florian Bidaud on 10 Oct 2023
Since you said your issue was both finding the variation and drawing the rectangle.
To find the variations, I would suggest using diff, and set a trigger value from which you decide this variation is interesting (that could be with a percentage from the others for example). This part, you need to think a bit about what variation is interesting.
Let's take the following example :
z = [1 3 5 6 11
1 3 6 5 8
1 2 3 9 7
1 1 2 3 6
1 1 1 2 4
];
zp = diff(z)
zp = 4×5
0 0 1 -1 -3 0 -1 -3 4 -1 0 -1 -1 -6 -1 0 0 -1 -1 -2
Here, if you set a limit of variation to 2, you can, for example find the x and y coordinate of your limits:
var_lim = 2;
x = find(sum(diff(z)<-var_lim),1)
x = 3
y = max(size(z))-find(sum(diff(z')>var_lim),1)
y = 4
Then from these, you can draw the rectangle as follows :
(in my case I chose 1 and 5 as I deal with indexes from 1 to 5, but this should be replace by your own limits on the plot)
f = figure;
hold on
s = surf(z);
p = plot([x 5 5 x x],[1 1 y y 1]);
p.ZData = [50 50 50 50 50]; %to put it on top
p.Color = 'black';
p.LineWidth = 3;
f.Children.View = [0 90];

Products

Community Treasure Hunt

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

Start Hunting!