# 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.```

example

````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`.```

example

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

example

````[___,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`.```

example

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

example

````[___] = xspectrogram(___,Name,Value)` specifies additional options using name-value pair 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

collapse all

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.

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.

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 ```

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.

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

collapse all

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, 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`

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`

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,2p), where p = ⌈log2 Nw and

• Nw = `window` if `window` is a scalar.

• Nw = `length(window)` if `window` is a vector.

Data Types: `single` | `double`

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`

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`

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`

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.

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.

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 Pair Arguments

Specify optional comma-separated pairs of `Name,Value` arguments. `Name` is the argument name and `Value` is the corresponding value. `Name` must appear inside quotes. You can specify several name and value pair arguments in any order as `Name1,Value1,...,NameN,ValueN`.

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.

Threshold, specified as the comma-separated pair consisting of `MinThreshold` and a real scalar expressed in decibels. `xspectrogram` sets to zero those elements of `s` such that 10 log10(`s`) ≤ `thresh`.

Output time dimension, specified as the comma-separated pair consisting of `OutputTimeDimension` and `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

collapse all

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`

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

Data Types: `double` | `single`

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`

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

Data Types: `double` | `single`

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., Ronald W. Schafer, and John R. Buck. Discrete-Time Signal Processing. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.