cpsd() - Signal Processing

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:ParameterDescriptionDefault Value
nfftFFT length which determines the frequencies at which
the PSD is estimatedFor 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
fsSampling frequency1
windowWindowing function and number of samples to use for each
sectionPeriodic Hamming window of length to obtain eight equal
sections of x and y
noverlapNumber of samples by which the sections overlapValue 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' string
in any position after the noverlap parameter.cpsd(...) plots the CPSD
versus frequency in the current figure window.

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(...)

Example

CPSD of Colored Noise SignalsOpen This Example
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)

Cross Spectrum Phase of Lagged SinusoidsOpen This ExampleGenerate 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')

Output / Return Value

Limitations

Alternatives / See Also

Reference