xspectrogram

Cross-spectrogram using short-time Fourier transforms

Syntax

s = xspectrogram(x,y)

s = xspectrogram(x,y,window)

s = xspectrogram(x,y,window,noverlap)

s = xspectrogram(x,y,window,noverlap,nfft)

[s,w,t] = xspectrogram(___)

[s,f,t] = xspectrogram(___,fs)

[s,w,t] = xspectrogram(x,y,window,noverlap,w)

[s,f,t] = xspectrogram(x,y,window,noverlap,f,fs)

[___,c] = xspectrogram(___)

[___] = xspectrogram(___,freqrange)

[___] = xspectrogram(___,Name,Value)

[___] = xspectrogram(___,spectrumtype)

xspectrogram(___)

xspectrogram(___,freqloc)

Description

s = xspectrogram(x,y) returns the cross-spectrogram of the signals specified by x and y. The input signals must be vectors with the same number of elements. Each column of s contains an estimate of the short-term, time localized frequency content common to x and y.

s = xspectrogram(x,y,window) uses window to divide x and y into segments and perform windowing.

s = xspectrogram(x,y,window,noverlap) uses noverlap samples of overlap between adjoining segments.

s = xspectrogram(x,y,window,noverlap,nfft) uses nfft sampling points to calculate the discrete Fourier transform.

[s,w,t] = xspectrogram(___) returns a vector of normalized frequencies, w, and a vector of time instants, t, at which the cross-spectrogram is computed. This syntax can include any combination of input arguments from previous syntaxes.

[s,f,t] = xspectrogram(___,fs) returns a vector of frequencies, f, expressed in terms of fs, the sample rate. fs must be the sixth input to xspectrogram. To input a sample rate and still use the default values of the preceding optional arguments, specify these arguments as empty, [].

[s,w,t] = xspectrogram(x,y,window,noverlap,w) returns the cross-spectrogram at the normalized frequencies specified in w.

[s,f,t] = xspectrogram(x,y,window,noverlap,f,fs) returns the cross-spectrogram at the frequencies specified in f.

[___,c] = xspectrogram(___) also returns a matrix, c, containing an estimate of the time-varying complex cross-spectrum of the input signals. The cross-spectrogram, s, is the magnitude of c.

[___] = xspectrogram(___,freqrange) returns the cross-spectrogram over the frequency range specified by freqrange. Valid options for freqrange are 'onesided', 'twosided', and 'centered'.

[___] = xspectrogram(___,Name,Value) specifies additional options using name-value arguments. Options include the minimum threshold and output time dimension.

[___] = xspectrogram(___,spectrumtype) returns short-term cross power spectral density estimates if spectrumtype is specified as 'psd' and returns short-term cross power spectrum estimates if spectrumtype is specified as 'power'.

xspectrogram(___) with no output arguments plots the cross-spectrogram in the current figure window.

xspectrogram(___,freqloc) specifies the axis on which to plot the frequency. Specify freqloc as either 'xaxis' or 'yaxis'.

Examples

Cross-Spectrogram of Linear Chirps

Generate two linear chirps sampled at 1 MHz for 10 milliseconds.

• The first chirp has an initial frequency of 150 kHz that increases to 350 kHz by the end of the measurement.

• The second chirp has an initial frequency of 200 kHz that increases to 300 kHz by the end of the measurement.

Add white Gaussian noise such that the signal-to-noise ratio is 40 dB.

nSamp = 10000;
Fs = 1000e3;
SNR = 40;
t = (0:nSamp-1)'/Fs;

x1 = chirp(t,150e3,t(end),350e3);
x1 = x1+randn(size(x1))*std(x1)/db2mag(SNR);

x2 = chirp(t,200e3,t(end),300e3);
x2 = x2+randn(size(x2))*std(x2)/db2mag(SNR);

Compute and plot the cross-spectrogram of the two chirps. Divide the signals into 200-sample segments and window each segment with a Hamming window. Specify 80 samples of overlap between adjoining segments and a DFT length of 1024 samples.

xspectrogram(x1,x2,hamming(200),80,1024,Fs,'yaxis')

Modify the second chirp so that the frequency rises from 50 kHz to 350 kHz during the measurement. Use a 500-sample Kaiser window with shape factor $$\beta =5$$ to window the segments. Specify 450 samples of overlap and a DFT length of 256. Compute and plot the cross-spectrogram.

x2 = chirp(t,50e3,t(end),350e3);
x2 = x2+randn(size(x2))*std(x2)/db2mag(SNR);

