Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adding idft #129

Merged
merged 35 commits into from
Jan 28, 2021
Merged
Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
35 commits
Select commit Hold shift + click to select a range
256824f
Merge pull request #1 from xgcm/master
lanougue Nov 10, 2020
e7d8f2b
Merge pull request #2 from xgcm/master
lanougue Nov 17, 2020
7849305
Merge pull request #3 from xgcm/master
lanougue Dec 4, 2020
bb59a8a
updated idft
Dec 4, 2020
148fb26
update idft + fft + ifft
Dec 10, 2020
0ee8ad5
debug spacing_tol
Dec 11, 2020
0515c57
adding tests and warning
Dec 11, 2020
5bf333a
adding dft parseval tests + some typos
Dec 11, 2020
eeec439
moving test to parseval function
Dec 11, 2020
8b97038
correction idft with real
Dec 13, 2020
0502b21
simplification power spectrum
Dec 13, 2020
792cc2b
simplification spectrum and adding true_flags
Dec 14, 2020
6333630
simplification cross phase
Dec 14, 2020
d9f3d87
debug spectrum with False density
Dec 14, 2020
5ff68a7
debug cross-spectrum test
Dec 14, 2020
1e33e06
debug test cross_phase
Dec 14, 2020
f8a04f0
debug test cross_phase 2
Dec 14, 2020
b49f132
adding cross phase tests
Dec 14, 2020
febecb5
adding cross phase tests 2
Dec 14, 2020
b7684d1
adding test_real_dft_true_phase + scaling parameter/density deprecation
Dec 15, 2020
ac511fc
rm old code
Dec 15, 2020
886850e
typo + unknown sclaing
Dec 15, 2020
c4bc754
debug test_isotropize
Dec 15, 2020
c325c54
flag warning + adding test
Dec 15, 2020
7af46ce
update test density error
Dec 15, 2020
a8275f9
test_coordinates + idft default behaviour
Dec 15, 2020
644c731
restoring default dft true_flags
Dec 15, 2020
186a4aa
correct indentation
Dec 15, 2020
5ca6ae7
adding warning
Dec 15, 2020
5beb1aa
user warning
Dec 15, 2020
fbb9d06
density and scaling flag handling
Dec 15, 2020
9a3caaa
restoring phase warning in cross_spectrum
Dec 15, 2020
dcb40ee
function signature
Dec 15, 2020
41ec321
function signature2 + rm comments
Dec 16, 2020
5cb5cd2
dask debug in test
Dec 16, 2020
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
102 changes: 52 additions & 50 deletions xrft/tests/test_xrft.py
Original file line number Diff line number Diff line change
Expand Up @@ -54,16 +54,16 @@ def time_data(request):
return cftime.num2date(np.arange(0, 10 * 365), units, request.param)


class TestDFTImag(object):
def test_dft_1d(self, test_data_1d):
class TestFFTImag(object):
def test_fft_1d(self, test_data_1d):
"""Test the discrete Fourier transform function on one-dimensional data."""

da = test_data_1d
Nx = len(da)
dx = float(da.x[1] - da.x[0]) if "x" in da.dims else 1

# defaults with no keyword args
ft = xrft.dft(da, detrend="constant")
ft = xrft.fft(da, detrend="constant")
# check that the frequency dimension was created properly
assert ft.dims == ("freq_x",)
# check that the coords are correct
Expand All @@ -81,28 +81,28 @@ def test_dft_1d(self, test_data_1d):
npt.assert_allclose(ft_data_expected, ft.values, atol=1e-14)

# redo without removing mean
ft = xrft.dft(da)
ft = xrft.fft(da)
ft_data_expected = np.fft.fftshift(np.fft.fft(da))
npt.assert_allclose(ft_data_expected, ft.values)

# redo with detrending linear least-square fit
ft = xrft.dft(da, detrend="linear")
ft = xrft.fft(da, detrend="linear")
da_prime = sps.detrend(da.values)
ft_data_expected = np.fft.fftshift(np.fft.fft(da_prime))
npt.assert_allclose(ft_data_expected, ft.values, atol=1e-14)

if "x" in da and not da.chunks:
da["x"].values[-1] *= 2
with pytest.raises(ValueError):
ft = xrft.dft(da)
ft = xrft.fft(da)

