pspectrum

Analyze signals in the frequency and time-frequency domains

Syntax

p = pspectrum(x)

p = pspectrum(x,fs)

p = pspectrum(x,t)

p = pspectrum(___,type)

p = pspectrum(___,Name,Value)

[p,f] = pspectrum(___)

[p,f,t] = pspectrum(___,'spectrogram')

[p,f,pwr] = pspectrum(___,'persistence')

pspectrum(___)

Description

p = pspectrum(x) returns the power spectrum of x.

• If x is a vector or a timetable with a vector of data, then it is treated as a single channel.

• If x is a matrix, a timetable with a matrix variable, or a timetable with multiple vector variables, then the spectrum is computed independently for each channel and stored in a separate column of p.

p = pspectrum(x,fs) returns the power spectrum of a vector or matrix signal sampled at a rate fs.

p = pspectrum(x,t) returns the power spectrum of a vector or matrix signal sampled at the time instants specified in t.

p = pspectrum(___,type) specifies the kind of spectral analysis performed by the function. Specify type as 'power', 'spectrogram', or 'persistence'. This syntax can include any combination of input arguments from previous syntaxes.

p = pspectrum(___,Name,Value) specifies additional options using name-value pair arguments. Options include the frequency resolution bandwidth and the percent overlap between adjoining segments.

[p,f] = pspectrum(___) returns the frequencies corresponding to the spectral estimates contained in p.

[p,f,t] = pspectrum(___,'spectrogram') also returns a vector of time instants corresponding to the centers of the windowed segments used to compute short-time power spectrum estimates.

[p,f,pwr] = pspectrum(___,'persistence') also returns a vector of power values corresponding to the estimates contained in a persistence spectrum.

pspectrum(___) with no output arguments plots the spectral estimate in the current figure window. For the plot, the function converts p to dB using 10 log₁₀(p).

Examples

Power Spectra of Sinusoids

Generate 128 samples of a two-channel complex sinusoid.

• The first channel has unit amplitude and a normalized sinusoid frequency of $\pi /4$ rad/sample

• The second channel has an amplitude of $1/\sqrt{2}$and a normalized frequency of $\pi /2$ rad/sample.

Compute and plot the power spectrum of each channel. Zoom in on the frequency range from $0.15\pi$ rad/sample to $0.6\pi$ rad/sample. pspectrum scales the spectrum so that, if the frequency content of a signal falls exactly within a bin, its amplitude in that bin is the true average power of the signal. For a complex exponential, the average power is the square of the amplitude. Verify by computing the discrete Fourier transform of the signal. For more details, see Measure Power of Deterministic Periodic Signals.

N = 128;
x = [1 1/sqrt(2)].*exp(1j*pi./[4;2]*(0:N-1)).';

[p,f] = pspectrum(x);

plot(f/pi,p)
hold on
stem(0:2/N:2-1/N,abs(fft(x)/N).^2)
hold off
axis([0.15 0.6 0 1.1])
legend("Channel "+[1;2]+", "+["pspectrum" "fft"])
grid

Generate a sinusoidal signal sampled at 1 kHz for 296 milliseconds and embedded in white Gaussian noise. Specify a sinusoid frequency of 200 Hz and a noise variance of 0.1². Store the signal and its time information in a MATLAB® timetable.

Fs = 1000;
t = (0:1/Fs:0.296)';
x = cos(2*pi*t*200)+0.1*randn(size(t));
xTable = timetable(seconds(t),x);

Compute the power spectrum of the signal. Express the spectrum in decibels and plot it.

[pxx,f] = pspectrum(xTable);

plot(f,pow2db(pxx))
grid on
xlabel('Frequency (Hz)')
ylabel('Power Spectrum (dB)')
title('Default Frequency Resolution')

Recompute the power spectrum of the sinusoid, but now use a coarser frequency resolution of 25 Hz. Plot the spectrum using the pspectrum function with no output arguments.

pspectrum(xTable,'FrequencyResolution',25)

Two-Sided Spectra

Generate a signal sampled at 3 kHz for 1 second. The signal is a convex quadratic chirp whose frequency increases from 300 Hz to 1300 Hz during the measurement. The chirp is embedded in white Gaussian noise.

fs = 3000;
t = 0:1/fs:1-1/fs;

