I am getting wrong results for the Continuous wavelet transform (CWT) on converting the frequency axis (Y Axis) from log scale to linear scale for the same frequency limits.

11 次查看(过去 30 天)
Hello All,
I have been trying to apply CWT on a signal. The function cwt(signal, Fs) gives the correct frequency when the Frequency scale is in the log scale but the incorrect frequency on converting the axis to linear scale.
Please suggest how to correct this problem. Will upgrading to newer version of MATLAB help?
Code -
fs = 20000; % Sampling frequency in Hz
t1 = 0:1/fs:10-1/fs; % Time vector for the first 10 seconds
t2 = 10:1/fs:20-1/fs; % Time vector for 10 to 20 seconds
t3 = 20:1/fs:30-1/fs; % Time vector for 20 to 30 seconds
% Create Signal
% 1. Sine wave of amplitude 10 and frequency 1 kHz for the first 10 seconds
signal1 = 10 * sin(2 * pi * 1000 * t1) + 3 * sin(2 * pi * 3000 * t1);
% 2. Combination of 2 sine waves:
signal2 = 4 * sin(2 * pi * 2000 * t2) + 6 * sin(2 * pi * 500 * t2);
% 3. Sine wave of amplitude 20 and frequency 300 Hz for time 20-30 seconds
signal3 = 8 * sin(2 * pi * 300 * t3) + 5 * sin(2 * pi * 1500 * t3);
% Concatenate all parts of the signal to form the complete signal
completeSignal = [signal1, signal2, signal3];
% Time vector for the complete signal
t = [t1, t2, t3];
% Wavelet Analysis
frequencyLimits = [0 10000]; % Frequency range from 0 Hz to 10,000 Hz
voicesPerOctave = 48; % Choosing 48 voices per octave for finer frequency resolution
% Perform Continuous Wavelet Transform
% [wt, f] = cwt(completeSignal, 'amor', fs, 'FrequencyLimits', frequencyLimits, 'VoicesPerOctave', voicesPerOctave);
cwt(completeSignal, 'amor', fs);
% Plot the results
% figure (2);
% imagesc(t, f, abs(wt));
% axis xy;
% xlabel('Time (s)');
% ylabel('Frequency (Hz)');
% title(['Continuous Wavelet Transform (CWT) of the Signal (0 to 10 kHz) with ', num2str(voicesPerOctave), ' Voices/Octave']);
% colorbar;
Plot on using log scale for the Frequency Axis (Y Axis) -
Plot on using linear scale for the Frequency Axis (Y Axis) -

采纳的回答

Voss
Voss 2024-9-18
Use a surface rather than an image. A surface allows you to use arbitrary x and y values, whereas an image requires linearly spaced x and y, which is not the case with your f vector since it is logarithmically spaced.
See below for how to (more-or-less) reproduce the plot created by cwt() on both a linear and log y-scale. (The code wouldn't run here in the forum environment when specifying 48 voicesPerOctave or when plotting all columns of t and wt, so I made modifications to get it to run.)
fs = 20000; % Sampling frequency in Hz
t1 = 0:1/fs:10-1/fs; % Time vector for the first 10 seconds
t2 = 10:1/fs:20-1/fs; % Time vector for 10 to 20 seconds
t3 = 20:1/fs:30-1/fs; % Time vector for 20 to 30 seconds
% Create Signal
% 1. Sine wave of amplitude 10 and frequency 1 kHz for the first 10 seconds
signal1 = 10 * sin(2 * pi * 1000 * t1) + 3 * sin(2 * pi * 3000 * t1);
% 2. Combination of 2 sine waves:
signal2 = 4 * sin(2 * pi * 2000 * t2) + 6 * sin(2 * pi * 500 * t2);
% 3. Sine wave of amplitude 20 and frequency 300 Hz for time 20-30 seconds
signal3 = 8 * sin(2 * pi * 300 * t3) + 5 * sin(2 * pi * 1500 * t3);
% Concatenate all parts of the signal to form the complete signal
completeSignal = [signal1, signal2, signal3];
% Time vector for the complete signal
t = [t1, t2, t3];
% Wavelet Analysis
frequencyLimits = [0 10000]; % Frequency range from 0 Hz to 10,000 Hz
voicesPerOctave = 48; % Choosing 48 voices per octave for finer frequency resolution
% Perform Continuous Wavelet Transform
[wt, f] = cwt(completeSignal, 'amor', fs, 'FrequencyLimits', frequencyLimits);%, 'VoicesPerOctave', voicesPerOctave);
cwt(completeSignal, 'amor', fs);
yl = ylim();
% Plot the results
figure();
ax = gca();
% surface(ax,t, f/1000, abs(wt), 'EdgeColor', 'none');
surface(ax,t(1:10:end), f/1000, abs(wt(:,1:10:end)), 'EdgeColor', 'none');
ax.YLim = yl;
ax.YScale = 'log';
ax.YTick = 10.^(-4:0);
xlabel(ax,'Time (s)');
ylabel(ax,'Frequency (kHz)');
title(ax,'reproduction of cwt() plot with log y-scale');
colorbar(ax)
% Plot the results
figure()
ax = gca();
% surface(ax,t, f/1000, abs(wt), 'EdgeColor', 'none');
surface(ax,t(1:10:end), f/1000, abs(wt(:,1:10:end)), 'EdgeColor', 'none');
ax.YLim = yl;
xlabel(ax,'Time (s)');
ylabel(ax,'Frequency (kHz)');
title(ax,'reproduction of cwt() plot with linear y-scale');
colorbar(ax)
  2 个评论
