Plot scalogram of Morlet wavelet transform

Hi there,

I am trying to plot a scalogram of time-series data from an Excel file using code instead of the Wavelet Time-Frequency Analyzer Toolbox. However, I am not obtaining the same result. Here is the code:

% scalogram_single_24h_FINAL.m

% CWT Scalogram – matches Wavelet Time-Frequency Analyzer app output exactly

%

% Rendering pipeline confirmed from official MathWorks documentation:

% pcolor(t, f, abs(wt)) + shading flat + YScale log

% COI: plot(t, coi, 'w--')

% Freq ticks: powers of 2 between min/max f (app default)

% Signal: z-score normalised before CWT

%

% Settings (from your toolbox screenshot):

% Wavelet : Morlet ('amor')

% VoicesPerOct : 24

% FreqLimits : [3.1261e-05, 3.6828] Hz

% ExtendSignal : Symmetric (true)

% Fs : 10 Hz

% ---------------------------------------------------------------

clearvars; close all;

%% ── 1. PARAMETERS ──────────────────────────────────────────────

Fs = 10;

hours_total = 24;

samplesPerHour = 3600 * Fs;

requiredSamples = hours_total * samplesPerHour;

voicesPerOct = 24;

waveletName = 'amor'; % Morlet (your toolbox selection)

fLow = 3.1261e-05; % Hz

fHigh = 3.6828; % Hz

cmap = parula(256);

maxDisplayCols = 4000; % decimate for rendering speed

%% ── 2. ASK EXPORT PREFERENCE FIRST (avoids graphics timeout) ───

choice = questdlg('Export figure after plotting?', ...

'Export', 'Yes', 'No', 'No');

doExport = strcmp(choice, 'Yes');

outfn = 0; outpath = '';

if doExport

[outfn, outpath] = uiputfile( ...

{'*.png','PNG image'; '*.pdf','PDF vector'; '*.fig','MATLAB figure'}, ...

'Save scalogram as ...');

if isequal(outfn, 0)

doExport = false;

disp(' Export cancelled - figure will still be plotted.');

end

end

%% ── 3. LOAD EXCEL DATA ─────────────────────────────────────────

[fn, fpath] = uigetfile( ...

{'*.xlsx;*.xls','Excel Files (*.xlsx, *.xls)'; '*.*','All Files'}, ...

'Select Excel file - single data column');

if isequal(fn, 0), error('File selection cancelled.'); end

raw = readmatrix(fullfile(fpath, fn));

if isempty(raw), error('No numeric data found.'); end

x = raw(:,1);

x = x(:);

%% ── 4. LENGTH CHECK ────────────────────────────────────────────

N = numel(x);

if N < requiredSamples

warning('Data has %d samples (%.2f h); expected %d (24 h).', ...

N, N/samplesPerHour, requiredSamples);

else

x = x(1:requiredSamples);

N = requiredSamples;

end

t_hours = (0:N-1).' / Fs / 3600; % decimal hours

%% ── 5. Z-SCORE NORMALISE (matches toolbox magnitude range) ──────

x_norm = (x - mean(x,'omitnan')) / std(x,'omitnan');

fprintf('\n Normalised: mean=%.2e std=%.4f\n\n', mean(x_norm), std(x_norm));

%% ── 6. COMPUTE CWT ─────────────────────────────────────────────

fprintf(' Computing CWT for %d samples ... please wait.\n', N);

tic

[wt, f, coi] = cwt(x_norm, Fs, waveletName, ...

'VoicesPerOctave', voicesPerOct, ...

'FrequencyLimits', [fLow fHigh], ...

'ExtendSignal', true);

fprintf(' CWT done in %.1f s.\n\n', toc);

% NOTE: cwt returns f ordered HIGH to LOW

% pcolor needs f LOW to HIGH → flipud required

% coi is already a 1×N vector in Hz, plot directly

