Spectrogram output versus figure

When using spectrogram to calculate power spectral density for a large dataset, the plot produced by directly running the command with no output arguments is slightly different from what I get when I run it with an output argument and plot that output myself.

For example, the following plots spectrograms for a subset of my data and compares the direct result of spectrogram with the output raster:

load('exampledata.mat')
window = 10000;
noverlap = round(0.1*window);
fs = 10000;

figure(1);clf
t = tiledlayout(3,1, 'TileIndexing', 'columnmajor');

% Plot spectrogram directly from spectrogram function
s(1) = nexttile;
spectrogram(exampledata,window,noverlap,window,fs);
title('Original spectrogram plot (real data)')
PSDdB_orig = get(get(gca,'Children'),'CData');
caxis([min(PSDdB_orig,[],'all') max(PSDdB_orig,[],'all')])

% Calculate power spectral density and plot as decibels
[~,F,T,PSD] = spectrogram(exampledata,window,noverlap,window,fs);
PSDdB = (10.*log10(PSD))';

s(2) = nexttile;
imagesc(F./1000,T./60,PSDdB)
title('Spectrogram output raster (real data)')
c = colorbar;
caxis([min(PSDdB,[],'all') max(PSDdB,[],'all')])
set(get(c,'label'),'string','Power/frequency (dB/Hz)','Rotation',90);

% Calculate ratio of the two approaches above and plot that
PSDdiff = PSDdB_orig./PSDdB;

s(3) = nexttile;
imagesc(F./1000,T./60,PSDdiff)
title('Original plot data divided by output raster data')
c = colorbar;
caxis([min(PSDdiff,[],'all') max(PSDdiff,[],'all')])
set(get(c,'label'),'string','original/output','Rotation',90);
set(s(2:3),'YDir','normal')
xlabel(s(2:3), s(1).XLabel.String)
ylabel(s(2:3), s(1).YLabel.String)

From the figures above, it looks like they produce nearly the same result, except at high frequencies (the middle figure looks more different than the first mostly due to the different color scales).

Comparing the ranges and plotting a cross-section (example spectra at a single time) from each one shows that the difference seems to be that spectrogram's direct plotting method is cutting off local minima:

disp(['min/max of original spectrogram plot: ',num2str(min(PSDdB_orig,[],'all')),', ',num2str(max(PSDdB_orig,[],'all'))])

min/max of original spectrogram plot: -156.5351, -59.1828

disp(['min/max of spectrogram output raster: ',num2str(min(PSDdB,[],'all')),', ',num2str(max(PSDdB,[],'all'))])

min/max of spectrogram output raster: -196.3296, -59.1828

% Plot example spectra from each approach
figure(2);clf
subplot(2,1,1);
plot(F./1000,PSDdB(10,:),F./1000,PSDdB_orig(10,:),'--');
legend('Spectrogram output raster','Original spectrogram')
title('Example spectra')
xlabel('Frequency (kHz)')
ylabel('Power/frequency (dB/Hz)')

% zoom in on high frequency range where values differ
subplot(2,1,2,copyobj(gca,gcf))
xlim([4.75 4.85])
title('high frequency close-up')

But weirdly, I can't reproduce this with randomly generated synthetic data:

x = randn(length(exampledata),1);

figure(3);clf

t = tiledlayout(3,1, 'TileIndexing', 'columnmajor');

s(1) = nexttile;
spectrogram(x,window,noverlap,window,fs);
title('Original spectrogram plot (synthetic data)')
PSDdB_orig = get(get(gca,'Children'),'CData');
caxis([min(PSDdB_orig,[],'all') max(PSDdB_orig,[],'all')])

[~,F,T,PSD] = spectrogram(x,window,noverlap,window,fs);
PSDdB = (10.*log10(PSD))';
PSDdiff = PSDdB_orig./PSDdB;

s(2) = nexttile;
imagesc(F./1000,T./60,PSDdB)
title('Spectrogram output raster (synthetic data)')
c = colorbar;
caxis([min(PSDdB,[],'all') max(PSDdB,[],'all')])
set(get(c,'label'),'string','Power/frequency (dB/Hz)','Rotation',90);