x1 = chirp(t,300,t(end),1300,'quadratic',0,'convex') + ...
    randn(size(t))/100;

Compute and plot the two-sided power spectrum of the signal using a rectangular window. For real signals, pspectrum plots a one-sided spectrum by default. To plot a two-sided spectrum, set TwoSided to true.

pspectrum(x1,fs,'Leakage',1,'TwoSided',true)

Generate a complex-valued signal with the same duration and sample rate. The signal is a chirp with sinusoidally varying frequency content and embedded in white noise. Compute the spectrogram of the signal and display it as a waterfall plot. For complex-valued signals, the spectrogram is two-sided by default.

x2 = exp(2j*pi*100*cos(2*pi*2*t)) + randn(size(t))/100;

[p,f,t] = pspectrum(x2,fs,'spectrogram');

waterfall(f,t,p')
xlabel('Frequency (Hz)')
ylabel('Time (seconds)')
wtf = gca;
wtf.XDir = 'reverse';
view([30 45])

Window Leakage and Tone Resolution

Generate a two-channel signal sampled at 100 Hz for 2 seconds.

1. The first channel consists of a 20 Hz tone and a 21 Hz tone. Both tones have unit amplitude.

2. The second channel also has two tones. One tone has unit amplitude and a frequency of 20 Hz. The other tone has an amplitude of 1/100 and a frequency of 30 Hz.

fs = 100;
t = (0:1/fs:2-1/fs)';

x = sin(2*pi*[20 20].*t) + [1 1/100].*sin(2*pi*[21 30].*t);

Embed the signal in white noise. Specify a signal-to-noise ratio of 40 dB. Plot the signals.

x = x + randn(size(x)).*std(x)/db2mag(40);

plot(t,x)

Compute the spectra of the two channels and display them.

pspectrum(x,t)

The default value for the spectral leakage, 0.5, corresponds to a resolution bandwidth of about 1.29 Hz. The two tones in the first channel are not resolved. The 30 Hz tone in the second channel is visible, despite being much weaker than the other one.

Increase the leakage to 0.85, equivalent to a resolution of about 0.74 Hz. The weak tone in the second channel is clearly visible.

pspectrum(x,t,'Leakage',0.85)

Increase the leakage to the maximum value. The resolution bandwidth is approximately 0.5 Hz. The two tones in the first channel are resolved. The weak tone in the second channel is masked by the large window sidelobes.

pspectrum(x,t,'Leakage',1)

Compare spectrogram and pspectrum Functions

Generate a signal that consists of a voltage-controlled 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(-(x-mu).^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*(tx-1).^2),[0.1 0.4]*fs,fs);
x = s+x';

Short-Time 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 two-sided or centered frequency range. However, for one-sided 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 whereas pspectrum 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 segment-by-segment 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 one-sided transforms, pspectrum adds an extra factor of 2 to the spectrogram.

g = kaiser(M,40*(1-lk));

k = (length(x)-L)/(M-L);
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).^2-S)))

maxd = 
2.4419e-08

Power Spectra and Convenience Plots

The spectrogram function has a fourth argument that corresponds to the segment-by-segment 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 one-sided 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.4419e-08

When called with no output arguments, both pspectrum and spectrogram plot the spectrogram of the signal in decibels. Include the factor of 2 for one-sided spectrograms. Set the colormaps to be the same for both plots. Set the x-limits 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 y-axis.

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

Persistence Spectrum of Transient Signal

Visualize an interference narrowband signal embedded within a broadband signal.

Generate a chirp sampled at 1 kHz for 500 seconds. The frequency of the chirp increases from 180 Hz to 220 Hz during the measurement.

fs = 1000;
t = (0:1/fs:500)';

x = chirp(t,180,t(end),220) + 0.15*randn(size(t));

The signal also contains a 210 Hz sinusoid. The sinusoid has an amplitude of 0.05 and is present only for 1/6 of the total signal duration.

idx = floor(length(x)/6);
x(1:idx) = x(1:idx) + 0.05*cos(2*pi*t(1:idx)*210);

Compute the spectrogram of the signal. Restrict the frequency range from 100 Hz to 290 Hz. Specify a time resolution of 1 second. Both signal components are visible.

pspectrum(x,fs,'spectrogram', ...
    'FrequencyLimits',[100 290],'TimeResolution',1)

