Read about our upcoming Code of Conduct on this issue

This instance will be upgraded to Heptapod 0.28.1 on 2022-01-26 at 16:00 UTC+1 (a few minutes of down time)

Commit c5d20ade authored by Pierre Augier's avatar Pierre Augier
Browse files

spatiotemporal_spectra: _compute_spectrum

parent b2ff3e87d581
......@@ -530,6 +530,19 @@ class SpatioTemporalSpectra(SpecificOutput):
return series
def _compute_spectrum(self, data):
if not hasattr(self, "f_sample"):
# TODO: Fix this value (can be read from a data file)
self.f_sample = 1
self.domega = 2 * pi * self.f_sample / data.shape[-1]
# TODO: I'm not sure if detrend=False is good in prod, but it's much
# better for testing
freq, spectrum = signal.periodogram(
data, fs=self.f_sample, scaling="spectrum", detrend=False
)
return freq, spectrum / self.domega
def compute_spectra(self, tmin=0, tmax=None, dtype=None):
"""compute spatiotemporal spectra from files"""
if tmax is None:
......@@ -538,10 +551,6 @@ class SpatioTemporalSpectra(SpecificOutput):
# load time series as state_spect arrays + times
series = self.load_time_series(tmin=tmin, tmax=tmax, dtype=dtype)
# get the sampling frequency
times = series["times"]
f_sample = 1 / np.mean(times[1:] - times[:-1])
# compute spectra
print("computing temporal spectra...")
......@@ -551,7 +560,7 @@ class SpatioTemporalSpectra(SpecificOutput):
if "_Fourier" not in key:
continue
key_spectrum = "spectrum_" + key.split("_Fourier")[0]
freq, spectrum = signal.periodogram(data, fs=f_sample)
freq, spectrum = self._compute_spectrum(data)
spectra[key_spectrum] = spectrum
spectra["omegas"] = 2 * pi * freq
......
......@@ -388,18 +388,18 @@ class SpatioTemporalSpectraNS3D(SpatioTemporalSpectra):
del series
# ud
spectra["spectrum_Khd"] = np.zeros(udx_fft.shape, dtype=dtype)
freq, spectrum = signal.periodogram(udx_fft, fs=f_sample)
spectra["spectrum_Khd"] += 0.5 * spectrum
freq, spectrum = signal.periodogram(udy_fft, fs=f_sample)
spectra["spectrum_Khd"] += 0.5 * spectrum
spectra["spectrum_Khd"] = Khd = np.zeros(udx_fft.shape, dtype=dtype)
freq, spectrum = self._compute_spectrum(udx_fft)
Khd += 0.5 * spectrum
freq, spectrum = self._compute_spectrum(udy_fft)
Khd += 0.5 * spectrum
# ur
spectra["spectrum_Khr"] = np.zeros(udx_fft.shape, dtype=dtype)
freq, spectrum = signal.periodogram(urx_fft, fs=f_sample)
spectra["spectrum_Khr"] += 0.5 * spectrum
freq, spectrum = signal.periodogram(ury_fft, fs=f_sample)
spectra["spectrum_Khr"] += 0.5 * spectrum
spectra["spectrum_Khr"] = Khr = np.zeros(udx_fft.shape, dtype=dtype)
freq, spectrum = self._compute_spectrum(urx_fft)
Khr += 0.5 * spectrum
freq, spectrum = self._compute_spectrum(ury_fft)
Khr += 0.5 * spectrum
spectra["omegas"] = 2 * pi * freq
spectra["dims_order"] = order
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment