Functions specgram and spectrogram - what does the spectogram' 3rd dimension mean and what unit it has?

10 visualizaciones (últimos 30 días)
Peter
Peter el 29 de Jul. de 2023
Comentada: Star Strider el 3 de Ag. de 2023
I use two Matlab functions to plot a spectrogram of signal X: specgram and spectrogram - see the fragments of my code listed below. Both functions plot the same result. To show the legend of spectogram' 3rd dimension I use colorbar function. It is clear that the data of 3rd dimension are in format of decibels (dB).
What signifies the data of 3rd dimension: amplitude of signal X, or power of signal X, or power spectral density (P.S.D.) of signal X, or something else? And what is the unit (or units) of the data of spectogram' 3rd dimension? Thank You in advance.
fs = 10000;
WINDOW = 1024;
OVERLAP = WINDOW * 0.5;
...
figure('Position', scr_size);
specgram(X, WINDOW, fs);
colorbar;
...
[ss, ff, tt] = spectrogram(X, WINDOW, OVERLAP, WINDOW, fs);
figure('Position', scr_size);
contourf(tt, ff, (20 * log10(abs(ss))));
colorbar;

Iniciar sesión para responder a esta pregunta.

Respuesta aceptada

Star Strider
Star Strider el 29 de Jul. de 2023
The spectrogram output is the power spectral density, so in units of power/frequency, or dB/Hz.
It differs in that respect from pspectrum (with the 'spectrogram' option) that plots in terms of power (dB).
  6 comentarios
Mostrar 4 comentarios más antiguos
Peter
Peter el 3 de Ag. de 2023
Thanks a lot for your answer!
I tried this code:
t = 0:0.001:2;
hz1 = 10;
hz2 = 50;
s = (1 * cos(2*pi*t*hz1)) + (100 * cos(2*pi*t*hz2));
plot(t, s);
B = specgram(s, 1024, (1/0.001));
figure
surf(abs(B));
colorbar;
figure
surf(mag2db(abs(B)));
colorbar;
This is most likely true:

Iniciar sesión para comentar.

Más respuestas (1)

