Расчет квад-спектров из scipy.signal.csd - PullRequest
0 голосов
/ 17 марта 2020

Я пытаюсь рассчитать ко- и квадраспектры из трех 30-минутных наборов данных временных рядов смещений поверхности моря (полученных из буя Datawell Waverider) в вертикальной, северной и западной плоскостях.

Я столкнулся с функцией scipy.signal.csd (которая, кажется, предлагает простой способ сделать это). Эта функция особенно привлекательна тем, что она использует метод Уэлча для оценки поперечной спектральной плотности мощности.

Несмотря на то, что она хороша для обеспечения ко-спектров, я, кажется, не могу получить квад-спектры из него.

Может ли кто-нибудь сообщить, действительно ли функция предоставляет средства для вычисления квад-спектров, и если да, то как можно go узнать об этом?

def Do_cross_power_spectral_density(x,y):
    """
    Estimate the cross power spectral density, Pxy, using Welch’s method.
    See https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.csd.html
    """
    fs = 2.56
    wk = np.hanning(512)
    wnorm = np.sqrt(fs*sum(wk**2))
    window = wk/wnorm

    return(signal.csd(x, y, fs=2.56, window=window, nperseg=512, noverlap=256,\
                      return_onesided=False))    # Do_cross_power_spectral_density()

# Get Cross Periodograms
fhh,Phh = Do_cross_power_spectral_density(heave, heave)
fnn,Pnn = Do_cross_power_spectral_density(north, north)
fww,Pww = Do_cross_power_spectral_density(west, west)
fhn,Phn = Do_cross_power_spectral_density(heave, north)

fhw,Phw = Do_cross_power_spectral_density(heave, west)
fwh,Pwh = Do_cross_power_spectral_density(heave, west)

fnw,Pnw = Do_cross_power_spectral_density(north, west)
fwn,Pwn = Do_cross_power_spectral_density(north, west)

# Determine Co- and Quad- spectra
Chh = Phh[fhh>=0]*2
Cnn = Pnn[fnn>=0]*2
Cww = Pww[fww>=0]*2
Cnw = Pnw[fnw>=0].real*2
Qhn = Phn[fhn>=0].imag
Qhw = Phw[fhw>=0].imag
...