def test_dft_1d_time(self, time_data):
def test_fft_1d_time(self, time_data):
"""Test the discrete Fourier transform function on timeseries data."""
time = time_data
Nt = len(time)
da = xr.DataArray(np.random.rand(Nt), coords=[time], dims=["time"])

ft = xrft.dft(da, shift=False)
ft = xrft.fft(da, shift=False)

# check that frequencies are correct
if pd.api.types.is_datetime64_dtype(time):
Expand All @@ -112,16 +112,16 @@ def test_dft_1d_time(self, time_data):
freq_time_expected = np.fft.fftfreq(Nt, dt)
npt.assert_allclose(ft["freq_time"], freq_time_expected)

def test_dft_2d(self):
def test_fft_2d(self):
"""Test the discrete Fourier transform on 2D data"""
N = 16
da = xr.DataArray(
np.random.rand(N, N), dims=["x", "y"], coords={"x": range(N), "y": range(N)}
)
ft = xrft.dft(da, shift=False)
ft = xrft.fft(da, shift=False)
npt.assert_almost_equal(ft.values, np.fft.fftn(da.values))

ft = xrft.dft(da, shift=False, window=True, detrend="constant")
ft = xrft.fft(da, shift=False, window=True, detrend="constant")
dim = da.dims
window = np.hanning(N) * np.hanning(N)[:, np.newaxis]
da_prime = (da - da.mean(dim=dim)).values
Expand All @@ -134,19 +134,19 @@ def test_dft_2d(self):
)
assert (xrft.power_spectrum(da, shift=False, density=True) >= 0.0).all()

def test_dim_dft(self):
def test_dim_fft(self):
N = 16
da = xr.DataArray(
np.random.rand(N, N), dims=["x", "y"], coords={"x": range(N), "y": range(N)}
)
npt.assert_array_equal(
xrft.dft(da, dim="y", shift=False).values,
xrft.dft(da, dim=["y"], shift=False).values,
xrft.fft(da, dim="y", shift=False).values,
xrft.fft(da, dim=["y"], shift=False).values,
)
assert xrft.dft(da, dim="y").dims == ("x", "freq_y")
assert xrft.fft(da, dim="y").dims == ("x", "freq_y")

@pytest.mark.parametrize("dask", [False, True])
def test_dft_3d_dask(self, dask):
def test_fft_3d_dask(self, dask):
"""Test the discrete Fourier transform on 3D dask array data"""
N = 16
da = xr.DataArray(
Expand All @@ -156,40 +156,40 @@ def test_dft_3d_dask(self, dask):
)
if dask:
da = da.chunk({"time": 1})
daft = xrft.dft(da, dim=["x", "y"], shift=False)
daft = xrft.fft(da, dim=["x", "y"], shift=False)
npt.assert_almost_equal(
daft.values, np.fft.fftn(da.chunk({"time": 1}).values, axes=[1, 2])
)
da = da.chunk({"x": 1})
with pytest.raises(ValueError):
xrft.dft(da, dim=["x"])
xrft.fft(da, dim=["x"])
with pytest.raises(ValueError):
xrft.dft(da, dim="x")
xrft.fft(da, dim="x")

da = da.chunk({"time": N})
daft = xrft.dft(da, dim=["time"], shift=False, detrend="linear")
daft = xrft.fft(da, dim=["time"], shift=False, detrend="linear")
da_prime = sps.detrend(da, axis=0)
npt.assert_almost_equal(daft.values, np.fft.fftn(da_prime, axes=[0]))
npt.assert_array_equal(
daft.values, xrft.dft(da, dim="time", shift=False, detrend="linear")
daft.values, xrft.fft(da, dim="time", shift=False, detrend="linear")
)

@pytest.mark.skip(reason="3D detrending not implemented")
def test_dft_4d(self):
def test_fft_4d(self):
"""Test the discrete Fourier transform on 2D data"""
N = 16
da = xr.DataArray(
np.random.rand(N, N, N, N),
dims=["time", "z", "y", "x"],
coords={"time": range(N), "z": range(N), "y": range(N), "x": range(N)},
)
ft = xrft.dft(da, shift=False)
ft = xrft.fft(da, shift=False)
npt.assert_almost_equal(ft.values, np.fft.fftn(da.values))