Compute the power spectrum of the signal. The weak sinusoid is obscured by the chirp.

pspectrum(x,fs,'FrequencyLimits',[100 290])

Compute the persistence spectrum of the signal. Now both signal components are clearly visible.

pspectrum(x,fs,'persistence', ...
    'FrequencyLimits',[100 290],'TimeResolution',1)

Spectrogram and Reassigned Spectrogram of Chirp

Generate a quadratic chirp sampled at 1 kHz for 2 seconds. The chirp has an initial frequency of 100 Hz that increases to 200 Hz at t = 1 second. Compute the spectrogram using the default settings of the pspectrum function. Use the waterfall function to plot the spectrogram.

fs = 1e3;
t = 0:1/fs:2;
y = chirp(t,100,1,200,"quadratic");

[sp,fp,tp] = pspectrum(y,fs,"spectrogram");

waterfall(fp,tp,sp')
set(gca,XDir="reverse",View=[60 60])
ylabel("Time (s)")
xlabel("Frequency (Hz)")

Compute and display the reassigned spectrogram.

[sr,fr,tr] = pspectrum(y,fs,"spectrogram",Reassign=true);

waterfall(fr,tr,sr')
set(gca,XDir="reverse",View=[60 60])
ylabel("Time (s)")
xlabel("Frequency (Hz)")

Recompute the spectrogram using a time resolution of 0.2 second. Visualize the result using the pspectrum function with no output arguments.

pspectrum(y,fs,"spectrogram",TimeResolution=0.2)

Compute the reassigned spectrogram using the same time resolution.

pspectrum(y,fs,"spectrogram",TimeResolution=0.2,Reassign=true)

Spectrogram of Dial Tone Signal

Create a signal, sampled at 4 kHz, that resembles pressing all the keys of a digital telephone. Save the signal as a MATLAB® timetable.

fs = 4e3;
t = 0:1/fs:0.5-1/fs;

ver = [697 770 852 941];
hor = [1209 1336 1477];

tones = [];

for k = 1:length(ver)
    for l = 1:length(hor)
        tone = sum(sin(2*pi*[ver(k);hor(l)].*t))';
        tones = [tones;tone;zeros(size(tone))];
    end
end

% To hear, type soundsc(tones,fs)

S = timetable(seconds(0:length(tones)-1)'/fs,tones);

Compute the spectrogram of the signal. Specify a time resolution of 0.5 second and zero overlap between adjoining segments. Specify the leakage as 0.85, which is approximately equivalent to windowing the data with a Hann window.

pspectrum(S,'spectrogram', ...
    'TimeResolution',0.5,'OverlapPercent',0,'Leakage',0.85)

The spectrogram shows that each key is pressed for half a second, with half-second silent pauses between keys. The first tone has a frequency content concentrated around 697 Hz and 1209 Hz, corresponding to the digit '1' in the DTMF standard.

Input Arguments

x — Input signal

vector | matrix | timetable

Input signal, specified as a vector, a matrix, or a MATLAB^® timetable.

• If x is a timetable, then it must contain increasing finite row times.

  Note

  If a timetable has missing or duplicate time points, you can fix it using the tips in Clean Timetable with Missing, Duplicate, or Nonuniform Times.

• If x is a timetable representing a multichannel signal, then it must have either a single variable containing a matrix or multiple variables consisting of vectors.

If x is nonuniformly sampled, then pspectrum interpolates the signal to a uniform grid to compute spectral estimates. The function uses linear interpolation and assumes a sample time equal to the median of the differences between adjacent time points. For a nonuniformly sampled signal to be supported, the median time interval and the mean time interval must obey

$$\frac{1}{100}<\frac{\text{Median time interval}}{\text{Mean time interval}}<100.$$

Example: cos(pi./[4;2]*(0:159))'+randn(160,2) is a two-channel signal consisting of sinusoids embedded in white noise.

Example: timetable(seconds(0:4)',rand(5,2)) specifies a two-channel random variable sampled at 1 Hz for 4 seconds.

Example: timetable(seconds(0:4)',rand(5,1),rand(5,1)) specifies a two-channel random variable sampled at 1 Hz for 4 seconds.

Data Types: single | double
Complex Number Support: Yes

fs — Sample rate

2π (default) | positive numeric scalar

Sample rate, specified as a positive numeric scalar.

t — Time values

vector | datetime array | duration array | duration scalar

Time values, specified as a vector, a datetime or duration array, or a duration scalar representing the time interval between samples.

Example: seconds(0:1/100:1) is a duration array representing 1 second of sampling at 100 Hz.

Example: seconds(1) is a duration scalar representing a 1-second time difference between consecutive signal samples.

type — Type of spectrum to compute

'power' (default) | 'spectrogram' | 'persistence'

Type of spectrum to compute, specified as 'power', 'spectrogram', or 'persistence':

• 'power' — Compute the power spectrum of the input. Use this option to analyze the frequency content of a stationary signal. For more information, Spectrum Computation.

• 'spectrogram' — Compute the spectrogram of the input. Use this option to analyze how the frequency content of a signal changes over time. For more information, see Spectrogram Computation.

• 'persistence' — Compute the persistence power spectrum of the input. Use this option to visualize the fraction of time that a particular frequency component is present in a signal. For more information, see Persistence Spectrum Computation.

Note

The 'spectrogram' and 'persistence' options do not support multichannel input.

Name-Value Arguments

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Before R2021a, use commas to separate each name and value, and enclose Name in quotes.

Example: 'Leakage',1,'Reassign',true,'MinThreshold',-35 windows the data using a rectangular window, computes a reassigned spectrum estimate, and sets all values smaller than –35 dB to zero.

FrequencyLimits — Frequency band limits

[0 fs/2] (default) | two-element numeric vector

Frequency band limits, specified as the comma-separated pair consisting of 'FrequencyLimits' and a two-element numeric vector:

• If the input contains time information, then the frequency band is expressed in Hz.

• If the input does not contain time information, then the frequency band is expressed in normalized units of rad/sample.

By default, pspectrum computes the spectrum over the whole Nyquist range:

• If the specified frequency band contains a region that falls outside the Nyquist range, then pspectrum truncates the frequency band.

• If the specified frequency band lies completely outside of the Nyquist range, then pspectrum throws an error.

See Spectrum Computation for more information about the Nyquist range.

If x is nonuniformly sampled, then pspectrum linearly interpolates the signal to a uniform grid and defines an effective sample rate equal to the inverse of the median of the differences between adjacent time points. Express 'FrequencyLimits' in terms of the effective sample rate.

Example: [0.2*pi 0.7*pi] computes the spectrum between 0.2π rad/sample and 0.7π rad/sample of a signal with no time information.

FrequencyResolution — Frequency resolution bandwidth

real numeric scalar

Frequency resolution bandwidth, specified as the comma-separated pair consisting of 'FrequencyResolution' and a real numeric scalar, expressed in Hz if the input contains time information, or in normalized units of rad/sample if not. This argument cannot be specified simultaneously with 'TimeResolution'. The default value of this argument depends on the size of the input data. See Spectrogram Computation for details.

Example: pi/100 computes the spectrum of a signal with no time information using a frequency resolution of π/100 rad/sample.

Leakage — Spectral leakage

0.5 (default) | real numeric scalar between 0 and 1

Spectral leakage, specified as the comma-separated pair consisting of 'Leakage' and a real numeric scalar between 0 and 1. 'Leakage' controls the Kaiser window sidelobe attenuation relative to the mainlobe width, compromising between improving resolution and decreasing leakage:

• A large leakage value resolves closely spaced tones, but masks nearby weak tones.

• A small leakage value finds small tones in the vicinity of larger tones, but smears close frequencies together.

Example: 'Leakage',0 reduces leakage to a minimum at the expense of spectral resolution.

Example: 'Leakage',0.85 approximates windowing the data with a Hann window.

Example: 'Leakage',1 is equivalent to windowing the data with a rectangular window, maximizing leakage but improving spectral resolution.

MinThreshold — Lower bound for nonzero values

-Inf (default) | real scalar

Lower bound for nonzero values, specified as the comma-separated pair consisting of 'MinThreshold' and a real scalar. pspectrum implements 'MinThreshold' differently based on the value of the type argument:

• 'power' or 'spectrogram' — pspectrum sets those elements of p such that 10 log₁₀(p) ≤ 'MinThreshold' to zero. Specify 'MinThreshold' in decibels.

• 'persistence' — pspectrum sets those elements of p smaller than 'MinThreshold' to zero. Specify 'MinThreshold' between 0 and 100%.

NumPowerBins — Number of power bins for persistence spectrum

256 (default) | integer between 20 and 1024

Number of power bins for persistence spectrum, specified as the comma-separated pair consisting of 'NumPowerBins' and an integer between 20 and 1024.

OverlapPercent — Overlap between adjoining segments

real scalar in the interval [0, 100)

Overlap between adjoining segments for spectrogram or persistence spectrum, specified as the comma-separated pair consisting of 'OverlapPercent' and a real scalar in the interval [0, 100). The default value of this argument depends on the spectral window. See Spectrogram Computation for details.

Reassign — Reassignment option

false (default) | true

Reassignment option, specified as the comma-separated pair consisting of 'Reassign' and a logical value. If this option is set to true, then pspectrum sharpens the localization of spectral estimates by performing time and frequency reassignment. The reassignment technique produces periodograms and spectrograms that are easier to read and interpret. This technique reassigns each spectral estimate to the center of energy of its bin instead of the bin's geometric center. The technique provides exact localization for chirps and impulses.

TimeResolution — Time resolution of spectrogram or persistence spectrum

real scalar

Time resolution of spectrogram or persistence spectrum, specified as the comma-separated pair consisting of 'TimeResolution' and a real scalar, expressed in seconds if the input contains time information, or as an integer number of samples if not. This argument controls the duration of the segments used to compute the short-time power spectra that form spectrogram or persistence spectrum estimates. 'TimeResolution' cannot be specified simultaneously with 'FrequencyResolution'. The default value of this argument depends on the size of the input data and, if it was specified, the frequency resolution. See Spectrogram Computation for details.

TwoSided — Two-sided spectral estimate

false | true

Two-sided spectral estimate, specified as the comma-separated pair consisting of 'TwoSided' and a logical value.

• If this option is true, the function computes centered, two-sided spectrum estimates over [–π, π]. If the input has time information, the estimates are computed over [–f_s/2, f_s/2], where f_s is the effective sample rate.

• If this option is false, the function computes one-sided spectrum estimates over the Nyquist range [0, π]. If the input has time information, the estimates are computed over [0, f_s/2], where f_s is the effective sample rate. To conserve the total power, the function multiplies the power by 2 at all frequencies except 0 and the Nyquist frequency. This option is valid only for real signals.

If not specified, 'TwoSided' defaults to false for real input signals and to true for complex input signals.

Output Arguments

p — Spectrum

vector | matrix

Spectrum, returned as a vector or a matrix. The type and size of the spectrum depends on the value of the type argument:

• 'power' — p contains the power spectrum estimate of each channel of x. In this case, p is of size N_f × N_ch, where N_f is the length of f and N_ch is the number of channels of x. pspectrum scales the spectrum so that, if the frequency content of a signal falls exactly within a bin, its amplitude in that bin is the true average power of the signal. For example, the average power of a sinusoid is one-half the square of the sinusoid amplitude. For more details, see Measure Power of Deterministic Periodic Signals.

• 'spectrogram' — p contains an estimate of the short-term, time-localized power spectrum of x. In this case, p is of size N_f × N_t, where N_f is the length of f and N_t is the length of t.

• 'persistence' — p contains, expressed as percentages, the probabilities that the signal has components of a given power level at a given time and frequency location. In this case, p is of size N_pwr × N_f, where N_pwr is the length of pwr and N_f is the length of f.

f — Spectrum frequencies

vector

Spectrum frequencies, returned as a vector. If the input signal contains time information, then f contains frequencies expressed in Hz. If the input signal does not contain time information, then the frequencies are in normalized units of rad/sample.

t — Time values of spectrogram

vector | datetime array | duration array

Time values of spectrogram, returned as a vector of time values in seconds or a duration array. If the input does not have time information, then t contains sample numbers. t contains the time values corresponding to the centers of the data segments used to compute short-time power spectrum estimates.

• If the input to pspectrum is a timetable, then t has the same format as the time values of the input timetable.

• If the input to pspectrum is a numeric vector sampled at a set of time instants specified by a numeric, duration, or datetime array, then t has the same type and format as the input time values.

• If the input to pspectrum is a numeric vector with a specified time difference between consecutive samples, then t is a duration array.

pwr — Power values of persistence spectrum

vector

Power values of persistence spectrum, returned as a vector.

More About

Spectrum Computation

To compute signal spectra, pspectrum finds a compromise between the spectral resolution achievable with the entire length of the signal and the performance limitations that result from computing large FFTs:

• If possible, the function computes a single modified periodogram of the whole signal using a Kaiser window.

• If it is not possible to compute a single modified periodogram in a reasonable amount of time, the function computes a Welch periodogram: It divides the signal into overlapping segments, windows each segment using a Kaiser window, and averages the periodograms of the segments.

Spectral Windowing

Any real-world signal is measurable only for a finite length of time. This fact introduces nonnegligible effects into Fourier analysis, which assumes that signals are either periodic or infinitely long. Spectral windowing, which assigns different weights to different signal samples, deals systematically with finite-size effects.

The simplest way to window a signal is to assume that it is identically zero outside of the measurement interval and that all samples are equally significant. This "rectangular window" has discontinuous jumps at both ends that result in spectral ringing. All other spectral windows taper at both ends to lessen this effect by assigning smaller weights to samples close to the signal edges.

The windowing process always involves a compromise between conflicting aims: improving resolution and decreasing leakage:

• Resolution is the ability to know precisely how the signal energy is distributed in the frequency space. A spectrum analyzer with ideal resolution can distinguish two different tones (pure sinusoids) present in the signal, no matter how close in frequency. Quantitatively, this ability relates to the mainlobe width of the transform of the window.

• Leakage is the fact that, in a finite signal, every frequency component projects energy content throughout the complete frequency span. The amount of leakage in a spectrum can be measured by the ability to detect a weak tone from noise in the presence of a neighboring strong tone. Quantitatively, this ability relates to the sidelobe level of the frequency transform of the window.

• The spectrum is normalized so that a pure tone within that bandwidth, if perfectly centered, has the correct amplitude.

The better the resolution, the higher the leakage, and vice versa. At one end of the range, a rectangular window has the narrowest possible mainlobe and the highest sidelobes. This window can resolve closely spaced tones if they have similar energy content, but it fails to find the weaker one if they do not. At the other end, a window with high sidelobe suppression has a wide mainlobe in which close frequencies are smeared together.

pspectrum uses Kaiser windows to carry out windowing. For Kaiser windows, the fraction of the signal energy captured by the mainlobe depends most importantly on an adjustable shape factor, β. pspectrum uses shape factors ranging from β = 0, which corresponds to a rectangular window, to β = 40, where a wide mainlobe captures essentially all the spectral energy representable in double precision. An intermediate value of β ≈ 6 approximates a Hann window quite closely. To control β, use the 'Leakage' name-value pair. If you set 'Leakage' to ℓ, then ℓ and β are related by β = 40(1 – ℓ). See kaiser for more details.

51-point Hann window and 51-point Kaiser window with β = 5.7 in the time domain	51-point Hann window and 51-point Kaiser window with β = 5.7 in the frequency domain

Parameter and Algorithm Selection

To compute signal spectra, pspectrum initially determines the resolution bandwidth, which measures how close two tones can be and still be resolved. The resolution bandwidth has a theoretical value of

$${\text{RBW}}_{\text{theory}}=\frac{\text{ENBW}}{t_{\max }-t_{\min }}.$$

• t_max – t_min, the record length, is the time-domain duration of the selected signal region.

• ENBW is the equivalent noise bandwidth of the spectral window. See enbw for more details.

  Use the 'Leakage' name-value pair to control the ENBW. The minimum value of the argument corresponds to a Kaiser window with β = 40. The maximum value corresponds to a Kaiser window with β = 0.

In practice, however, pspectrum might lower the resolution. Lowering the resolution makes it possible to compute the spectrum in a reasonable amount of time and to display it with a finite number of pixels. For these practical reasons, the lowest resolution bandwidth pspectrum can use is

$${\text{RBW}}_{\text{performance}}=\frac{f_{\text{span}}}{4096-1},$$

where f_span is the width of the frequency band specified using 'FrequencyLimits'. If 'FrequencyLimits' is not specified, then pspectrum uses the sample rate as f_span. RBW_performance cannot be adjusted.

To compute the spectrum of a signal, the function chooses the larger of the two values, called the target resolution bandwidth:

$$\text{RBW}=\max ({\text{RBW}}_{\text{theory}},{\text{RBW}}_{\text{performance}}).$$

• If the resolution bandwidth is RBW_theory, then pspectrum computes a single modified periodogram for the whole signal. The function uses a Kaiser window with shape factor controlled by the 'Leakage' name-value pair. See periodogram for more details.

• If the resolution bandwidth is RBW_performance, then pspectrum computes a Welch periodogram for the signal. The function:

  1. Divides the signals into overlapping segments.

  2. Windows each segment separately using a Kaiser window with the specified shape factor.

  3. Averages the periodograms of all the segments.

  Welch’s procedure is designed to reduce the variance of the spectrum estimate by averaging different “realizations” of the signals, given by the overlapping sections, and using the window to remove redundant data. See pwelch for more details.

  • The length of each segment (or, equivalently, of the window) is computed using

    $$\text{Segment length}=\frac{f_{\text{Nyquist}}\times \text{ENBW}}{\text{RBW}},$$

    where f_Nyquist is the Nyquist frequency. (If there is no aliasing, the Nyquist frequency is one-half the effective sample rate, defined as the inverse of the median of the differences between adjacent time points. The Nyquist range is [0, f_Nyquist] for real signals and [–f_Nyquist, f_Nyquist] for complex signals.)

  • The stride length is found by adjusting an initial estimate,

    $$\text{Stride length}=\frac{\text{Window length}}{2\times \text{ENBW}-1},$$

    so that the first window starts exactly on the first sample of the first segment and the last window ends exactly on the last sample of the last segment.

Spectrogram Computation

To compute the time-dependent spectrum of a nonstationary signal, pspectrum divides the signal into overlapping segments, windows each segment with a Kaiser window, computes the short-time Fourier transform, and then concatenates the transforms to form a matrix. For more information, see Spectrogram Computation with Signal Processing Toolbox.

A nonstationary signal is a signal whose frequency content changes with time. The spectrogram of a nonstationary signal is an estimate of the time evolution of its frequency content. To construct the spectrogram of a nonstationary signal, pspectrum follows these steps:

1. Divide the signal into equal-length segments. The segments must be short enough that the frequency content of the signal does not change appreciably within a segment. The segments may or may not overlap.

2. Window each segment and compute its spectrum to get the short-time Fourier transform.

3. Use the segment spectra to construct the spectrogram:

   • If called with output arguments, concatenate the spectra to form a matrix.

   • If called with no output arguments, display the power of each spectrum in decibels segment by segment. Depict the magnitudes side-by-side as an image with magnitude-dependent colormap.

The function can compute the spectrogram only for single-channel signals.

Divide Signal into Segments

To construct a spectrogram, first divide the signal into possibly overlapping segments. With the pspectrum function, you can control the length of the segments and the amount of overlap between adjoining segments using the 'TimeResolution' and 'OverlapPercent' name-value pair arguments. If you do not specify the length and overlap, the function chooses a length based on the entire length of the signal and an overlap percentage given by

$$\left(1-\frac{\text{1}}{2\times \text{ENBW}-1}\right)\times 100,$$

where ENBW is the equivalent noise bandwidth of the spectral window. See enbw and Spectrum Computation for more information.

Specified Time Resolution

• If the signal does not have time information, specify the time resolution (segment length) in samples. The time resolution must be an integer greater than or equal to 1 and smaller than or equal to the signal length.

  If the signal has time information, specify the time resolution in seconds. The function converts the result into a number of samples and rounds it to the nearest integer that is less than or equal to the number but not smaller than 1. The time resolution must be smaller than or equal to the signal duration.

• Specify the overlap as a percentage of the segment length. The function converts the result into a number of samples and rounds it to the nearest integer that is less than or equal to the number.

Default Time Resolution

If you do not specify a time resolution, then pspectrum uses the length of the entire signal to choose the length of the segments. The function sets the time resolution as ⌈N/d⌉ samples, where the ⌈⌉ symbols denote the ceiling function, N is the length of the signal, and d is a divisor that depends on N:

Signal Length (N)	Divisor (d)	Segment Length
2 samples – 63 samples	2	1 sample – 32 samples
64 samples – 255 samples	8	8 samples – 32 samples
256 samples – 2047 samples	8	32 samples – 256 samples
2048 samples – 4095 samples	16	128 samples – 256 samples
4096 samples – 8191 samples	32	128 samples – 256 samples
8192 samples – 16383 samples	64	128 samples – 256 samples
16384 samples – N samples	128	128 samples – ⌈N / 128⌉ samples

You can still specify the overlap between adjoining segments. Specifying the overlap changes the number of segments. Segments that extend beyond the signal endpoint are zero-padded.

Consider the seven-sample signal [s0 s1 s2 s3 s4 s5 s6]. Because ⌈7/2⌉ = ⌈3.5⌉ = 4, the function divides the signal into two segments of length four when there is no overlap. The number of segments changes as the overlap increases.

Number of Overlapping Samples	Resulting Segments
0	s0 s1 s2 s3 s4 s5 s6 0
1	s0 s1 s2 s3 s3 s4 s5 s6
2	s0 s1 s2 s3 s2 s3 s4 s5 s4 s5 s6 0
3	s0 s1 s2 s3 s1 s2 s3 s4 s2 s3 s4 s5 s3 s4 s5 s6

pspectrum zero-pads the signal if the last segment extends beyond the signal endpoint. The function returns t, a vector of time instants corresponding to the centers of the segments.

Window the Segments and Compute Spectra

After pspectrum divides the signal into overlapping segments, the function windows each segment with a Kaiser window. The shape factor β of the window, and therefore the leakage, can be adjusted using the 'Leakage' name-value pair. The function then computes the spectrum of each segment and concatenates the spectra to form the spectrogram matrix. To compute the segment spectra, pspectrum follows the procedure described in Spectrum Computation, except that the lower limit of the resolution bandwidth is

$${\text{RBW}}_{\text{performance}}=4\times \frac{f_{\text{span}}}{1024-1}.$$

Display Spectrum Power

If called with no output arguments, the function displays the power of the short-time Fourier transform in decibels, using a color bar with the default MATLAB colormap. The color bar comprises the full power range of the spectrogram.

Persistence Spectrum Computation

The persistence spectrum is a time-frequency view that shows the percentage of the time that a given frequency is present in a signal. The persistence spectrum is a histogram in power-frequency space. The longer a particular frequency persists in a signal as the signal evolves, the higher its time percentage and thus the brighter or "hotter" its color in the display. Use the persistence spectrum to identify signals hidden in other signals.

To compute the persistence spectrum, pspectrum performs these steps:

1. Compute the spectrogram using the specified leakage, time resolution, and overlap. See Spectrogram Computation for more details.

2. Partition the power and frequency values into 2-D bins. (Use the 'NumPowerBins' name-value pair to specify the number of power bins.)

3. For each time value, compute a bivariate histogram of the logarithm of the power spectrum. For every power-frequency bin where there is signal energy at that instant, increase the corresponding matrix element by 1. Sum the histograms for all the time values.

4. Plot the accumulated histogram against the power and the frequency, with the color proportional to the logarithm of the histogram counts expressed as normalized percentages. To represent zero values, use one-half of the smallest possible magnitude.

Power Spectra

Histograms

Accumulated Histogram

References

[1] harris, fredric j. “On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform.” Proceedings of the IEEE^®. Vol. 66, January 1978, pp. 51–83.

[2] Welch, Peter D. "The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms." IEEE Transactions on Audio and Electroacoustics. Vol. 15, June 1967, pp. 70–73.

Extended Capabilities

C/C++ Code Generation
Generate C and C++ code using MATLAB® Coder™.

GPU Arrays

Accelerate code by running on a graphics processing unit (GPU) using Parallel Computing Toolbox™.

The pspectrum function supports GPU array input with these usage notes and limitations:

• Persistence spectrum is not supported.

• Reassigned spectrum or spectrogram is not supported.

For more information, see Run MATLAB Functions on a GPU (Parallel Computing Toolbox).

Version History

Introduced in R2017b

R2023b: Use timetables for code generation

The pspectrum function supports MATLAB timetables for code generation.

R2023a: Visualize function outputs using Create Plot Live Editor task

You can now use the Create Plot Live Editor task to visualize the output of pspectrum interactively. You can select different chart types and set optional parameters. The task also automatically generates code that becomes part of your live script.