# Documentation

### This is machine translation

Translated by
Mouse over text to see original. Click the button below to return to the English verison of the page.

# cpsd

Cross power spectral density

## Syntax

`Pxy = cpsd(x,y)Pxy = cpsd(x,y,window)Pxy = cpsd(x,y,window,noverlap)[Pxy,W] = cpsd(x,y,window,noverlap,nfft)[Pxy,F] = cpsd(x,y,window,noverlap,nfft,fs)[...] = cpsd(...,'twosided')cpsd(...)`

## Description

`Pxy = cpsd(x,y)` estimates the cross power spectral density, `Pxy`, of two discrete-time signals, `x` and `y`, using Welch's averaged, modified periodogram method of spectral estimation.

The input signals may be either vectors or two-dimensional matrices. If both are vectors, they must have the same length. If both are matrices, they must have the same size, and `cpsd` operates columnwise: `Pxy(:,n) = cpsd(x(:,n),y(:,n))`. If one is a matrix and the other is a vector, then the vector is converted to a column vector and internally expanded so both inputs have the same number of columns.

For real `x` and `y`, `cpsd` returns a one-sided CPSD and for complex `x` or `y`, it returns a two-sided CPSD.

`cpsd` uses the following default values:

Parameter

Description

Default Value

`nfft`

FFT length which determines the frequencies at which the PSD is estimated

For real `x` and `y`, the length of `Pxy` is (`nfft`/2+1) if `nfft` is even or (`nfft`+1)/2 if `nfft` is odd. For complex `x` or `y`, the length of `Pxy` is `nfft`.

If `nfft` is greater than the signal length, the data is zero-padded. If `nfft` is less than the signal length, the segment is wrapped so that the length is equal to `nfft`.

Maximum of 256 or the next power of 2 greater than the length of each section of `x` or `y`

`fs`

Sampling frequency

1

`window`

Windowing function and number of samples to use for each section

Periodic Hamming window of length to obtain eight equal sections of `x` and `y`

`noverlap`

Number of samples by which the sections overlap

Value to obtain 50% overlap

 Note   You can use the empty matrix, `[]`, to specify the default value for any input argument except `x` or `y`. For example, `Pxy = cpsd(x,y,[],[],128)` uses a Hamming window, default `noverlap` to obtain 50% overlap, and the specified 128 `nfft`.

`Pxy = cpsd(x,y,window)` specifies a windowing function, divides `x` and `y` into overlapping sections of the specified window length, and windows each section using the specified window function. If you supply a scalar for `window`, `Pxy` uses a Hamming window of that length. `x` and `y` are divided into eight equal sections of that length. If the signal cannot be sectioned evenly with 50% overlap, it is truncated.

`Pxy = cpsd(x,y,window,noverlap)` overlaps the sections of `x` by `noverlap` samples. `noverlap` must be an integer smaller than the length of `window`.

`[Pxy,W] = cpsd(x,y,window,noverlap,nfft)` uses the specified FFT length `nfft` in estimating the CPSD. It also returns `W`, which is the vector of normalized frequencies (in rad/sample) at which the CPSD is estimated. For real signals, the range of `W` is [0, π] when `nfft` is even and [0, π) when `nfft` is odd. For complex signals, the range of `W` is [0, 2π).

`[Pxy,F] = cpsd(x,y,window,noverlap,nfft,fs)` returns `Pxy` as a function of frequency and a vector `F` of frequencies at which the CPSD is estimated. `fs` is the sampling frequency in Hz. For real signals, the range of `F` is [0, `fs`/2] when `nfft` is even and [0, `fs`/2) when `nfft` is odd. For complex signals, the range of `F` is [0, `fs`).

`[...] = cpsd(...,'twosided')` returns the two-sided CPSD of real signals `x` and `y`. The length of the resulting `Pxy` is `nfft` and its range is [0, 2π) if you do not specify `fs`. If you specify `fs`, the range is [0, `fs`). Entering `'onesided'` for a real signal produces the default. You can place the `'onesided'` or `'twosided'` option in any position after the `noverlap` parameter.

`cpsd(...)` plots the CPSD versus frequency in the current figure window.

## Examples

collapse all

Generate two colored noise signals and plot their CPSD. Specify a length-1024 FFT and a 500-point triangular window with no overlap.

```rng default h = fir1(30,0.2,rectwin(31)); h1 = ones(1,10)/sqrt(10); r = randn(16384,1); x = filter(h1,1,r); y = filter(h,1,x); cpsd(x,y,triang(500),250,1024) ```

Generate two 100 Hz sinusoidal signals sampled at 1 kHz for 296 ms. One of the sinusoids lags the other by 2.5 ms, equivalent to a phase lag of π/2. Both signals are embedded in white Gaussian noise of variance 1/42. Reset the random number generator for reproducible results.

```rng('default') Fs = 1000; t = 0:1/Fs:0.296; x = cos(2*pi*t*100)+0.25*randn(size(t)); tau = 1/400; y = cos(2*pi*100*(t-tau))+0.25*randn(size(t)); ```

Compute and plot the magnitude of the cross power spectral density. Use the default settings for `cpsd`. The magnitude peaks at the frequency where there is significant coherence between the signals.

```cpsd(x,y,[],[],[],Fs) ```

Plot the phase of the cross spectrum. The ordinate at the high-coherence frequency corresponds to the phase lag between the sinusoids.

```[Pxy,F] = cpsd(x,y,[],[],[],Fs); plot(F,angle(Pxy)) hold on plot(F,2*pi*100*tau*ones(size(F)),'--') hold off xlabel('Hz') ylabel('\Theta(f)') title('Cross Spectrum Phase') ```

collapse all

### Cross Power Spectral Density

The cross power spectral density is the distribution of power per unit frequency and is defined as

`${P}_{xy}\left(\omega \right)=\sum _{m=-\infty }^{\infty }{R}_{xy}\left(m\right){e}^{-j\omega m}.$`

The cross-correlation sequence is defined as

`${R}_{xy}\left(m\right)=E\left\{{x}_{n+m}{y}_{n}^{\ast }\right\}=E\left\{{x}_{n}{y}_{n-m}^{\ast }\right\},$`

where xn and yn are jointly stationary random processes, $-\infty , and E {· } is the expected value operator.

### Algorithms

`cpsd` uses Welch's averaged periodogram method. See the references listed below.

## References

[1] Rabiner, Lawrence R., and B. Gold. Theory and Application of Digital Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1975, pp. 414–419.

[2] Welch, Peter D. "The Use of the 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. AU-15, June 1967, pp. 70–73.

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