Paul
Paul el 1 de Ag. de 2023
Ran in:
Hi Peter,
Until I read this question, I had never really thought about the outputs of spectrogram in an absolute sense. I only look at the relative magnitudes across frequency and time to identify features of the underlying signal. But this question piqued my curiosity.
In summary, referring to this command, and assuming the fs input to spectrogram is specified:
% [s,f,t,ps] = spectrogram(....)
The first output, s, is the normalized DFT of the input, so its units would be the same units of the input (call that U). Note that the magnitude of s is the same regardless of the setting for the freqrange input.
The fourth output, ps, has units of power/Hz, or U^2/Hz, if the spectrumtype is 'psd' (the default), and units of U^2 if the spectrum type is 'power.'
When ps is the psd, ps relates to the power estimate of the time domain samples in the expected way. When ps is power, it's just scaled from the psd by a constant, and I can't figure out how it relates to the time domain samples (except for a rectangular window)
When spectrogram is called without output arguments, it makes a nice auto-plot, but actually plots db(ps,'power'). When ps is the psd in units of U^2/Hz, it seems to me that the auto-plot should have units of dB(power/Hz), not dB(power)/Hz as shown.
Here are examples to illustrate the points in the summary.
Define a simple signal
rng(100);
fs = 100;
N = 2048;
t = (0:N-1)/fs;
x = 5*sin(2*pi*10*t) + randn(1,N);
x = x - mean(x);
The estimated power in x[n] based on this realization of a subset of its samples is
estimatedpower = sum(abs(x).^2)/N
estimatedpower = 13.6252
in units of U^2, where U is the unit of x.
Compute its DFT normalized by its length and plot its amplitude.
X = fft(x)/N;
figure;
plot((0:N-1)/N*fs,abs(X));grid
Call spectrogram using a rectangular window, no overlap or padding. So this is basically a fancy way to compute the DFT
win = ones(1,N);
[s,f,t,psd] = spectrogram(x,win,0,N,fs);
[~,~,~,pow] = spectrogram(x,win,0,N,fs,[],'power');
figure
plot(f,abs(s)/sum(win));grid
We see that the default freqrange is one-sided, and that s, when normalized by the window, has the same amplitude as the normalized DFT. So, even though the freqrange is one-sided, spectrogram is not doing the mulitply-by-2 on the DFT.
For the rectangular window, the power estimates derived from psd and pow match the time domain estimate
sum(psd)*f(2) % f(2) is dF
ans = 13.6252
sum(pow)
ans = 13.6252
from which it's clear psd is in U^2/Hz and pow is in U^2, and that neither have been converted to dB (more discussion on this below).
Now do the same for a non-rectangular window
win = hamming(N,'periodic').';
[s,f,t,psd] = spectrogram(x,win,0,N,fs);
[~,~,~,pow] = spectrogram(x,win,0,N,fs,[],'power');
figure
plot(f,abs(s)/sum(win));grid
Unsurprisingly, the window does have an effect on the DFT, s, but the normalization is still based on the window.
The power estimate from the psd is reasonable
sum(psd)*f(2)
ans = 13.7374
However, the power estimate from pow is
sum(pow)
ans = 18.7217
The normalization factors used to compute pow and psd are based on the window. Here's the calculation for psd compared to the output of spectrogram
Xw = fft(x.*win);
cpsd = abs(Xw).^2./sum(win.^2)/fs;
cpsd = cpsd(1:numel(f)); % for onesided
% matlab doesn't multiply DC or Nyquist by 2 for onesided
cpsd = abs([cpsd(1) 2*cpsd(2:end-1) cpsd(end)]);
plot(f,psd,f,cpsd(1:numel(f))),grid,legend('spectrogram','verification')
xlim([9 11])
sum(cpsd(1:numel(f)))*f(2)
ans = 13.7374
The pow calculation is
cpow = abs(Xw).^2./sum(win)^2;
cpow = cpow(1:numel(f)); % for onesided
% matlab doesn't multiply DC or Nyquist by 2 for onesided
cpow = abs([cpow(1) 2*cpow(2:end-1) cpow(end)]);
plot(f,pow,f,cpow(1:numel(f))),grid,legend('spectrogram','verification')
xlim([9 11])
sum(cpow)
ans = 18.7217
We see that pow and psd are the same to within a scale factor across all frequencies.
sum(win)^2/(sum(win.^2)*fs)
ans = 15.0276
[max(psd./pow) min(psd./pow)]
ans = 1×2
15.0276 15.0276
I'm not exactly sure what pow actually represents in this case; maybe someone else can explain. I did find a reference for the pow calculation (reference), but I wasn't able to determine how exactly pow is related to the actual 'power' in the signal. However, the calculation of psd, which is is also in that same paper, seems to be easy to find in the literature (well, at least two other sources that I had handy).
Finally, let's look at the auto-plot for the psd, using the default parameters for spectrogram
figure
spectrogram(x,[],[],[],fs);
clim = get(gca,'CLim');
Now get the outputs
[~,f,t,psd] = spectrogram(x,[],[],[],fs);
The auto-plot can be recreated by converting psd to dB and using the same color limits
figure
pcolor(f,t,db(max(psd.',eps),'power'));shading flat,colorbar
set(gca,'CLim',clim)
Those plots look nearly (if not exactly) identical. But if psd is already in units of power/Hz, I don't understand how db(psd,'power') = 10*log10(psd) can be in units of dB/Hz as shown on the auto-plot.
As an aside, because the signal is stationary (or at least assumed to be so), the power estimated from the psd estimate of each window is roughly the estimated power in the signal
sum(psd).*f(2)
ans = 1×7
12.8010 13.9808 13.6328 14.6990 12.8732 13.3256 13.3747

Categorías

Más información sobre Time-Frequency Analysis en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2013b

Community Treasure Hunt

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

Start Hunting!

Translated by