Median-averaging of complex CSDs

[hadrev]

signal.csd does not median-average complex CSDs correctly because of unexpected behavior of np.median for complex-valued arrays (see the numpy issue). Line 589 of signal/spectral.py could be modified so that if the CSD is complex, its real and imaginary parts are separately median-averaged.

Reproducing code example:

import numpy as np
from scipy import signal

s = 0.3 * np.random.uniform(size=100000) # white noise
n1 = np.random.uniform(size=100000)
n2 = np.random.uniform(size=100000)
d1 = s + n1
d2 = s + n2
freq, Pss = signal.csd(s, s, nperseg=100) # mean average PSD of s
_, Pssmed = signal.csd(s, s, nperseg=100, average='median') # median average PSD of s
_, P12 = signal.csd(d1, d2, nperseg=100) # mean average CSD of d1 and d2
_, P12med = signal.csd(d1, d2, nperseg=100, average='median') # median average CSD of d1 and d2

The real and imaginary parts of P12 should have similar variances, but in the median-averaged case the variance of the imaginary portion is much too high.

Scipy/Numpy/Python version information:

import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
1.4.1 1.16.2 sys.version_info(major=3, minor=7, micro=3, releaselevel='final', serial=0)

[larsoner]

@hadrev agreed they should be processed separately (not 100% sure it's correct but seems reasonable to me), would you be up for making a PR for this?

[e-q]

Good catch @hadrev! It definitely looks like that the current behavior of np.median on complex arrays isn't what we want. There isn't a unique notion of a multivariate median, though, which is effectively what we're looking for in this case (See this wikipedia page).

I haven't looked at this carefully yet, but I think that for cross-spectra, we may in fact want something more like a geometric median rather than the medians of the real and imaginary parts of the array, as this may throw away correlations between the two.

[larsoner]

No consensus on this so I'm removing the 1.6.0 milestone