Question regarding description of cross-spectral density

[Xunius]

Hi,

Thanks for providing the notebook.
I have a question regarding the description of cross-spectral density (CSD) in this notebook:

The cross-spectral density (CSD) estimate is computed using the real valued PSD estimate of $i$ defined as ${\hat{f}}^{MT}{ii}\left(\omega \right)$ and the complex conjugate of the PSD estimate of $j$ defined as ${\hat{f}}^{*MT}{jj}\left(\omega \right)$ and given by:

If the PSD of i is real-valued, then that of j is also real since i and j are symmetrical. Then the CSD is real-valued, too. However, the cij output from the mtem() function is complex. That's why we can compute the phase from the cij.

Could you help clarify this point?

[mattijn]

Hi, it has been years, but the text states that j is complex.

[Xunius]

I still don't get it. The PSD of both i and j should be complex but with imaginary part being 0. So when multiplying one with the conjugate of another we still won't get any imaginary part.

[mattijn]

For me this video Coherence and the Cross Spectrum and Estimation of Coherence and Cross Spectra in combination with this oceanography course material on Test of coherence and phase lag calculations were useful

Maybe you did take the PSD directly on the input signal (then your input signal should be complex), or did you compute the CSD estimate using the PSD estimates?

[Xunius]

Thanks for the materials, I've watched/read them and they are helpful.

However, what I was asking is whether there is an error in the following equation (and the relevant description) from the notebook:

The PSD is always real-valued regardless whether the input signal is real or complex. Then the ${\hat{f}}^{MT}{ii}$ and ${\hat{f}}^{MT}{jj}$ are both real-valued. Then the computed ${\hat{f}}^{MT}{ij}$ term won't have any non-zero imaginary part.

Regarding the Test of coherence and phase lag caculations material, I'm not quite sure about the coherence (or magnitude-squared coherence to be consistent with Estimation of Coherence and Cross Spectra) part. In fact even the author mentioned

% I've commented out the emery/thomson formula, and (incorrectly?) used:
%fftcoh= abs(csdxy).^2 ./ (psdhx.*psdhy);
fftcoh= real(csdxy).^2 ./ (psdhx.*psdhy);

I think it is indeed incorrect because he is taking the square of the real part, rather than the abs. The author mentioned Matlab may have some precision issues but I doubt it is because the directly computed fftcoh will be strictly 1 according to explanation given in the Estimation of Coherence and Cross Spectra video.

[mattijn]

Thanks for the feedback. First of all, I'm open for improvements. Secondly, I'll have to do some reading by myself before I can answer your question with confidence (regardless right or wrong).

[mattijn]

I prepared a bit of code to show how it works. I think fii^MT and fjj^MT should be written as fii^k fjj^k in:

import numpy as np
import matplotlib.pyplot as plt

# prepare two signals
dt = 0.01
t = np.arange(0, 30, dt)
nse1 = np.random.randn(len(t))  # white noise 1
nse2 = np.random.randn(len(t))  # white noise 2
r = np.exp(-t / 0.05)

cnse1 = np.convolve(nse1, r, mode="same") * dt  # colored noise 1
cnse2 = np.convolve(nse2, r, mode="same") * dt  # colored noise 2

# two signals with a coherent part and a random part
s1 = 0.01 * np.sin(2 * np.pi * 10 * t) + cnse1
s2 = 0.01 * np.sin(2 * np.pi * 10 * t) + cnse2

# function for psd-estimates
def psd_est(s, Fs):
    nfft = 256
    num_freqs = nfft // 2 + 1
    freqs = np.fft.fftfreq(nfft, 1 / Fs)[:num_freqs]
    freqs[-1] *= -1

    shape = (nfft, (s.shape[-1]) // nfft)
    strides = (s.strides[0], nfft * s.strides[0])
    psd_ests = np.lib.stride_tricks.as_strided(
        s, shape=shape, strides=strides
    )  # windows of length nfft
    psd_ests = psd_ests - psd_ests.mean(0, keepdims=True)  # detrend
    windows = np.tile(np.hanning(nfft), (shape[1], 1)).T
    psd_ests = psd_ests * windows  # multiply tapers
    psd_ests = np.fft.fft(psd_ests, n=nfft, axis=0)[
        :num_freqs, :
    ]  # fourier on tapered windows

    return psd_ests, freqs

# compute psd estimates (still compex)
psi_ests, f = psd_est(s1, Fs=1 / dt)
psj_ests, f = psd_est(s2, Fs=1 / dt)

# derive csd (still complex)
csd_ests = np.conj(psi_ests) * psj_ests
csd = csd_ests.mean(axis=1)

# get psd's (real-valued)
psi = (np.conj(psi_ests) * psi_ests).mean(axis=1).real
psj = (np.conj(psj_ests) * psj_ests).mean(axis=1).real

# and plot psd and csd side-by-side
fig = plt.figure(figsize=(9, 2))
ax1 = fig.add_subplot(1, 2, 1)
ax1.set_title("Power Spectral Densities")
ax1.plot(f, psi)
ax1.plot(f, psj)

ax2 = fig.add_subplot(1, 2, 2)
ax2.set_title("Cross Spectral Density")
ax2.plot(f, np.real(csd))
ax2.plot(f, np.imag(csd))

[Xunius]

@mattijn Thanks for the update.

I guess my confusion lies in some of the notations. What confused me the most is the double subscripts like \hat{f}_{ii}, \hat{f}_{jj}, \hat{f}_{ij}. Do people use such double subscripts to denote the product of one Fourier coefficient and the conjugate of another? What notation should I use to denote just the Fourier coefficient it self, a single subscript, like \hat{f}_i?

Because according to this equation in the notebook:

fii^k and fjj^k, as you suggested, are the PSDs of i and j, and they are real-valued, consistent with your code.

So I agree that the CSD should be computed as csd_ests = np.conj(psi_ests) * psj_ests. I am just not sure how the math should be formally written.