xspectrogram(x1,x2,kaiser(500,5),450,256,Fs,'yaxis')

In both cases, the function highlights the frequency content that the two signals have in common.

Cross-Spectrogram of Speech Signals

Load a file containing two speech signals sampled at 44,100 Hz.

• The first signal is a recording of a female voice saying "transform function."

• The second signal is a recording of the same female voice saying "reform justice."

Plot the two signals. Offset the second signal vertically so both are visible.

load('voice.mat')

% To hear, type soundsc(transform,fs),pause(2),soundsc(reform,fs)

t = (0:length(reform)-1)/fs;

plot(t,transform,t,reform+0.3)
legend('"Transform function"','"Reform justice"')

Compute the cross-spectrogram of the two signals. Divide the signals into 1000-sample segments and window them with a Hamming window. Specify 800 samples of overlap between adjoining segments. Include only frequencies up to 4 kHz.

nwin = 1000;
nvlp = 800;
fint = 0:4000;

[s,f,t] = xspectrogram(transform,reform,hamming(nwin),nvlp,fint,fs);

mesh(t,f,20*log10(s))
view(2)
axis tight

The cross-spectrogram highlights the time intervals where the signals have more frequency content in common. The syllable "form" is particularly noticeable.

Phase Shift Between Two Quadratic Chirps

Generate two quadratic chirps, each sampled at 1 kHz for 2 seconds. Both chirps have an initial frequency of 100 Hz that increases to 200 Hz midway through the measurement. The second chirp has a phase difference of 23° compared to the first.

fs = 1e3;
t = 0:1/fs:2;

y1 = chirp(t,100,1,200,'quadratic',0);
y2 = chirp(t,100,1,200,'quadratic',23);

Compute the complex cross-spectrogram of the chirps to extract the phase shift between them. Divide the signals into 128-sample segments. Specify 120 samples of overlap between adjoining segments. Window each segment using a Kaiser window with shape factor β = 18 and specify a DFT length of 128 samples. Use the plotting functionality of xspectrogram to display the cross-spectrogram.

[~,f,xt,c] = xspectrogram(y1,y2,kaiser(128,18),120,128,fs);

xspectrogram(y1,y2,kaiser(128,18),120,128,fs,'yaxis')

Extract and display the maximum-energy time-frequency ridge of the cross-spectrogram.

[tfr,~,lridge] = tfridge(c,f);

hold on
plot(xt,tfr,'k','linewidth',2)
hold off

The phase shift is the ratio of imaginary part to real part of the time-varying cross-spectrum along the ridge. Compute the phase shift and express it in degrees. Display its mean value.

pshft = angle(c(lridge))*180/pi;

mean(pshft)

ans = -23.0000

Cross-Spectrogram of Complex Signals

Generate two signals, each sampled at 3 kHz for 1 second. The first signal is a quadratic chirp whose frequency increases from 300 Hz to 1300 Hz during the measurement. The chirp is embedded in white Gaussian noise. The second signal, also embedded in white noise, is a chirp with sinusoidally varying frequency content.

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

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

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

Compute and plot the cross-spectrogram of the two signals. Divide the signals into 256-sample segments with 255 samples of overlap between adjoining segments. Use a Kaiser window with shape factor β = 30 to window the segments. Use the default number of DFT points. Center the cross-spectrogram at zero frequency.

nwin = 256;

xspectrogram(x1,x2,kaiser(nwin,30),nwin-1,[],fs,'centered','yaxis')

Compute the power spectrum instead of the power spectral density. Set to zero the values smaller than –40 dB. Center the plot at the Nyquist frequency.

xspectrogram(x1,x2,kaiser(nwin,30),nwin-1,[],fs, ...
    'power','MinThreshold',-40,'yaxis')
title('Cross-Spectrogram of Quadratic Chirp and Complex Chirp')

The thresholding further highlights the regions of common frequency.

Cross-Spectrogram of Two Sequences

Compute and plot the cross-spectrogram of two sequences.

Specify each sequence to be 4096 samples long.

N = 4096;

To create the first sequence, generate a convex quadratic chirp embedded in white Gaussian noise and bandpass filter it.

• The chirp has an initial normalized frequency of 0.1π that increases to 0.8π by the end of the measurement.

• The 16th-order filter passes normalized frequencies between 0.2π and 0.4π rad/sample and has a stopband attenuation of 60 dB.

