How to plot scalogram in time period?

Question

Jan Ali 2021-4-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/812585-how-to-plot-scalogram-in-time-period

编辑： Jan Ali 2021-5-3

Hi eveyone!

I want to analyze an event in a signal in a certain time span, but using the cwt(x, fs) built in function plots in Scalogram at time axes from 0 as seen in the below plot.

How can I plot the scalogram in a time interval? I appreciate if anyone can help me with this.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Wayne King 2021-4-29

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/812585-how-to-plot-scalogram-in-time-period#answer_687820

在 MATLAB Online 中打开

Hi Jan, you do not mention the release you are using, but you can put the data in a timetable with the RowTimes equal to whatever they happened to be and those should be preserved in the convenience plot.

For example:

Fs = 1e3;
t = 0:1/Fs:1;
x = cos(2*pi*32*t).*(t>=0.1 & t<0.3)+sin(2*pi*64*t).*(t>0.7);
wgnNoise = 0.05*randn(size(t));
x = x+wgnNoise;
t = t+4;
% Here make the timetable
tt = timetable(x(:),'RowTimes',seconds(t'));
cwt(tt)

Now you see that the x-axis in the plot has the correct range. Here t was a duration array, but it can also be a datetime array.

Hope that helps,

Wayne

4 个评论
显示 2更早的评论隐藏 2更早的评论

Wayne King 2021-4-30

在 MATLAB Online 中打开

Hi Jan, that is because the convenience plot tries to scale the time axis in a way that is amenable to display. Since your orignal time axis has seconds in 5,000 range, it tries to be smart and thinks that hours is better.

However, you can have full control over the plot by returning the outputs of cwt() and plotting those yourself. Then a time table is not required.

For example:

Fs = 1e3;
t = 0:1/Fs:1;
x = cos(2*pi*32*t).*(t>=0.1 & t<0.3)+sin(2*pi*64*t).*(t>0.7);
wgnNoise = 0.05*randn(size(t));
x = x+wgnNoise;
t = t+5083;
[wt,f,coi] = cwt(x,Fs);
ax = newplot;
imagesc(ax,t,f,abs(wt));
ax.YScale = 'log';
ax.YDir = 'normal';
axis tight
hold on
plot(ax,t,coi,'w--')
xlabel('Seconds')
ylabel('Hz')
%% 
% Using cwtfilterbank is actually more efficient
fb = cwtfilterbank('SignalLength',numel(x),'SamplingFrequency',1e3);
[wt,f,coi] = fb.wt(x);
ax = newplot;
imagesc(ax,t,f,abs(wt));
ax.YScale = 'log';
ax.YDir = 'normal';
axis tight
hold on
plot(ax,t,coi,'w--')
xlabel('Seconds')
ylabel('Hz')

I've shown above that actually using cwtfilterbank and then the wt() method is more efficient than using cwt() because cwt() uses cwtfilterbank. However, if you like the functional API, you can use that.

Note I've left the shading of the region outside the "cone of influence" here but that is easy to add if you also want that. But that is just a visualize convenience. Also, you do not need to plot the cone of influence if you don't care about that. In that case, just don't use the plot(ax,t,coi,'w--') code.

Hope that helps,

Wayne

Jan Ali 2021-4-30

Thanks so much Wayne! I really appreciate your kind help.

Wish you all the bests!

Jan Ali 2021-5-3

编辑：Jan Ali 2021-5-3

Hi Wayne,

Now the y axis scales are impaired. It does not show the actual value. For exp, the 60 Hz frequency is now measured under 10 Hz. Eeven if I change the ax.YScale = 'linear'; still no actual value is plotted. Do you have any suggestion?

请先登录，再进行评论。

Answer 2

Wayne King 2021-5-3

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/812585-how-to-plot-scalogram-in-time-period#answer_690355

在 MATLAB Online 中打开

Hi Jan, can you try surf() instead of imagesc()? Alternatively, can you attach a sample time series along with the necessary information like sample rate? Note in the following example, if you put a datatip on the two localized sinusoids, the frequencies are correct.

Fs = 1e3;
t = 0:1/Fs:1;
x = cos(2*pi*32*t).*(t>=0.1 & t<0.3)+sin(2*pi*64*t).*(t>0.7);
wgnNoise = 0.05*randn(size(t));
x = x+wgnNoise;
[cfs,f] = cwt(x,1e3);
ax = newplot;
surf(ax,t,f,abs(cfs));
view(0,90); shading interp;
ax.YScale = 'log';
ax.YDir = 'normal';
axis tight
hold on
plot(ax,t,coi,'w--')
xlabel('Seconds')
ylabel('Hz')

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Jan Ali 2021-5-3

Thanks alot! Much appreciate it.

Sorry for distubing you again. I am wondering how to plot that Magnitude color bar on the right?

Additionally, do you have any suggestion on the same issue for the spectrogram using stft?

I also need the same measures for the time axis.

请先登录，再进行评论。

Answer 3

Wayne King 2021-5-3

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/812585-how-to-plot-scalogram-in-time-period#answer_690975

在 MATLAB Online 中打开

Hi Jan, you can do the following:

Fs = 1e3;
t = 0:1/Fs:1;
x = cos(2*pi*32*t).*(t>=0.1 & t<0.3)+sin(2*pi*64*t).*(t>0.7);
wgnNoise = 0.05*randn(size(t));
x = x+wgnNoise;
[cfs,f] = cwt(x,1e3);
ax = newplot;
surf(ax,t,f,abs(cfs));
ax.YScale = 'log';
view(0,90); shading interp;
axis tight;
cb = colorbar(ax);
cb.Title.String = 'magnitude';

As far as the spectrogram, that does not have the same time resolution as the original data. So you'll need the times from the stft() function and then shift those by t(1)

Fs = 1e3;
t = 0:1/Fs:10;
t = t+5000; % say it starts at t=5e3 seconds
x = cos(2*pi*32*t).*(t>=5001 & t<5003)+sin(2*pi*64*t).*(t>5005);
wgnNoise = 0.05*randn(size(t));
x = x+wgnNoise;
% Get the window center times from STFT
[S,F,T] = stft(x,1e3,'FrequencyRange','on');  %with whatever other options you want
% Shift by t(1)
T = T+t(1);
ax = newplot;
surf(T,F,20*log10(abs(S))); shading interp;
view(0,90)
axis tight;
cb = colorbar(ax);
cb.Title.String = 'Power(db)';
xlabel('Seconds')
ylabel('Hz')

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Jan Ali 2021-5-3

编辑：Jan Ali 2021-5-3

Many thanks Wayne!

Now the coi region is disappeared. Have a clue?

I really appreciate your kind support and the time you spent.

Best regards,

Jan

请先登录，再进行评论。

How to plot scalogram in time period?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

4 个评论
显示 2更早的评论隐藏 2更早的评论

更多回答（2 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

How to plot scalogram in time period?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

4 个评论 显示 2更早的评论隐藏 2更早的评论

更多回答（2 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论