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

7 次查看（过去 30 天）
Ricardo Duarte
Ricardo Duarte 2023-10-10
评论： Mathieu NOE 2023-10-12
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 个评论
显示 1更早的评论
Dyuman Joshi
Dyuman Joshi 2023-10-10
编辑：Dyuman Joshi 2023-10-10
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2023-10-10
Ran in:
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 个评论
显示 2更早的评论

请先登录，再进行评论。

更多回答（2 个）

Mathieu NOE
Mathieu NOE 2023-10-10
Ran in:
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 个评论
显示 1更早的评论
Mathieu NOE
Mathieu NOE 2023-10-10
Ran in:
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)

请先登录，再进行评论。

Florian Bidaud
Florian Bidaud 2023-10-10
编辑：Florian Bidaud 2023-10-10
Ran in:
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];

类别

在 Help Center 和 File Exchange 中查找有关 Just for fun 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by