rx = chirp(0:N-1,0.1/2,N,0.8/2,'quadratic',[],'convex')'+randn(N,1)/100;
dx = designfilt('bandpassiir','FilterOrder',16, ...
    'StopbandFrequency1',0.2,'StopbandFrequency2',0.4, ...
    'StopbandAttenuation',60);
x = filter(dx,rx);

To create the second sequence, generate a linear chirp embedded in white Gaussian noise and bandstop filter it.

• The chirp has an initial normalized frequency of 0.9π that decreases to 0.1π by the end of the measurement.

• The 16th-order filter stops normalized frequencies between 0.6π and 0.8π rad/sample and has a passband ripple of 1 dB.

ry = chirp(0:N-1,0.9/2,N,0.1/2)'+randn(N,1)/100;
dy = designfilt('bandstopiir','FilterOrder',16, ...
    'PassbandFrequency1',0.6,'PassbandFrequency2',0.8, ...
    'PassbandRipple',1);
y = filter(dy,ry);

Plot the two sequences. Offset the second sequence vertically so both are visible.

plot([x y+2])

Compute and plot the cross-spectrogram of x and y. Use a 512-sample Hamming window. Specify 500 samples of overlap between adjoining segments and 2048 DFT points.

xspectrogram(x,y,hamming(512),500,2048,'yaxis')

Set to zero the cross-spectrogram values smaller than –50 dB.

xspectrogram(x,y,hamming(512),500,2048,'MinThreshold',-50,'yaxis')

The spectrogram shows the frequency regions that are enhanced or suppressed by the filters.

Input Arguments

x, y — Input signals

vectors

Input signals, specified as vectors.

Example: cos(pi/4*(0:159))+randn(1,160) specifies a sinusoid embedded in white Gaussian noise.

Data Types: single | double
Complex Number Support: Yes

window — Window

integer | vector | []

Window, specified as an integer or as a row or column vector. Use window to divide the signals into segments.

• If window is an integer, then xspectrogram divides x and y into segments of length window and windows each segment with a Hamming window of that length.

• If window is a vector, then xspectrogram divides x and y into segments of the same length as the vector and windows each segment using window.

If the input signals cannot be divided exactly into an integer number of segments with noverlap overlapping samples, then they are truncated accordingly.

If you specify window as empty, then xspectrogram uses a Hamming window such that x and y are divided into eight segments with noverlap overlapping samples.

For a list of available windows, see Windows.

