View Raw SPL
/*****************************************************************************
* *
* CPSD.SPL Copyright (C) 2024 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: Welch cross power spectral density *
* *
* Revisions: 25 Feb 2024 RRR Creation *
* *
*****************************************************************************/
#include
#if @HELP_CPSD
CPSD
Purpose: Estimates the cross power spectral density using the Welch method.
Syntax: CPSD(x, y, win, olag, nfft, fs, detrend, zeropad, range)
x - A series or array, the input.
y - A series or array, the response.
win - An integer or series, the segment size or window.
If a series, each input segment is multiplied by
WIN with the segment size set to LENGTH(win).
Defaults to a Hamming window of where the length
is 1/8 the length of the input series.
olap - Optional. An integer, the number of points to
overlap. Defaults to length(WIN)/2.
nfft - Optional. An integer or series. If NFFT is an
integer, NFFT is the FFT length. If NFFT is a
series, the series contains the frequencies used
to evaluate the CPSD. Defaults to the integer
bestpow2(WIN).
fs - Optional. A real, the input series sample rate.
Defaults to RATE(s).
detrend - Optional. A string, the segment detrend mode:
"none" : no detrending (default)
"constant" : remove the mean of the segment
"linear" : remove the linear trend of the
segment
zeropad - Optional. A string, the short segment zero pad mode:
"none" : no zero padding, ignore the last column
if it is shorter than the other columns
(default)
"zeropad" : add zeros to the last column if it is
shorter than the other columns to force
the same length for all columns
range - Optional. A string, the frequency range of the CPSD.
"onesided" : one sided CPSD from 0 to Fs/2 (default)
"twosided" : full CPSD from 0 to Fs
"center" : centered CPSD from -Fs/2 to Fs/2
Returns: A complex series or array, the 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 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 in W5 and overplots the full log magnitude CPSD
in W4.
Remarks:
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.
The output of CPSD is complex.
By convention, the complex conjugate of X is used to compute
the CPSD.
The magnitude squared of the cross power spectral density
is determined by:
mag(CPSD(x, y)^2)
The phase of the cross power spectral density is determined by:
phase(CPSD(x, y))
If NFFT is an array of one or more target frequencies, CPSD
uses the Goertzel algorithm to compute the CPSD magnitudes
at each input frequency. The input frequencies can be any real
values. Frequency values less than zero or greater than rate(S)
are wrapped to the range between 0 and rate(S).
If DETREND is "constant" or "linear", the mean or linear
trend is computed and removed from each segment before FFT
processing.
By default, the RANGE value is "twosided", producing the full
CPSD 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:
{0, 1, 2, 3, 4, 5, -4, -3, -2, -1}
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:
{-4, -3, -2, -1, 0, 1, 2, 3, 4, 5}
If RANGE is "onesided", the CPSD is displayed from
0 to rate(X)/2.
See Also:
Hamming
Hanning
Kaiser
Mscohere
Psd
Pwelch
Tfestimate
References:
[1] Oppenheim & Shafer
Digital Signal Processing
Prentice Hall, 1975
pp. 548-556
[2] Oppenheim & Shafer
Discrete-Time Signal Processing
Prentice Hall, 1989
pp. 730-742
[3] Programs for Digital Signal Processing
IEEE Press, 1979
Section 2.1-1 - 2.1-10
[4] Press, Flannery, Teukolsky, Vetterling
Numerical Recipes in C
Cambridge Press, 1988
pp. 407-552
[5] S. Lawrence Marple, Jr.
Digital Spectral Analysis with Applications
Prentice Hall, 1987
pp. 152-158
#endif
/* welch cross power spectrum */
ITERATE cpsd(x, y, win, olap, nfft, fs, detrend, zeropad, range, output = "complex")
{
local wpxy, hlen, welch_flag, vu, hu;
if (argc < 2 || not(isarray(x)) || not(isarray(y)))
{
error(sprintf("%s - two input series required", __FUNC__));
}
/* parse arguments */
(win, olap, nfft, fs, dflag, range, output) = welch_parse_args(x, y, win, olap, nfft, fs, detrend, zeropad, range, output, "", __FUNC__);
/* detrend option returned as dflag */
welch_flag = (WELCH_F_CPSD + WELCH_F_FULL + WELCH_F_FFTS) + dflag;
if (isarray(nfft))
{
/* explicit frequencies, use Goertzel algorithm */
welch_flag += WELCH_F_GOERTZEL;
/* compute psd of each segment and average */
wpxy = welchseg(x, win, y, nfft, olap, length(nfft), fs, welch_flag);
/* xy */
wpxy = xy(nfft, wpxy);
}
else
{
/* compute FFTs of overlapping segments */
wpxy = welchseg(x, win, y, olap, nfft, fs, welch_flag);
}
if (range == "onesided" && not(isarray(nfft)))
{
/* normalization */
wpxy *= 2.0 / fs / mean(win * win) / length(win);
/* correct DC component */
wpxy[1] /= 2.0;
/* half length */
hlen = int(nfft / 2) + (nfft % 2 == 0);
/* extract first half */
wpxy = extract(wpxy, 1, hlen);
/* correct Nyquist component */
wpxy[hlen] /= 2.0;
}
else
{
/* normalization */
wpxy *= 1.0 / fs / mean(win * win) / length(win);
if (range == "centered")
{
wpxy = fftshift(wpxy);
}
}
/* units */
(vu, hu) = cpsd_units(x, y);
setvunits(wpxy, vu);
sethunits(wpxy, hu);
if (output == "magnitude")
{
wpxy = mag(wpxy);
}
else if (output == "mag2")
{
wpxy = mag(wpxy)^2;
}
if (outargc > 1)
{
/* cross PSD and frequency */
return(yvals(wpxy), xvals(wpxy));
}
else
{
return(wpxy);
}
}
/* units for cross psd */
cpsd_units(x, y)
{
local vu, hu, p;
p = psd(1..2);
setdeltax(p, deltax(x));
vu = getvunits(p * x[1..2] * y[1..2]);
hu = gethunits(psd(x[1..2]));
return(vu, hu);
}