Clear Filters
Clear Filters

Spectrogram output versus figure

17 views (last 30 days)
DLR
DLR on 5 Sep 2024
Commented: Paul on 13 Sep 2024 at 16:42
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!

Answers (1)

Paul
Paul on 5 Sep 2024
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)
  2 Comments
DLR
DLR on 13 Sep 2024 at 1:45
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
Paul on 13 Sep 2024 at 16:42
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.

Sign in to comment.

Products


Release

R2024a

Community Treasure Hunt

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

Start Hunting!