Example: hann(N+1) and (1-cos(2*pi*(0:N)'/N))/2 both specify a Hann window of length N + 1.

Data Types: single | double

noverlap — Number of overlapped samples

positive integer | []

Number of overlapped samples, specified as a positive integer.

• If window is scalar, then noverlap must be smaller than window.

• If window is a vector, then noverlap must be smaller than the length of window.

If you specify noverlap as empty, then xspectrogram uses a number that produces 50% overlap between segments. If the segment length is unspecified, the function sets noverlap to ⌊N/4.5⌋, where N is the length of the input signals.

Data Types: double | single

nfft — Number of DFT points

positive integer | []

Number of DFT points, specified as a positive integer scalar. If you specify nfft as empty, then xspectrogram sets the DFT length to max(256,2^p), where p = ⌈log₂ N_w⌉ and

• N_w = window if window is a scalar.

• N_w = length(window) if window is a vector.

Data Types: single | double

w — Normalized frequencies

vector

Normalized frequencies, specified as a vector. w must have at least two elements. Normalized frequencies are in rad/sample.

Example: pi./[2 4]

Data Types: double | single

f — Frequencies

vector

Frequencies, specified as a vector. f must have at least two elements. The units of f are specified by the sample rate, fs.

Data Types: double | single

fs — Sample rate

1 Hz (default) | positive scalar

Sample rate, specified as a positive scalar. The sample rate is the number of samples per unit time. If the unit of time is seconds, then the sample rate is in Hz.

Data Types: double | single

freqrange — Frequency range for cross-spectrum estimate

'onesided' | 'twosided' | 'centered'

Frequency range for the cross-spectrum estimate, specified as 'onesided', 'twosided', or 'centered'. For real-valued signals, the default is 'onesided'. For complex-valued signals, the default is 'twosided', and specifying 'onesided' results in an error.

• 'onesided' — Returns the one-sided cross-spectrogram of a real input signal. If nfft is even, then s has nfft/2 + 1 rows and is computed over the interval [0, π] rad/sample. If nfft is odd, then s has (nfft + 1)/2 rows and the interval is [0, π) rad/sample. If you specify fs, then the intervals are respectively [0, fs/2] cycles/unit time and [0, fs/2) cycles/unit time.

• 'twosided' — Returns the two-sided cross-spectrogram of a real or complex signal. s has nfft rows and is computed over the interval [0, 2π) rad/sample. If you specify fs, then the interval is [0, fs) cycles/unit time.

• 'centered' — Returns the centered two-sided cross-spectrogram for a real or complex signal. s has nfft rows. If nfft is even, then s is computed over the interval (–π, π] rad/sample. If nfft is odd, then s is computed over (–π, π) rad/sample. If you specify fs, then the intervals are respectively (–fs/2, fs/2] cycles/unit time and (–fs/2, fs/2) cycles/unit time.

spectrumtype — Cross power spectrum scaling

'psd' (default) | 'power'

Cross power spectrum scaling, specified as 'psd' or 'power'.

• Omitting spectrumtype, or specifying 'psd', returns the cross power spectral density.

• Specifying 'power' scales each estimate of the cross power spectral density by the resolution bandwidth, which depends on the equivalent noise bandwidth of the window and the segment duration. The result is an estimate of the power at each frequency.

freqloc — Frequency display axis

'xaxis' (default) | 'yaxis'

Frequency display axis, specified as 'xaxis' or 'yaxis'.

• 'xaxis' — Displays frequency on the x-axis and time on the y-axis.

• 'yaxis' — Displays frequency on the y-axis and time on the x-axis.

This argument is ignored if you call xspectrogram with output arguments.

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: xspectrogram(x,100,'OutputTimeDimension','downrows') divides x and y into segments of length 100 and windows each segment with a Hamming window of that length. The output of the spectrogram has time dimension down the rows.

MinThreshold — Threshold

-Inf (default) | real scalar

Threshold, specified as a real scalar expressed in decibels. xspectrogram sets to zero those elements of s such that 10 log₁₀(s) ≤ thresh.

OutputTimeDimension — Output time dimension

acrosscolumns (default) | downrows

Output time dimension, specified as acrosscolumns or downrows. Set this value to downrows, if you want the time dimension of s, ps, fc, and tc down the rows and the frequency dimension along the columns. Set this value to acrosscolumns, if you want the time dimension of s, ps, fc, and tc across the columns and frequency dimension along the rows. This input is ignored if the function is called without output arguments.

Output Arguments

s — Cross-spectrogram

matrix

Cross-spectrogram, returned as a matrix. Time increases across the columns of s and frequency increases down the rows, starting from zero.

• If the input signals x and y are of length N, then s has k columns, where:

  • k = ⌊(N – noverlap)/(window – noverlap)⌋ if window is a scalar.

  • k = ⌊(N – noverlap)/(length(window) – noverlap)⌋ if window is a vector.

• If the input signals are real and nfft is even, then s has (nfft/2 + 1) rows.

• If the input signals are real and nfft is odd, then s has (nfft + 1)/2 rows.

• If the input signals are complex, then s has nfft rows.

Data Types: double | single

w — Normalized frequencies

vector

Normalized frequencies, returned as a vector. w has a length equal to the number of rows of s.

Data Types: double | single

t — Time instants

vector

Time instants, returned as a vector. The time values in t correspond to the midpoint of each segment specified using window.

Data Types: double | single

f — Cyclical frequencies

vector

Cyclical frequencies, returned as a vector. f has a length equal to the number of rows of s.

Data Types: double | single

c — Time-varying complex cross-spectrum

matrix

Time-varying complex cross-spectrum, returned as a matrix. The cross-spectrogram, s, is the magnitude of c.

Data Types: double | single

References

[1] Mitra, Sanjit K. Digital Signal Processing: A Computer-Based Approach. 2nd Ed. New York: McGraw-Hill, 2001.

[2] Oppenheim, Alan V., and Ronald W. Schafer, with John R. Buck. Discrete-Time Signal Processing. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.

Extended Capabilities

C/C++ Code Generation

Generate C and C++ code using MATLAB® Coder™.

Usage notes and limitations:

• Arguments specified using name value pairs must be compile time constants.

• Window must be double precision.

Thread-Based Environment

Run code in the background using MATLAB® backgroundPool or accelerate code with Parallel Computing Toolbox™ ThreadPool.

Usage notes and limitations:

• The syntax with no output arguments is not supported.

• The frequency vector must be uniformly spaced.

For more information, see Run MATLAB Functions in Thread-Based Environment.

GPU Arrays

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

Usage notes and limitations:

• The window length must not be greater than nfft or the length of f.

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

Version History

Introduced in R2017a