MSCOHERE

Purpose:

Estimates the magnitude-squared coherence using the Welch method.

Syntax:

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

Returns:

A real series or array, the magnitude-squared of the coherence estimate of the input series.

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.

Example:

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

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

W3: mscohere(w1, w2, 128, 64, 1024, "complex")

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 coherence 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 coherence.

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

Remarks:

The complex coherence function C_xy( f ) between two series x[n] and y[n] is defined as:

where:

S_xy( f ) is the cross-power spectral density (CPSD) of x[n] and y[n] (can be complex).

S_xx( f ) and S_yy( f ) are the auto-power spectral densities of x[n] and y[n] (always real and non-negative).

The complex coherence is purely real only if the cross-power spectral density S_xy( f ) is purely real, i.e. the phase of the CPSD is 0 or π meaning x[n] and y[n] are perfectly in-phase or anti-phase.

The magnitude-squared coherence is always real and always lies between 0 and 1 inclusive.

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.

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 magnitude-squared coherence estimate provide a measure of the linear correlation between x and y at frequency f and should lie between 0 (no correlation) and 1 (perfect correlation) 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(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 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:

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

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:

{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 magnitude-squared of the complex coherence estimate is returned. However, if output is "complex", the complex coherence estimate is returned. Use MAGNITUDE and PHASE to compute the magnitude and phase components from the complex coherence. For example, the magnitude of the complex coherence function is computed with:

mag(MSCOHERE(x, y, "complex"))^2

and phase of the complex coherence function is computed with:

phase(MSCOHERE(x, y, "complex"))

See CPSD to estimate the complex cross-power spectral density.

See TFESTIMATE to estimate the complex transfer function.

See PWELCH to estimate the auto-power spectral density.

See DADiSP/FFTXL to optimize the underlying FFT computation.

See Also:

CPSD

HAMMING

HANNING

KAISER

FFT

GOERTZELDFT

MAGNITUDE

PHASE

PHASEDIFF

PSD

PWELCH

SPECTRUM

TFESTIMATE