Spectrogram Computation with Signal Processing Toolbox
Signal Processing Toolbox™ provides three functions that compute the spectrogram of a nonstationary signal. Each of the functions has different input arguments, default values, and outputs. The best choice for you depends on your particular application.
A nonstationary signal is a signal whose frequency content changes over time. The shorttime Fourier transform (STFT) is used to analyze how this frequency content changes as the signal evolves. The magnitude squared of the STFT is known as the spectrogram timefrequency representation of the signal.
The spectrogram is only one of several possible timefrequency representations. For an
overview of other timefrequency representations available in Signal Processing Toolbox and Wavelet Toolbox™, see TimeFrequency Gallery. For a treatment of
stationary signals using the periodogram
function, see Power Spectral Density Estimates Using FFT.
Functions for Spectrogram Computation
Signal Processing Toolbox has these functions that can be used to compute the spectrogram:
spectrogram
— Designed for maximum flexibility. Computes STFT and segmentbysegment power spectral densities or power spectra. Supports reassignment.stft
— Designed for invertibility and maximum control. Computes STFT. Used bydlstft
andstftLayer
. Supports multichannel input. Its counterpart, theistft
function, computes the inverse STFT.pspectrum
— Designed for ease of use. Computes power spectra. Used in the analysis scripts generated by Signal Analyzer. Can compute spectra of stationary signals and persistence spectra. Supports reassignment.
The spectrogram
function is used as reference in the
discussion that follows.
Category  Parameters  Function  

spectrogram  stft  pspectrum  
Input  Signal  Vector with N_{x} elements 


Window, g(n)  Second positional argument (Default: Hamming window) 
(Default: Periodic Hann window)  Kaiser window only  
Window length, M  Specified as number of samples (Default: ⌊N_{x}/4.5⌋)  Specified as number of samples (Default: 128)  TimeResolution namevalue argument  
Leakage, ℓ 

 Leakage namevalue argument related to Kaiser window
β shape factor: ℓ = 1 – β/40  
Overlap, L  Number of samples specified as third positional argument (Default: 50% of window length)  Number of samples specified with (Default: 75% of window length)  Percentage of segment length specified with (Default: $$\left(1\frac{\text{1}}{2\times \text{ENBW}1}\right)\times 100,$$ where ENBW is the equivalent noise bandwidth of the window)  
Number of DFT points, N_{DFT}  Fourth positional argument (Default: max(256,2^{⌈log2M⌉})) 
(Default: 128)  Always 1024  
Time information  Sample rate specified as fifth positional argument  Sample rate or time vector specified as second positional argument  Sample rate or time vector specified as second positional argument  
Function call 
fs = 100; x = exp(2j*pi*20* ... (0:1/fs:21/fs)); M = 200; lk = 0.5; g = kaiser(M,40*(1lk)); L = 100; Ndft = 1024; 
sps = abs( ... spectrogram(x,g,L, ... Ndft,fs,"centered") ... ).^2; 
sts = abs( ... stft(x,fs,Window=g, ... OverlapLength=L, ... FFTLength=Ndft) ... ).^2; 
pss = pspectrum(x,fs, ... "spectrogram", ... TimeResolution=M/fs, ... Leakage=lk, ... OverlapPercent=L/M*100 ... )*sum(g)^2; 
Convenience plot 
fs = 6e2; ts = 0:1/fs:2.05; x = vco(sin(2*pi*ts).* ... exp(ts),[0.1 0.4]*fs,fs); M = 32; lk = 0.9; g = kaiser(M,40*(1lk)); L = 22; Ndft = 1024;




Output  Frequency range 
 Controlled using
 Controlled using

Time interval 


 
Normalization 
 First output argument is STFT. Its magnitude squared is the spectrogram. 
 
PSD and power spectrum 
 Output argument is STFT. 
 