Sumanta Kumar Dutta
编辑:Sumanta Kumar Dutta 2024-9-19
Thank you for your answer. The plots using the surface has a better resolution than what I obtained by applying interpolation and using image.
Do you know how to generate the cone of influence when we use imagesc or surface plot? Is it something that only appears when using the cwt() function to plot directly?

请先登录,再进行评论。

更多回答(1 个)

Jonas
Jonas 2024-9-18
编辑:Jonas 2024-9-19
you could interpolate the image pixels, but this makes the upper frequency values barely visible
fs = 10000; % Sampling frequency in Hz
t1 = 0:1/fs:3-1/fs; % Time vector for the first 10 seconds
t2 = 3:1/fs:6-1/fs; % Time vector for 10 to 20 seconds
t3 = 6:1/fs:9-1/fs; % Time vector for 20 to 30 seconds
% Create Signal
% 1. Sine wave of amplitude 10 and frequency 1 kHz for the first 10 seconds
signal1 = 10 * sin(2 * pi * 1000 * t1) + 3 * sin(2 * pi * 3000 * t1);
% 2. Combination of 2 sine waves:
signal2 = 4 * sin(2 * pi * 2000 * t2) + 6 * sin(2 * pi * 500 * t2);
% 3. Sine wave of amplitude 20 and frequency 300 Hz for time 20-30 seconds
signal3 = 8 * sin(2 * pi * 300 * t3) + 5 * sin(2 * pi * 1500 * t3);
% Concatenate all parts of the signal to form the complete signal
completeSignal = [signal1, signal2, signal3];
% Time vector for the complete signal
t = [t1, t2, t3];
% Wavelet Analysis
frequencyLimits = [0 fs/2];
voicesPerOctave = 48; % Choosing 48 voices per octave for finer frequency resolution
% Perform Continuous Wavelet Transform
[wt, f] = cwt(completeSignal, 'amor', fs, 'FrequencyLimits', frequencyLimits, 'VoicesPerOctave', voicesPerOctave);
wt=abs(wt);
% Plot the results
figure;
cwt(completeSignal, 'amor', fs, 'FrequencyLimits', frequencyLimits, 'VoicesPerOctave', voicesPerOctave);
f_linear = min(f):10:max(f); % to reduce memory a bit we use step of 10 here
% Interpolate the wavelet transform to the linear frequency scale
wt_linear = interp1(f, abs(wt), f_linear, 'makima');
t = (0:length(completeSignal)-1)/fs;
% Plot the scalogram with linear frequency scale
imagesc(t, f_linear, wt_linear);
set(gca,'YDir','normal')
  2 个评论
Sumanta Kumar Dutta
Thank you for replying. I too tried to interpolate but faced the same issues.
Do you know how to generate the cone of influence when we use imagesc or surface plot? Is it something that only appears when using the cwt() function to plot directly?

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Continuous Wavelet Transforms 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by