csd (Signal Processing Toolbox)

Estimate the cross spectral density (CSD) of two signals

Syntax

Pxy = csd(x,y)
Pxy = csd(x,y,nfft)
[Pxy,f] = csd(x,y,nfft,fs)
Pxy = csd(x,y,nfft,fs,window)
Pxy = csd(x,y,nfft,fs,window,numoverlap)
Pxy = csd(x,y,...,'dflag')
[Pxy,Pxyc,f] = csd(x,y,nfft,fs,window,numoverlap,p)
csd(x,y,...)

Description

Pxy = csd(x,y) estimates the cross spectral density of the length n sequences x and y using the Welch method of spectral estimation. Pxy = csd(x,y) uses the following default values:

nfft = min(256,length(x))
fs = 2
window is a periodic Hann window of length nfft.
numoverlap = 0

nfft specifies the FFT length that csd uses. This value determines the frequencies at which the cross spectrum is estimated. fs is a scalar that specifies the sampling frequency. window specifies a windowing function and the number of samples csd uses in its sectioning of the x and y vectors. numoverlap is the number of samples by which the sections overlap. Any arguments omitted from the end of the parameter list use the default values shown above.

If x and y are real, csd estimates the cross spectral density at positive frequencies only; in this case, the output Pxy is a column vector of length nfft/2 + 1 for nfft even and (nfft + 1)/2 for nfft odd. If x or y is complex, csd estimates the cross spectral density at both positive and negative frequencies and Pxy has length nfft.

Pxy = csd(x,y,nfft) uses the FFT length nfft in estimating the cross spectral density of x and y. Specify nfft as a power of 2 for fastest execution.

[Pxy,f] = csd(x,y,nfft,fs) returns a vector f of frequencies at which the function evaluates the CSD. f is the same size as Pxy, so plot(f,Pxy) plots the spectrum versus properly scaled frequency. fs has no effect on the output Pxy; it is a frequency scaling multiplier.

Pxy = csd(x,y,nfft,fs,window) specifies a windowing function and the number of samples per section of the x vector. If you supply a scalar for window, csd uses a Hann window of that length. The length of the window must be less than or equal to nfft; csd zero pads the sections if the length of the window is less than nfft. csd returns an error if the length of the window is greater than nfft.

Pxy = csd(x,y,nfft,fs,window,numoverlap) overlaps the sections of x and y by numoverlap samples.

You can use the empty matrix [] to specify the default value for any input argument except x or y. For example,

```
csd(x,y,[],10000)
```

is equivalent to

```
csd(x)
```

but with a sampling frequency of 10,000 Hz instead of the default of 2 Hz.

Pxy = csd(x,y,...,'dflag') specifies a detrend option, where the string 'dflag' is one of the following:

'linear', to remove the best straight-line fit from the prewindowed sections of x and y
'mean', to remove the mean from the prewindowed sections of x and y
'none', for no detrending (default)

The dflag parameter must appear last in the list of input arguments. csd recognizes a dflag string no matter how many intermediate arguments are omitted.

[Pxy,Pxyc,f] = csd(x,y,nfft,fs,window,numoverlap,p) where p is a positive scalar between 0 and 1 returns a vector Pxyc that contains an estimate of the p*100 percent confidence interval for Pxy. Pxyc is a two-column matrix the same length as Pxy. The interval [Pxyc(:,1),Pxyc(:,2)] covers the true CSD with probability p. plot(f,[Pxy Pxyc]) plots the cross spectrum inside the p*100 percent confidence interval. If unspecified, p defaults to 0.95.

csd(x,y,...) plots the CSD versus frequency in the current figure window. If the p parameter is specified, the plot includes the confidence interval.

Examples

Generate two colored noise signals and plot their CSD 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);
csd(x,y,1024,10000,triang(500),0,[])

Algorithm

csd implements the Welch method of spectral density estimation (see references [1] and [2]):

It applies the window specified by the window vector to each successive detrended section.
It transforms each section with an nfft-point FFT.
It forms the periodogram of each section by scaling the product of the transform of the y section and the conjugate of the transformed x section.
It averages the periodograms of the successive overlapping sections to form Pxy, the cross spectral density of x and y.

The number of sections that csd averages is k, where k is

fix((length(x)-numoverlap)/(length(window)-numoverlap))

Diagnostics

An appropriate diagnostic message is displayed when incorrect arguments to csd are used:

Requires window's length to be no greater than the FFT length.
Requires NOVERLAP to be strictly less than the window length.
Requires positive integer values for NFFT and NOVERLAP.
Requires vector (either row or column) input.
Requires inputs X and Y to have the same length.
Requires confidence parameter to be a scalar between 0 and 1.

See Also

cohere, pburg, pmtm, pmusic, pwelch, pyulear, tfe

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.

cremez czt