CPSD

Purpose:

Estimates the cross-power spectral density using the Welch method.

Syntax:

CPSD(x, y, win, overlap, nfft, fs, detrend, zeropad, range, output)

Returns:

A complex series or array, the complex cross-power spectral density estimate.

Example:

W1: gsin(10000, 1/10000, 3000);setvunits("V")

W2: gnorm(length(w1), deltax(w1)) + w1;setvunits("V")

W3: cpsd(W1, W2)

W4: 10*log10(mag(W3))

W1 contains 10000 samples of a 3 kHz sine sampled at 10 kHz.

W2 adds random noise to the sine.

W3 estimates the cross-power spectral density function by dividing W1 and W2 into 8 overlapping segments multiplied by a Hamming window and averaging the FFTs of the segments.

W4 displays the dB log magnitude of the complex cross-power spectral density. The spectral peak at 3 kHz is clearly visible.

Example:

W1: gsin(10000, 1/10000, 3000);setvunits("V")

W2: gnorm(length(w1), deltax(w1)) + w1;setvunits("V")

W3: cpsd(W1, W2, hamming(1000), 500, 1024)

W4: 10*log10(mag(W3))

Same as above, except the series are divided into segments of 1000 points overlapped by 500 points and the 1024 point FFT is computed for each segment.

Example:

W1: gsin(10000, 1/10000, 3000);setvunits("V")

W2: gnorm(length(w1), deltax(w1)) + w1;setvunits("V")

W3: cpsd(W1, W2, hamming(1000), 500, 1024, "constant")

W4: 10*log10(mag(W3))

Same as above, except the mean value of each segment is subtracted from the segment before FFT processing.

Example:

W1: gsin(10000, 1/10000, 3000);setvunits("V")

W2: gnorm(length(w1), deltax(w1)) + w1;setvunits("V")

W3: cpsd(W1, W2, hamming(1000), 500, 2000)

W4: 10*log10(mag(W3))

W5: cpsd(W1, W2, hamming(1000, 4), 500, 2940..3060)

W6: 10*log10(mag(W5));overp(W4)

Similar to above, except W4 displays the log magnitude of the CPSD. W5 computes the CPSD over the frequency range from 2940 Hz to 3060 Hz. W6 shows the log magnitude of the CPSD contained in W5 and overplots the full log magnitude CPSD contained in W4.

Example:

W1: gsin(10000, 1/10000, 3000);setvunits("V")

W2: gnorm(length(w1), deltax(w1)) + w1;setvunits("V")

W3: 10*log10(cpsd(W1, W2, hamming(1000), 500, 2000, "mag"))

W4: 10*log10(cpsd(W1, W2, hamming(1000, 4), 500, 2940..3060, "mag"));overp(w3)

Similar to above, except the magnitudes of the cross-power spectral densities are computed directly.

Example:

W1: gcos(10240, 1/1024, 256, pi/2)

W2: gcos(10240, 1/1024, 256, pi/3)

W3: cpsd(w1, w2, 128, 64, 1024)

W4: phase(w3)

p256 = w4[xtoidx(w4, 256)];

W1 and W2 contain 256 Hz cosine series sampled at 1024 Hz. The cosine in W1 has a phase shift of π/2 radians and the cosine in W2 has a phase shift of π/3 radians. The phase difference between the two series is π/3 - π/2 = -π/6 radians.

W3 computes the complex cross-power spectral density by dividing the series into segments of 128 samples overlapped by 64 samples. The 1024 FFTs of each segment are averaged.

W4 computes the phase of the complex cross-power spectral density.

The variable p256 returns the phase of the CPSD at 256 Hz with a value of -0.523599, essentially -π/6.

Remarks:

The cross-power spectral density, S_xy( f ), is the discrete-time Fourier transform (DTFT) of the cross-correlation R_xy( k ) of two series x[n] and y[n]. The cross-correlation for a random process is defined as:

where E is the expected value operator, x[n] and y[n] are a stationary random processes and * indicates complex conjugate. In practice, the cross-correlation is estimated because only a finite sample of an infinite duration random process is available. The estimate of the cross-correlation function for series of length N is defined as:

The digital CPSD is the DTFT of the cross-correlation estimate:

The CPSD can be computed from the DTFTs of the input series X( f ) and Y( f ):

CPSD estimates the complex cross-power spectral density between an input series x and a corresponding output series y using the Welch method (see PWELCH). The input and output series are divided into overlapping segments and the FFTs of the windowed segments are averaged.

By convention, the complex conjugate of x is used to compute the CPSD.

If nfft is an array of one or more target frequencies, CPSD uses the GOERTZEL algorithm to compute the CPSD at each input frequency. The input frequencies can be any real values. Frequency values less than zero or greater than rate(x) are wrapped to the range between 0 and rate(x).

If detrend is "constant" or "linear", the mean or linear trend is computed and removed from each segment before FFT processing.

If zeropad is "zeropad", the last column is padded with zeros to make it the same length as the other columns. If zeropad is "none" or empty, the last column is ignored if it is shorter than the others.

If range is "twosided", the full CPSD is displayed from 0 to rate(x) in wrap-around order. For example, a 10 point full CPSD with values at frequencies:

{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}

is mapped to frequencies:

{-4, -3, -2, -1, 0, 1, 2, 3, 4, 5}

If range is "center", the full CPSD is displayed from -rate(x)/2 to rate(x)/2 with the zero frequency in the middle. For example, a 10 point full CPSD with values at frequencies:

{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}

is mapped to frequencies:

{0, 1, 2, 3, 4, 5, -4, -3, -2, -1}

If both inputs are real, range defaults to "onesided", else if any input is complex, range defaults to "twosided".

By default, the output of CPSD is complex. Use MAGNITUDE and PHASE to compute the magnitude and phase components from the complex CPSD. For example, the magnitude-squared of the CPSD is computed with by:

mag(CPSD(x, y))^2

and phase of the CPSD is computed with:

phase(CPSD(x, y))

By default, the complex CPSD is returned. However, if output is "mag2" or "mag", CPSD returns the magnitude-squared or the magnitude of the complex spectral density directly.

See MSCOHERE to estimate the magnitude-squared coherence function.

See TFESTIMATE to estimate the complex transfer function.

See PWELCH to estimate the power spectral density.

See DADiSP/FFTXL to optimize the underlying FFT computation.

See Also:

FFT

GOERTZELDFT

HAMMING

HANNING

KAISER

MSCOHERE

PSD

PWELCH

SPECTRUM

TFESTIMATE

References: