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.

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]

@Sumanta Kumar Dutta: Looks like your images accidentally got deleted when I edited the question to apply code formatting. Sorry about that; add them again if necessary.

[Sumanta Kumar Dutta]

@Voss No problem. The kind people on this forum have already suggested some solutions for my problem.

Answer 1

[Voss]

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)

[Sumanta Kumar Dutta]

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?

[Voss]

You're welcome!

The cone of influence is the third output from cwt, as described here:

https://www.mathworks.com/help/wavelet/ref/cwt.html#bvaytl3-coi

You can plot it on top of the surface using area, fill, or patch.

See also:

https://www.mathworks.com/help/wavelet/ug/boundary-effects-and-the-cone-of-influence.html

Answer 2

[Jonas]

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')

[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?

[Jonas]

you can get the cone of influence y coordinates from the third output of the cwt function