dim = ["time", "y", "x"]
da_prime = xrft.detrend(da[:, 0], dim) # cubic detrend over time, y, and x
npt.assert_almost_equal(
xrft.dft(
xrft.fft(
da[:, 0].drop("z"),
dim=dim,
shift=False,
Expand All @@ -199,8 +199,8 @@ def test_dft_4d(self):
)


class TestDFTReal(object):
def test_dft_real_1d(self, test_data_1d):
class TestfftReal(object):
def test_fft_real_1d(self, test_data_1d):
"""
Test the discrete Fourier transform function on one-dimensional data.
"""
Expand All @@ -209,7 +209,7 @@ def test_dft_real_1d(self, test_data_1d):
dx = float(da.x[1] - da.x[0]) if "x" in da.dims else 1

# defaults with no keyword args
ft = xrft.dft(da, real="x", detrend="constant")
ft = xrft.fft(da, real="x", detrend="constant")
# check that the frequency dimension was created properly
assert ft.dims == ("freq_x",)
# check that the coords are correct
Expand All @@ -227,9 +227,9 @@ def test_dft_real_1d(self, test_data_1d):
npt.assert_allclose(ft_data_expected, ft.values, atol=1e-14)

with pytest.raises(ValueError):
xrft.dft(da, real="y", detrend="constant")
xrft.fft(da, real="y", detrend="constant")

def test_dft_real_2d(self):
def test_fft_real_2d(self):
"""
Test the real discrete Fourier transform function on one-dimensional
data. Non-trivial because we need to keep only some of the negative
Expand All @@ -244,11 +244,11 @@ def test_dft_real_2d(self):
dx = float(da.x[1] - da.x[0])
dy = float(da.y[1] - da.y[0])

daft = xrft.dft(da, real="x")
daft = xrft.fft(da, real="x")
npt.assert_almost_equal(
daft.values, np.fft.rfftn(da.transpose("y", "x")).transpose()
)
npt.assert_almost_equal(daft.values, xrft.dft(da, dim=["y"], real="x"))
npt.assert_almost_equal(daft.values, xrft.fft(da, dim=["y"], real="x"))

actual_freq_x = daft.coords["freq_x"].values
expected_freq_x = np.fft.rfftfreq(Nx, dx)
Expand All @@ -268,20 +268,20 @@ def test_chunks_to_segments():
)

with pytest.raises(ValueError):
xrft.dft(
xrft.fft(
da.chunk(chunks=((20, N, N), (N - 20, N, N))),
dim=["time"],
detrend="linear",
chunks_to_segments=True,
)

ft = xrft.dft(
ft = xrft.fft(
da.chunk({"time": 16}), dim=["time"], shift=False, chunks_to_segments=True
)
assert ft.dims == ("time_segment", "freq_time", "y", "x")
data = da.chunk({"time": 16}).data.reshape((2, 16, N, N))
npt.assert_almost_equal(ft.values, dsar.fft.fftn(data, axes=[1]), decimal=7)
ft = xrft.dft(
ft = xrft.fft(
da.chunk({"y": 16, "x": 16}),
dim=["y", "x"],
shift=False,
Expand All @@ -306,7 +306,7 @@ def test_chunks_to_segments():
dims=["time", "y", "x"],
coords={"time": range(N), "y": range(N), "x": range(N)},
)
ft2 = xrft.dft(
ft2 = xrft.fft(
da2.chunk({"y": 16, "x": 16}),
dim=["y", "x"],
shift=False,
Expand All @@ -326,11 +326,11 @@ def test_chunks_to_segments():
)


def test_dft_nocoords():
def test_fft_nocoords():
# Julius' example
# https://github.com/rabernat/xrft/issues/17
data = xr.DataArray(np.random.random([20, 30, 100]), dims=["time", "lat", "lon"])
dft = xrft.dft(data, dim=["time"])
dft = xrft.fft(data, dim=["time"])
ps = xrft.power_spectrum(data, dim=["time"])


Expand Down Expand Up @@ -367,19 +367,19 @@ def test_power_spectrum(self, dask):
ps = xrft.power_spectrum(
da, dim=["y", "x"], window=True, density=False, detrend="constant"
)
daft = xrft.dft(da, dim=["y", "x"], detrend="constant", window=True)
daft = xrft.fft(da, dim=["y", "x"], detrend="constant", window=True)
npt.assert_almost_equal(ps.values, np.real(daft * np.conj(daft)))
npt.assert_almost_equal(np.ma.masked_invalid(ps).mask.sum(), 0.0)

ps = xrft.power_spectrum(
da, dim=["y"], real="x", window=True, density=False, detrend="constant"
)
daft = xrft.dft(da, dim=["y"], real="x", detrend="constant", window=True)
daft = xrft.fft(da, dim=["y"], real="x", detrend="constant", window=True)
npt.assert_almost_equal(ps.values, np.real(daft * np.conj(daft)))

### Normalized
ps = xrft.power_spectrum(da, dim=["y", "x"], window=True, detrend="constant")
daft = xrft.dft(da, dim=["y", "x"], window=True, detrend="constant")
daft = xrft.fft(da, dim=["y", "x"], window=True, detrend="constant")
test = np.real(daft * np.conj(daft)) / N ** 4
dk = np.diff(np.fft.fftfreq(N, 1.0))[0]
test /= dk ** 2
Expand All @@ -390,7 +390,7 @@ def test_power_spectrum(self, dask):
ps = xrft.power_spectrum(
da, dim=["y", "x"], window=True, density=False, detrend="linear"
)
daft = xrft.dft(da, dim=["y", "x"], window=True, detrend="linear")
daft = xrft.fft(da, dim=["y", "x"], window=True, detrend="linear")
npt.assert_almost_equal(ps.values, np.real(daft * np.conj(daft)))
npt.assert_almost_equal(np.ma.masked_invalid(ps).mask.sum(), 0.0)

Expand Down Expand Up @@ -421,8 +421,8 @@ def test_cross_spectrum(self, dask):
da = da.chunk({"time": 1})
da2 = da2.chunk({"time": 1})

daft = xrft.dft(da, dim=dim, shift=True, detrend="constant", window=True)
daft2 = xrft.dft(da2, dim=dim, shift=True, detrend="constant", window=True)
daft = xrft.fft(da, dim=dim, shift=True, detrend="constant", window=True)
daft2 = xrft.fft(da2, dim=dim, shift=True, detrend="constant", window=True)
cs = xrft.cross_spectrum(
da, da2, dim=dim, window=True, density=False, detrend="constant"
)
Expand Down Expand Up @@ -891,19 +891,19 @@ def test_spacing_tol(test_data_1d):
da3 = xr.DataArray(np.random.rand(Nx), coords=[x], dims=["x"])

# This shouldn't raise an error
xrft.dft(da3, spacing_tol=1e-1)
xrft.fft(da3, spacing_tol=1e-1)
# But this should
with pytest.raises(ValueError):
xrft.dft(da3, spacing_tol=1e-4)
xrft.fft(da3, spacing_tol=1e-4)


def test_spacing_tol_float_value(test_data_1d):
da = test_data_1d
with pytest.raises(TypeError):
xrft.dft(da, spacing_tol="string")
xrft.fft(da, spacing_tol="string")


@pytest.mark.parametrize("func", ("dft", "power_spectrum"))
@pytest.mark.parametrize("func", ("fft", "power_spectrum"))
@pytest.mark.parametrize("dim", ["time"])
def test_keep_coords(sample_data_3d, func, dim):
"""Test whether xrft keeps multi-dim coords from rasm sample data."""
Expand Down Expand Up @@ -969,7 +969,9 @@ def test_true_phase():
f = np.fft.fftfreq(len(x), dx)
expected = np.fft.fft(np.fft.ifftshift(y)) * np.exp(-1j * 2.0 * np.pi * f * lag)
expected = xr.DataArray(expected, dims="freq_x", coords={"freq_x": f})
output = xrft.dft(s, dim="x", true_phase=True, shift=False, prefix="freq_")
output = xrft.dft(
s, dim="x", true_phase=True, true_amplitude=False, shift=False, prefix="freq_"
)
xrt.assert_allclose(expected, output)


Expand All @@ -982,8 +984,8 @@ def test_theoretical_matching(rtol=1e-8, atol=1e-3):
y = np.cos(2.0 * np.pi * f0 * x)
y[np.abs(x) >= (T / 2.0)] = 0.0
s = xr.DataArray(y, dims=("x",), coords={"x": x})
S = (
xrft.dft(s, dim="x", true_phase=True) * dx
S = xrft.dft(
s, dim="x", true_phase=True, true_amplitude=True
) # Fast Fourier Transform of original signal
f = S.freq_x # Frequency axis
TF_s = xr.DataArray(
Expand Down
Loading