%% ── 7. SCALOGRAM (linear magnitude – no log10) ──────────────────

scalogram = abs(wt); % exact variable name from app script

%% ── 8. DECIMATE FOR DISPLAY ────────────────────────────────────

if N > maxDisplayCols

step = floor(N / maxDisplayCols);

idx = 1 : step : N;

else

idx = 1 : N;

end

t_disp = t_hours(idx);

sc_disp = scalogram(:, idx);

coi_disp = coi(idx);

%% ── 9. FLIP f AND scalogram FOR pcolor (low freq → bottom) ──────

% pcolor + YScale log requires f low→high for correct orientation

f_plot = flipud(f);

sc_plot = flipud(sc_disp);

%% ── 10. FIGURE ─────────────────────────────────────────────────

fprintf(' Rendering figure ...\n');

fig = figure( ...

'Name', 'Magnitude Scalogram', ...

'Units', 'normalized', ...

'Position', [0.03 0.05 0.94 0.84], ...

'Color', 'w');

ax = axes(fig);

%% ── 11. PCOLOR (app-generated script recipe, confirmed in docs) ─

% From wavelettimefrequencyanalyzer-app.html:

% pcolor(dataTimes, frequency, scalogram)

% shading flat

% set(gca, YScale="log")

hp = pcolor(ax, t_disp, f_plot, sc_plot);

set(hp, 'EdgeColor', 'none');

shading(ax, 'flat');

set(ax, 'YScale', 'log');

colormap(ax, cmap);

% NO clim/caxis – CDataMapping auto, exactly as app

xlim(ax, [0 t_hours(end)]);

ylim(ax, [fLow fHigh]);

%% ── 12. Y-AXIS TICKS (powers of 2 – confirmed in docs) ─────────

% From cwt.html and app docs:

% minf = frequency(end); % f is high→low, so end = minimum

% maxf = frequency(1);

% freq = 2.^(round(log2(minf)):round(log2(maxf)));

minf = f(end); % f high→low: end is minimum

maxf = f(1); % f high→low: start is maximum

yTicks = 2 .^ (round(log2(minf)) : round(log2(maxf)));

yTicks = yTicks(yTicks >= fLow & yTicks <= fHigh);

set(ax, 'YTick', yTicks, ...

'YTickLabelMode', 'auto', ... % let MATLAB format naturally

'TickDir', 'out', ...

'FontSize', 9, ...

'Layer', 'top');

ylabel(ax, 'Frequency (Hz)', 'FontSize', 10);

%% ── 13. X-AXIS: numeric hours, tick every 5 h ───────────────────

xlabel(ax, 'Time (hrs)', 'FontSize', 10);

xTicks = 0 : 5 : hours_total;

set(ax, 'XTick', xTicks, ...

'XTickLabel', arrayfun(@num2str, xTicks, 'UniformOutput', false), ...

'XTickLabelRotation', 0);

%% ── 14. TITLE ───────────────────────────────────────────────────

title(ax, 'Magnitude Scalogram', 'FontSize', 12, 'FontWeight', 'bold');

%% ── 15. COI (app recipe: plot directly in Hz on log axis) ──────

% From app docs:

% plot(dataTimes, coi, "w--", LineWidth=2)

% coi is already in Hz – no clamping, no remapping needed

hold(ax, 'on');

plot(ax, t_disp, coi_disp, '--w', 'LineWidth', 2);

hold(ax, 'off');

%% ── 16. COLORBAR ────────────────────────────────────────────────

cb = colorbar(ax, 'eastoutside');

cLims = cb.Limits;

nTicks = 6;

cbTicks = linspace(cLims(1), cLims(2), nTicks);

magRange = cLims(2) - cLims(1);

if magRange < 0.01, fmt = '%.4f';

elseif magRange < 0.1, fmt = '%.3f';

else, fmt = '%.2f';

end

cb.Ticks = cbTicks;

