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.
10 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Sumanta Kumar Dutta
am 18 Sep. 2024
Kommentiert: Voss
am 19 Sep. 2024
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) -
2 Kommentare
Voss
am 18 Sep. 2024
@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.
Akzeptierte Antwort
Voss
am 18 Sep. 2024
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 Kommentare
Voss
am 19 Sep. 2024
You're welcome!
The cone of influence is the third output from cwt, as described here:
See also:
Weitere Antworten (1)
Jonas
am 18 Sep. 2024
Bearbeitet: Jonas
am 19 Sep. 2024
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 Kommentare
Jonas
am 19 Sep. 2024
you can get the cone of influence y coordinates from the third output of the cwt function
Siehe auch
Kategorien
Mehr zu Continuous Wavelet Transforms finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!