View Raw SPL
/*****************************************************************************
* *
* MSCOHERE.SPL Copyright (C) 2021 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: Welch coherence calculation *
* *
* Revisions: 27 Oct 2021 RRR Creation *
* *
*****************************************************************************/
#include
#if @HELP_MSCOHERE
MSCOHERE
Purpose: Estimates the magnitude squared coherence using the Welch method.
Syntax: MSCOHERE(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 coherence. 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 COH.
"onesided" : one sided COH from 0 to Fs/2 (default)
"twosided" : full COH from 0 to Fs
"center" : centered COH from -Fs/2 to Fs/2
Returns: A series or array.
Example:
W1: gsin(10000, 1/10000, 3000)
W2: gnorm(length(w1), deltax(w1)) + w1
W3: mscohere(W1, W2)
W1 contains 10000 samples of a 3 kHz sine sampled at 10 kHz.
W2 adds random noise to the sine.
W3 estimates the magnitude squared coherence function by
dividing W1 and W2 into 8 overlapping segments multiplied
by a Hamming window and averaging the FFTs of the
segments. The spectral peak at 3 kHz is clearly visible.
Example:
W1: gsin(10000, 1/10000, 3000)
W2: gnorm(length(w1), deltax(w1)) + w1
W3: mscohere(W1, W2, hamming(1000), 500, 1024)
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)
W2: cumsum(gnorm(length(w1), deltax(w1))) + w1
W3: mscohere(W1, W2, hamming(1000), 500, 1024)
W4: mscohere(W1, W2, hamming(1000), 500, 1024, "constant")
Similar to above except the output noise is not normally
distributed. W4 computes the coherence in the same manner
as W3 except the mean of each segment is removed before
FFT processing.
Example:
W1: gsin(10000, 1/10000, 3000);setvunits("V")
W2: gnorm(length(w1), deltax(w1)) + w1;setvunits("V")
W3: mscohere(W1, W2, hamming(1000), 500, 2000)
W4: 10*log10(mag(W3))
W5: mscohere(W1, W2, hamming(1000, 4), 500, 2940..3060)
W6: 10*log10(mag(W5));overp(W4)
Similar to above, except W4 displays the log magnitude
squared of the coherence. W5 computes the coherence over the
frequency range from 2940 Hz to 3060 Hz. W6 shows the log
magnitude squared of the coherence in W5 and overplots the
full log magnitude squared coherence in W4.
Remarks:
MSCOHERE estimates the magnitude squared coherence function
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 magnitude
squared FFTs of the segments are averaged.
The magnitude squared coherence function is determined by:
mag(WPXY)^2 / (WPXX * WPYY)
Series X is often a white noise process and Y is the
response of the system to the noise process.
By convention, the complex conjugate of X is used to compute
the coherence.
The amplitudes of the coherence estimate should lie between
0 and 1 inclusive.
If NFFT is an array of one or more target frequencies, MSCOHERE
uses the Goertzel algorithm to compute the coherence 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.
If RANGE is "twosided", the full coherence is displayed from
0 to rate(X) in wrap around order. For example, a 10 point
full coherence 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 coherence is displayed from
-rate(X)/2 to rate(X)/2 with the zero frequency in the
middle. For example, a 10 point full coherence 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}
See Also:
Cpsd
Hamming
Hanning
Kaiser
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 coherence */
ITERATE mscohere(x, y, win, olap, nfft, fs, detrend, zeropad, range, output)
{
local wcoh, wpxx, wpyy, wpxy, hlen, dflag, welch_flag, 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_WCOH + 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, wpxx, wpyy) = welchseg(x, win, y, nfft, olap, length(nfft), fs, welch_flag);
}
else
{
(wpxy, wpxx, wpyy) = welchseg(x, win, y, olap, nfft, fs, welch_flag);
}
/* full coherence */
if (output == "complex")
{
/* complex */
wcoh = wpxy / sqrt(wpxx * wpyy);
}
else if (output == "magnitude")
{
/* magnitude */
wcoh = mag(wpxy) / sqrt(wpxx * wpyy);
}
else
{
/* magnitude squared */
wcoh = mag(wpxy)^2 / (wpxx * wpyy);
}
if (isarray(nfft))
{
wcoh = xy(nfft, wcoh);
}
else if (range == "onesided")
{
/* half length */
hlen = int(nfft / 2) + (nfft % 2 == 0);
/* extract first half */
wcoh = extract(wcoh, 1, hlen);
}
else if (range == "centered")
{
wcoh = welchshift(wcoh);
}
/* x units */
hu = gethunits(psd(x[1..2]));
sethunits(wcoh, hu);
if (outargc > 1)
{
/* coherence and frequency */
return(yvals(wcoh), xvals(wcoh));
}
else
{
return(wcoh);
}
}