cpsd (Signal Processing Toolbox)

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(x,y,...)

Description

Pxy = cpsd(x,y) estimates the cross power spectral density Pxy of the discrete-time signals x and y using the Welch's averaged, modified periodogram method of spectral estimation. The cross power spectral density is the distribution of power per unit frequency.

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 power spectrum 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 using datawrap 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

Parameter	Description	Default Value
`nfft`	FFT length which determines the frequencies at which the power spectrum 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 using datawrap 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. The x and y vectors 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, pi] when nfft is even and [0, pi) when nfft is odd. For complex signals, the range of W is [0, 2*pi).

[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*pi) 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' string in any position after the noverlap parameter.

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

Examples

Generate two colored noise signals and plot their CPSD with a confidence interval of 95%. Specify a length 1024 FFT, a 500 point triangular window with no overlap, and a sampling frequency of 10 Hz:

randn('state',0);
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)

Algorithm

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

See Also

dspdata.psd, mscohere, pburg, pcov, peig, periodogram, pmcov, pmtm, pmusic, pwelch, pyulear, spectrum.welch, tfestimate

References

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

[2] Welch, P.D. "The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms." IEEE Trans. Audio Electroacoust. Vol. AU-15 (June 1967). Pgs. 70-73.

[3] Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing. Upper Saddle River, NJ: Prentice-Hall, 1999, pp. 737.

cplxpair czt