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

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.

[Florian Bidaud]

Is your issue drawing the rectangle or finding the right coordinates of the variation ? (Or both ?)

[Ricardo Duarte]

Hello @Florian Bidaud,

Thank you for your fast answering. In fact is both.

[Dyuman Joshi]

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.

Answer 1

[Star Strider]

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)

.

[Ricardo Duarte]

Hello @Star Strider,

Thank you for your help.

When I use this part of the code

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

I get the following error

"Unable to perform assignment because brace indexing is not supported for variables of this type."

[Star Strider]

My pleasure!

That probably means that you have defined ‘x’ or ‘y’ earlier as something other than a cell array (otherwise it would over-write the cell array). One option is to change the variable names, for example:

for k = 1:numel(idx) % Get Contour (x,y) Coordinates

clen = h(2,idx(k));

xc{k} = h(1,idx(k)+1:idx(k)+clen);

yc{k} = h(2,idx(k)+1:idx(k)+clen);

end

[xmax,xmin] = bounds(cat(2,xc{:})) % Extreme Values

[ymax,ymin] = bounds(cat(2,yc{:})) % Extreme Values

Remember to change any other references to them later in the code, if you use them.

I could not plot the lines on the spectrogram output plot itself, so I used the matrix output from it and plotted the lines on the surf plot of it instead.

Also, use:

axis('tight')

to eliminate the white areas around the surf plot if there are any in your plot. (I did not use that here, simply because I forgot to.)

.

[Ricardo Duarte]

Perfect! Thank you @Star Strider, this solve the problem.

[Star Strider]

As always, my pleasure!

Answer 2

[Mathieu NOE]

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)

[Ricardo Duarte]

Thank you @Mathieu NOE,

I know how to draw the rectangle.

The ploblem is that I want to draw automatically a rectangle when the energy "jump" happens/starts in time and in frequency and then finish the rectangle when the energy comes to the "base level".

Thank you

[Mathieu NOE]

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)

[Mathieu NOE]

huh, so my suggestion didn't please you ? (just kidding)

Answer 3

[Florian Bidaud]

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