Правильно преобразовать выход спектра мощности в дисперсию и волновое число - PullRequest
Новая
       128

Правильно преобразовать выход спектра мощности в дисперсию и волновое число

0 голосов
/ 01 мая 2020

Я использую многоточечную функцию из пакета python spectrum для расчета спектра мощности на основе разрезов спутниковых данных о температуре поверхности моря. Мне любопытно, правильно ли я делаю преобразование в волновое число. Я основал свое преобразование вывода из функции pmtm из того, как они это делают в пакете спектра здесь , но эта функция не часто используется для пространственных данных, поэтому мне любопытно, если мой подход, заключающийся в том, чтобы взять бины частоты / волнового числа, разделив их на NFFT, затем умножив их на 1 / dx (фактически пространственное разрешение пикселей в разрезе). Затем я беру только первую половину бинов, потому что вторая половина - просто отрицательный набор, и, наконец, умножаю волновое число каждого бина на отношение NFFT к фактической длине набора данных.

Вывод из это кажется разумным для моих реальных данных, я включил здесь фиктивный набор данных только для тестирования, но я не уверен на 100%, правильно ли я делаю преобразование в волновое число / частоту. 

from spectrum import *
import numpy as np
import matplotlib.pyplot as plt

def pmtm_spec_and_freqs(input_samples, dx_of_pix, NW=4, show=False):
  Sk_complex, weights, eigenvalues = pmtm(input_samples, NW=NW, method='adapt')

  Sk = abs(Sk_complex)**2
  Sk = np.mean(Sk.T * weights, axis=1)

  # calculate frequency bins based on f = x / NFFT * (sampling rate in Hz or 1/dx of pixels))
  NFFT = max(256, 2**nextpow2(len(input_samples)))
  bins = np.linspace(1,NFFT,NFFT)
  f = bins / NFFT * (1/dx_of_pix)
  # then we need to adjust them by the diff of the transect count and NFFT/2
  freq = f[:int(Sk.shape[0]/2)] * 1/(input_samples.shape[0]/int(Sk.shape[0]/2))
  # only take the first half of the values because the second half is the negative set
  real_var = Sk[:int(Sk.shape[0]/2)]

  if show:
    # this is the region of the spectrum that is viable based on min wavenumber = dx*7
    # and max wavenumber is total length / 8
    plt.axvline(1/(len(input_samples)*dx_of_pix/8), ls='--', color='green')
    plt.axvline(1/(dx_of_pix*7), ls='--', color='green')

    plt.plot(freq, real_var)
    plt.xscale('log')
    plt.yscale('log')

  return(real_var, freq)

# run the function with some made up data
data = data_cosine(N=2048, A=0.1, sampling=1024, freq=1/2)

# set the dx to be 2 which means the resolution of the data is 2 units (km)
dx_of_pix=2

pmtm_spec_and_freqs(data, dx_of_pix, NW=2.5, show=True)

function output image

...