s(3) = nexttile;
imagesc(F./1000,T./60,PSDdiff)
title('Original plot data divided by output raster data')
c = colorbar;
caxis([min(PSDdiff,[],'all') max(PSDdiff,[],'all')])
set(get(c,'label'),'string','original/output','Rotation',90);

set(s(2:3),'YDir','normal')
xlabel(s(2:3), s(1).XLabel.String)
ylabel(s(2:3), s(1).YLabel.String)

disp(['min/max of original spectrogram plot: ',num2str(min(PSDdB_orig,[],'all')),', ',num2str(max(PSDdB_orig,[],'all'))])

min/max of original spectrogram plot: -91.3499, -26.1202

disp(['min/max of spectrogram output raster: ',num2str(min(PSDdB,[],'all')),', ',num2str(max(PSDdB,[],'all'))])

min/max of spectrogram output raster: -91.3499, -26.1202

I do notice that the minimum PSD value for my data (-196 dB) is significantly lower than the minimum of the synthetic data. Is there a lower limit where spectrogram starts to censor spectra (say, around -156.5351 dB)? Or is something else going on here? Any help would be much appreciated!! Thanks!

Answer 1

[Paul]

Hi DLR,

Stepping through spectrogram shows that, when output argruments aren't specified, the PSD data that's plotted is offset by eps("double") when creating the plot. See modified code below.

load('exampledata.mat')
window = 10000;
noverlap = round(0.1*window);
fs = 10000;
figure(1);clf
t = tiledlayout(3,1, 'TileIndexing', 'columnmajor');
% Plot spectrogram directly from spectrogram function
s(1) = nexttile;
spectrogram(exampledata,window,noverlap,window,fs);
title('Original spectrogram plot (real data)')
PSDdB_orig = get(get(gca,'Children'),'CData');
caxis([min(PSDdB_orig,[],'all') max(PSDdB_orig,[],'all')])
% Calculate power spectral density and plot as decibels
[~,F,T,PSD] = spectrogram(exampledata,window,noverlap,window,fs);

%PSDdB = (10.*log10(PSD))';
PSDdB = 10*log10(PSD + eps("double"))'; % what's actually plotted when output arguments not specified.


s(2) = nexttile;
imagesc(F./1000,T./60,PSDdB)
title('Spectrogram output raster (real data)')
c = colorbar;
caxis([min(PSDdB,[],'all') max(PSDdB,[],'all')])
set(get(c,'label'),'string','Power/frequency (dB/Hz)','Rotation',90);
% Calculate ratio of the two approaches above and plot that
PSDdiff = PSDdB_orig./PSDdB;
s(3) = nexttile;
imagesc(F./1000,T./60,PSDdiff)
title('Original plot data divided by output raster data')
c = colorbar;
% caxis([min(PSDdiff,[],'all') max(PSDdiff,[],'all')])
set(get(c,'label'),'string','original/output','Rotation',90);
set(s(2:3),'YDir','normal')
xlabel(s(2:3), s(1).XLabel.String)
ylabel(s(2:3), s(1).YLabel.String)

[DLR]

Hi Paul, thanks so much for the answer! Just to be sure that I'm clear on what's happening here after reading the eps documentation... it sounds like this is just intended to avoid any numbers below the double precision limit, correct? Do you know whether this implies that working with values below 10^-16 should generally be avoided? I was under the impression that double precision limited the number of digits, not the absolute range of values--my data generally involves very small values (microstrain), but not a very large range, so I didn't think I needed to worry about precision. Am I mistaken? Thanks for any insights!

[Paul]

eps is not the double precision limit, at least that's not how I would state it.

As best I can tell, the purpose of adding eps to PSD for purposes of plotting is just to limit the range of the colors on the plot. On the other hand, there may have been better ways to do that besides adding a small value to every single element of PSD, so maybe there really is something more to it. But I doubt it.

"Do you know whether this implies that working with values below 10^-16 should generally be avoided?"

No I don't. Double precision floating point can represent much small numbers:

realmin("double")

ans = 2.2251e-308

Whether or not it's sensible to use such small numbers is beyond my knowledge.

Maybe you can change the units on your data or use the outputs from spectrogram and roll your own plots if needed.