Examples 
STFT and Spectrogram Definitions
The STFT of a signal is computed by sliding an analysis window g(n) of length M over the signal and calculating the discrete Fourier transform (DFT) of each segment of windowed data. The window hops over the original signal at intervals of R samples, equivalent to L = M – R samples of overlap between adjoining segments. Most window functions taper off at the edges to avoid spectral ringing. The DFT of each windowed segment is added to a complexvalued matrix that contains the magnitude and phase for each point in time and frequency. The STFT matrix has
$$k=\lfloor \frac{{N}_{x}L}{ML}\rfloor $$
columns, where N_{x} is the length of the signal x(n) and the ⌊⌋ symbols denote the floor function. The number of rows in the matrix equals N_{DFT}, the number of DFT points, for centered and twosided transforms and an odd number close to N_{DFT}/2 for onesided transforms of realvalued signals.
The mth column of the STFT matrix $$X(f)=\left[\begin{array}{ccccc}{X}_{1}(f)& {X}_{2}(f)& {X}_{3}(f)& \cdots & {X}_{k}(f)\end{array}\right]$$ contains the DFT of the windowed data centered about time mR:
$${X}_{m}(f)={\displaystyle \sum _{n=\infty}^{\infty}x(n)\text{\hspace{0.17em}}g(nmR)\text{\hspace{0.17em}}{e}^{j2\pi fn}}.$$
Compare spectrogram
Function and STFT Definition
Generate a signal that consists of a complexvalued convex quadratic chirp sampled at 600 Hz for 2 seconds. The chirp has an initial frequency of 250 Hz and a final frequency of 50 Hz.
fs = 6e2; ts = 0:1/fs:2; x = chirp(ts,250,ts(end),50,"quadratic",0,"convex","complex");
spectrogram
Function
Use the spectrogram
function to compute the STFT of the signal.
Divide the signal into segments, each $\mathit{M}=49$ samples long.
Specify $\mathit{L}=11$ samples of overlap between adjoining segments.
Discard the final, shorter segment.
Window each segment with a Bartlett window.
Evaluate the discrete Fourier transform of each segment at ${\mathit{N}}_{\mathrm{DFT}}=1024$ points. By default,
spectrogram
computes twosided transforms for complexvalued signals.
M = 49; L = 11; g = bartlett(M); Ndft = 1024; [s,f,t] = spectrogram(x,g,L,Ndft,fs);
Use the waterplot
function to compute and display the spectrogram, defined as the magnitude squared of the STFT.
waterplot(s,f,t)
STFT Definition
Compute the STFT of the ${\mathit{N}}_{\mathit{x}}$sample signal using the definition. Divide the signal into $$\lfloor \frac{{N}_{x}L}{ML}\rfloor $$ overlapping segments. Window each segment and evaluate its discrete Fourier transform at ${\mathit{N}}_{\mathrm{DFT}}$ points.
[segs,~] = buffer(1:length(x),M,L,"nodelay");
X = fft(x(segs).*g,Ndft);
Compute the time and frequency ranges for the STFT.
To find the time values, divide the time vector into overlapping segments. The time values are the midpoints of the segments, with each segment treated as an interval open at the lower end.
To find the frequency values, specify a Nyquist interval closed at zero frequency and open at the lower end.
tbuf = ts(segs); tint = mean(tbuf(2:end,:)); fint = 0:fs/Ndft:fsfs/Ndft;
Compare the output of spectrogram
to the definition. Display the spectrogram.
maxdiff = max(max(abs(sX)))
maxdiff = 0
waterplot(X,fint,tint)
function waterplot(s,f,t) % Waterfall plot of spectrogram waterfall(f,t,abs(s)'.^2) set(gca,XDir="reverse",View=[30 50]) xlabel("Frequency (Hz)") ylabel("Time (s)") end
Compare spectrogram
and stft
Functions
Generate a signal consisting of a chirp sampled at 1.4 kHz for 2 seconds. The frequency of the chirp decreases linearly from 600 Hz to 100 Hz during the measurement time.
fs = 1400; x = chirp(0:1/fs:2,600,2,100);
stft
Defaults
Compute the STFT of the signal using the spectrogram
and stft
functions. Use the default values of the stft
function:
Divide the signal into 128sample segments and window each segment with a periodic Hann window.
Specify 96 samples of overlap between adjoining segments. This length is equivalent to 75% of the window length.
Specify 128 DFT points and center the STFT at zero frequency, with the frequency expressed in hertz.
Verify that the two results are equal.
M = 128; g = hann(M,"periodic"); L = 0.75*M; Ndft = 128; [sp,fp,tp] = spectrogram(x,g,L,Ndft,fs,"centered"); [s,f,t] = stft(x,fs); dff = max(max(abs(sps)))
dff = 0
Use the mesh
function to plot the two outputs.
nexttile mesh(tp,fp,abs(sp).^2) title("spectrogram") view(2), axis tight nexttile mesh(t,f,abs(s).^2) title("stft") view(2), axis tight
spectrogram
Defaults
Repeat the computation using the default values of the spectrogram
function:
Divide the signal into segments of length $\mathit{M}=\lfloor {\mathit{N}}_{\mathit{x}}/4.5\rfloor $, where ${\mathit{N}}_{\mathit{x}}$ is the length of the signal. Window each segment with a Hamming window.
Specify 50% overlap between segments.
To compute the FFT, use $\mathrm{max}\left(256,{2}^{\lceil {\mathrm{log}}_{2}\mathit{M}\rceil}\right)$ points. Compute the spectrogram only for positive normalized frequencies.
M = floor(length(x)/4.5); g = hamming(M); L = floor(M/2); Ndft = max(256,2^nextpow2(M)); [sx,fx,tx] = spectrogram(x); [st,ft,tt] = stft(x,Window=g,OverlapLength=L, ... FFTLength=Ndft,FrequencyRange="onesided"); dff = max(max(sxst))
dff = 0
Use the waterplot
function to plot the two outputs. Divide the frequency axis by $\pi $ in both cases. For the stft
output, divide the sample numbers by the effective sample rate, $2\pi $.
figure nexttile waterplot(sx,fx/pi,tx) title("spectrogram") nexttile waterplot(st,ft/pi,tt/(2*pi)) title("stft")
function waterplot(s,f,t) % Waterfall plot of spectrogram waterfall(f,t,abs(s)'.^2) set(gca,XDir="reverse",View=[30 50]) xlabel("Frequency/\pi") ylabel("Samples") end
Compare spectrogram
and pspectrum
Functions
Generate a signal that consists of a voltagecontrolled oscillator and three Gaussian atoms. The signal is sampled at ${\mathit{f}}_{\mathit{s}}=2$ kHz for 2 seconds.
fs = 2000; tx = 0:1/fs:2; gaussFun = @(A,x,mu,f) exp((xmu).^2/(2*0.03^2)).*sin(2*pi*f.*x)*A'; s = gaussFun([1 1 1],tx',[0.1 0.65 1],[2 6 2]*100)*1.5; x = vco(chirp(tx+.1,0,tx(end),3).*exp(2*(tx1).^2),[0.1 0.4]*fs,fs); x = s+x';
ShortTime Fourier Transforms
Use the pspectrum
function to compute the STFT.
Divide the ${\mathit{N}}_{\mathit{x}}$sample signal into segments of length $\mathit{M}=80$ samples, corresponding to a time resolution of $80/2000=40$ milliseconds.
Specify $\mathit{L}=16$ samples or 20% of overlap between adjoining segments.
Window each segment with a Kaiser window and specify a leakage $\ell =0.7$.
M = 80; L = 16; lk = 0.7; [S,F,T] = pspectrum(x,fs,"spectrogram", ... TimeResolution=M/fs,OverlapPercent=L/M*100, ... Leakage=lk);
Compare to the result obtained with the spectrogram
function.
Specify the window length and overlap directly in samples.
pspectrum
always uses a Kaiser window as $\mathit{g}\left(\mathit{n}\right)$. The leakage $\ell $ and the shape factor $\beta $ of the window are related by $\beta =40\times \left(1\ell \right)$.pspectrum
always uses ${\mathit{N}}_{\mathrm{DFT}}=1024$ points when computing the discrete Fourier transform. You can specify this number if you want to compute the transform over a twosided or centered frequency range. However, for onesided transforms, which are the default for real signals,spectrogram
uses $1024/2+1=513$ points. Alternatively, you can specify the vector of frequencies at which you want to compute the transform, as in this example.If a signal cannot be divided exactly into $\mathit{k}=\lfloor \frac{{\mathit{N}}_{\mathit{x}}\mathit{L}}{\mathit{M}\mathit{L}}\rfloor $ segments,
spectrogram
truncates the signal whereaspspectrum
pads the signal with zeros to create an extra segment. To make the outputs equivalent, remove the final segment and the final element of the time vector.spectrogram
returns the STFT, whose magnitude squared is the spectrogram.pspectrum
returns the segmentbysegment power spectrum, which is already squared but is divided by a factor of ${\sum}_{\mathit{n}}\mathit{g}\left(\mathit{n}\right)$ before squaring.For onesided transforms,
pspectrum
adds an extra factor of 2 to the spectrogram.
g = kaiser(M,40*(1lk)); k = (length(x)L)/(ML); if k~=floor(k) S = S(:,1:floor(k)); T = T(1:floor(k)); end [s,f,t] = spectrogram(x/sum(g)*sqrt(2),g,L,F,fs);
Use the waterplot
function to display the spectrograms computed by the two functions.
subplot(2,1,1) waterplot(sqrt(S),F,T) title("pspectrum") subplot(2,1,2) waterplot(s,f,t) title("spectrogram")
maxd = max(max(abs(abs(s).^2S)))
maxd = 2.4419e08
Power Spectra and Convenience Plots
The spectrogram
function has a fourth argument that corresponds to the segmentbysegment power spectrum or power spectral density. Similar to the output of pspectrum
, the ps
argument is already squared and includes the normalization factor ${\sum}_{\mathit{n}}\mathit{g}\left(\mathit{n}\right)$. For onesided spectrograms of real signals, you still have to include the extra factor of 2. Set the scaling argument of the function to "power"
.
[~,~,~,ps] = spectrogram(x*sqrt(2),g,L,F,fs,"power");
max(abs(S(:)ps(:)))
ans = 2.4419e08
When called with no output arguments, both pspectrum
and spectrogram
plot the spectrogram of the signal in decibels. Include the factor of 2 for onesided spectrograms. Set the colormaps to be the same for both plots. Set the xlimits to the same values to make visible the extra segment at the end of the pspectrum
plot. In the spectrogram
plot, display the frequency on the yaxis.
subplot(2,1,1) pspectrum(x,fs,"spectrogram", ... TimeResolution=M/fs,OverlapPercent=L/M*100, ... Leakage=lk) title("pspectrum") cc = clim; xl = xlim; subplot(2,1,2) spectrogram(x*sqrt(2),g,L,F,fs,"power","yaxis") title("spectrogram") clim(cc) xlim(xl)
function waterplot(s,f,t) % Waterfall plot of spectrogram waterfall(f,t,abs(s)'.^2) set(gca,XDir="reverse",View=[30 50]) xlabel("Frequency (Hz)") ylabel("Time (s)") end
Compute Centered and OneSided Spectrograms
Generate a signal that consists of a realvalued chirp sampled at 2 kHz for 2 seconds.
fs = 2000;
tx = 0:1/fs:2;
x = vco(chirp(tx,0,tx(end),2).*exp(3*(tx1).^2), ...
[0.1 0.4]*fs,fs).*hann(length(tx))';
TwoSided Spectrogram
Compute and plot the twosided STFT of the signal.
Divide the signal into segments, each $\mathit{M}=73$ samples long.
Specify $\mathit{L}=24$ samples of overlap between adjoining segments.
Discard the final, shorter segment.
Window each segment with a flattop window.
Evaluate the discrete Fourier transform of each segment at ${\mathit{N}}_{\mathrm{DFT}}=895$ points, noting that it is an odd number.
M = 73;
L = 24;
g = flattopwin(M);
Ndft = 895;
neven = ~mod(Ndft,2);
[stwo,f,t] = spectrogram(x,g,L,Ndft,fs,"twosided");
Use the spectrogram
function with no output arguments to plot the twosided spectrogram.
spectrogram(x,g,L,Ndft,fs,"twosided","power","yaxis");
Compute the twosided spectrogram using the definition. Divide the signal into $\mathit{M}$sample segments with $\mathit{L}$ samples of overlap between adjoining segments. Window each segment and compute its discrete Fourier transform at ${\mathit{N}}_{\mathrm{DFT}}$ points.
[segs,~] = buffer(1:length(x),M,L,"nodelay");
Xtwo = fft(x(segs).*g,Ndft);
Compute the time and frequency ranges.
To find the time values, divide the time vector into overlapping segments. The time values are the midpoints of the segments, with each segment treated as an interval open at the lower end.
To find the frequency values, specify a Nyquist interval closed at zero frequency and open at the upper end.
tbuf = tx(segs); ttwo = mean(tbuf(2:end,:)); ftwo = 0:fs/Ndft:fs*(11/Ndft);
Compare the outputs of spectrogram
to the definitions. Use the waterplot
function to display the spectrograms.
diffs = [max(max(abs(stwoXtwo))) max(abs(fftwo')) max(abs(tttwo))]
diffs = 1×3
10^{12} ×
0 0.2274 0.0002
figure nexttile waterplot(Xtwo,ftwo,ttwo) title("TwoSided, Definition") nexttile waterplot(stwo,f,t) title("TwoSided, spectrogram Function")
Centered Spectrogram
Compute the centered spectrogram of the signal.
Use the same time values that you used for the twosided STFT.
Use the
fftshift
function to shift the zerofrequency component of the STFT to the center of the spectrum.For oddvalued ${\mathit{N}}_{\mathrm{DFT}}$, the frequency interval is open at both ends. For evenvalued ${\mathit{N}}_{\mathrm{DFT}}$, the frequency interval is open at the lower end and closed at the upper end.
Compare the outputs and display the spectrograms.
tcen = ttwo; if ~neven Xcen = fftshift(Xtwo,1); fcen = fs/2*(11/Ndft):fs/Ndft:fs/2; else Xcen = fftshift(circshift(Xtwo,1),1); fcen = (fs/2*(11/Ndft):fs/Ndft:fs/2)+fs/Ndft/2; end [scen,f,t] = spectrogram(x,g,L,Ndft,fs,"centered"); diffs = [max(max(abs(scenXcen))) max(abs(ffcen')) max(abs(ttcen))]
diffs = 1×3
10^{12} ×
0 0.2274 0.0002
figure nexttile waterplot(Xcen,fcen,tcen) title("Centered, Definition") nexttile waterplot(scen,f,t) title("Centered, spectrogram Function")
OneSided Spectrogram
Compute the onesided spectrogram of the signal.
Use the same time values that you used for the twosided STFT.
For oddvalued ${\mathit{N}}_{\mathrm{DFT}}$, the onesided STFT consists of the first $\left({\mathit{N}}_{\mathrm{DFT}}+1\right)/2$ rows of the twosided STFT. For evenvalued ${\mathit{N}}_{\mathrm{DFT}}$, the onesided STFT consists of the first ${\mathit{N}}_{\mathrm{DFT}}/2+1$ rows of the twosided STFT.
For oddvalued ${\mathit{N}}_{\mathrm{DFT}}$, the frequency interval is closed at zero frequency and open at the Nyquist frequency. For evenvalued ${\mathit{N}}_{\mathrm{DFT}}$, the frequency interval is closed at both ends.
Compare the outputs and display the spectrograms. For realvalued signals, the "onesided"
argument is optional.
tone = ttwo; if ~neven Xone = Xtwo(1:(Ndft+1)/2,:); else Xone = Xtwo(1:Ndft/2+1,:); end fone = 0:fs/Ndft:fs/2; [sone,f,t] = spectrogram(x,g,L,Ndft,fs); diffs = [max(max(abs(soneXone))) max(abs(ffone')) max(abs(ttone))]
diffs = 1×3
10^{12} ×
0 0.1137 0.0002
figure nexttile waterplot(Xone,fone,tone) title("OneSided, Definition") nexttile waterplot(sone,f,t) title("OneSided, spectrogram Function")
function waterplot(s,f,t) % Waterfall plot of spectrogram waterfall(f,t,abs(s)'.^2) set(gca,XDir="reverse",View=[30 50]) xlabel("Frequency (Hz)") ylabel("Time (s)") end
Compute Segment PSDs and Power Spectra
The spectrogram
function has a matrix containing either the power spectral density (PSD) or the power spectrum of each segment as the fourth output argument. The power spectrum is equal to the PSD multiplied by the equivalent noise bandwidth (ENBW) of the window.
Generate a signal that consists of a logarithmic chirp sampled at 1 kHz for 1 second. The chirp has an initial frequency of 400 Hz that decreases to 10 Hz by the end of the measurement.
fs = 1000;
tt = 0:1/fs:11/fs;
y = chirp(tt,400,tt(end),10,"logarithmic");
Segment PSDs and Power Spectra with Sample Rate
Divide the signal into 102sample segments and window each segment with a Hann window. Specify 12 samples of overlap between adjoining segments and 1024 DFT points.
M = 102; g = hann(M); L = 12; Ndft = 1024;
Compute the spectrogram of the signal with the default PSD spectrum type. Output the STFT and the array of segment power spectral densities.
[s,f,t,p] = spectrogram(y,g,L,Ndft,fs);
Repeat the computation with the spectrum type specified as "power"
. Output the STFT and the array of segment power spectra.
[r,~,~,q] = spectrogram(y,g,L,Ndft,fs,"power");
Verify that the spectrogram is the same in both cases. Plot the spectrogram using a logarithmic scale for the frequency.
max(max(abs(s).^2abs(r).^2))
ans = 0
waterfall(f,t,abs(s)'.^2) set(gca,XScale="log",... XDir="reverse",View=[30 50])
Verify that the power spectra are equal to the power spectral densities multiplied by the ENBW of the window.
max(max(abs(qp*enbw(g,fs))))
ans = 1.1102e16
Verify that the matrix of segment power spectra is proportional to the spectrogram. The proportionality factor is the square of the sum of the window elements.
max(max(abs(s).^2q*sum(g)^2))
ans = 3.4694e18
Segment PSDs and Power Spectra with Normalized Frequencies
Repeat the computation, but now work in normalized frequencies. The results are the same when you specify the sample rate as $2\pi $.
[~,~,~,pn] = spectrogram(y,g,L,Ndft);
[~,~,~,qn] = spectrogram(y,g,L,Ndft,"power");
max(max(abs(qnpn*enbw(g,2*pi))))
ans = 1.1102e16
See Also
Apps
Functions
pspectrum
spectrogram
stft
istft
xspectrogram