cb.TickLabels = arrayfun(@(v) sprintf(fmt,v), cbTicks, 'UniformOutput', false);

cb.Label.String = 'Magnitude';

cb.Label.FontSize = 9;

cb.Label.Rotation = 270;

cb.Label.VerticalAlignment = 'bottom';

cb.FontSize = 8.5;

fprintf(' Magnitude range: %.4e ... %.4e\n', cLims(1), cLims(2));

fprintf(' Figure ready.\n\n');

%% ── 17. EXPORT ───────────────────────────────────────────────────

if doExport

outfull = fullfile(outpath, outfn);

try

exportgraphics(fig, outfull, 'Resolution', 300);

catch

saveas(fig, outfull);

end

fprintf(' Figure saved -> %s\n', outfull);

end

Attached will find the data.

Thanks in advance for your cooperation.

Regards.

[Chuguang Pan]

@Angel Lozada. You can examine the source code of cwt to see how the cwt plot the scalogram when it is called without outputs.

help cwt

 cwt - Continuous 1-D wavelet transform
    This MATLAB function returns the continuous wavelet transform (CWT) of
    x.

    Syntax
      wt = cwt(x)
      wt = cwt(x,wname)

      [wt,f] = cwt(___,fs)
      [wt,period] = cwt(___,ts)
      [wt,f,coi] = cwt(___)
      [wt,period,coi] = cwt(___,ts)

      [___,coi,fb] = cwt(___)
      [___,fb,scalingcfs] = cwt(___)
      [___] = cwt(___,Name=Value)

      cwt(___)

    Input Arguments
      x - Input signal
        vector | timetable
      wname - Analytic wavelet
        "morse" (default) | "amor" | "bump"
      fs - Sampling frequency
        positive scalar
      ts - Sampling period
        scalar duration

    Name-Value Arguments
      ExtendSignal - Extend input signal symmetrically
        true or 1 (default) | false or 0
      FrequencyLimits - Frequency limits
        two-element scalar vector
      PeriodLimits - Period limits
        two-element duration array
      VoicesPerOctave - Number of voices per octave
        10 (default) | integer from 1 to 48
      TimeBandwidth - Time-bandwidth product of the Morse wavelet
        60 (default) |
        scalar greater than or equal to 3 and less than or equal to 120
      WaveletParameters - Symmetry and time-bandwidth product of the Morse wavelet
        [3,60] (default) | two-element vector of scalars
      FilterBank - CWT filter bank
        cwtfilterbank object

    Output Arguments
      wt - Continuous wavelet transform
        matrix
      f - Scale-to-frequency conversions
        vector
      period - Scale-to-period conversions
        array
      coi - Cone of influence
        array of real numbers | array of durations
      fb - CWT filter bank
        cwtfilterbank object
      scalingcfs - Scaling coefficients
        real- or complex-valued vector

    Examples
      openExample('wavelet/CWTdefaultOfSpeechSignalExample')
      openExample('wavelet/CWTofSpeechSignalUsingMorletWaveletExample')
      openExample('wavelet/CWTofKobeEarthquakeDataExample')
      openExample('wavelet/CWTofComplexExponentialsExample')
      openExample('wavelet/SinusoidAndWaveletCoefficientAmplitudesExample')
      openExample('wavelet/UsingCWTFilterBankOnMultipleTimeSeriesExample')
      openExample('wavelet/CUDACodeFromCWTExample')
      openExample('wavelet/ChangeDefaultFrequencyAxisLabelsExample')
      openExample('wavelet/ChangeScalogramColorationExample')
      openExample('wavelet/ChangingTheTimeBandwidthExample')
      openExample('wavelet/PlotCWTScalogramInSubplotExample')

    See also Wavelet Time-Frequency Analyzer, dlcwt, icwt, cwtfreqbounds,
      cwtfilterbank, cwtLayer

    Introduced in Wavelet Toolbox in R2016b
    Documentation for cwt
       doc cwt