From a86450998aaf92f35c85bc77501a9358f1292226 Mon Sep 17 00:00:00 2001 From: Shane Elipot Date: Tue, 8 Aug 2023 10:20:41 -0700 Subject: [PATCH] rotary transform function (#225) * rotary transform function * Update clouddrift/signal.py Co-authored-by: Milan Curcic * Update clouddrift/signal.py Co-authored-by: Milan Curcic * @milancurcic and @philippemiron suggestions --------- Co-authored-by: Milan Curcic --- clouddrift/signal.py | 68 ++++++++++++++++++++++++++++++++++++++++++- tests/signal_tests.py | 35 ++++++++++++++++++++++ 2 files changed, 102 insertions(+), 1 deletion(-) diff --git a/clouddrift/signal.py b/clouddrift/signal.py index f535e109..8bc983ca 100644 --- a/clouddrift/signal.py +++ b/clouddrift/signal.py @@ -13,7 +13,7 @@ def analytic_transform( x: Union[list, np.ndarray, xr.DataArray, pd.Series], boundary: Optional[str] = "mirror", ) -> np.ndarray: - """returns the analytic part of a real-valued signal or of a complex-valued + """Return the analytic part of a real-valued signal or of a complex-valued signal. To obtain the anti-analytic part of a complex-valued signal apply analytic_transform to the conjugate of the input. Analytic_transform removes the mean of the input signals. @@ -86,3 +86,69 @@ def analytic_transform( z = z[int(m0 + 1) - 1 : int(2 * m0 + 1) - 1] return z + + +def rotary_transform( + u: Union[list, np.ndarray, xr.DataArray, pd.Series], + v: Union[list, np.ndarray, xr.DataArray, pd.Series], + boundary: Optional[str] = "mirror", +) -> Tuple[np.ndarray, np.ndarray]: + """Return time-domain rotary components time series (zp,zn) from Cartesian components time series (u,v). + Note that zp and zn are both analytic time series which implies that the complex-valued time series + u+1j*v is recovered by zp+np.conj(zn). The mean of the original complex signal is split evenly between + the two rotary components. + + If up is the analytic transform of u, and vp the analytic transform of v, then the counterclockwise and + clockwise components are defined by zp = 0.5*(up+1j*vp), zp = 0.5*(up-1j*vp). + See as an example Lilly and Olhede (2010), doi: 10.1109/TSP.2009.2031729. + + Parameters + ---------- + u : np.ndarray + Real-valued signal, first Cartesian component (zonal, east-west) + v : np.ndarray + Real-valued signal, second Cartesian component (meridional, north-south) + boundary : str, optional ["mirror", "zeros", "periodic"] optionally specifies the + boundary condition to be imposed at the edges of the time series for the underlying analytic + transform. Default is "mirror". + + Returns + ------- + zp : np.ndarray + Time-domain complex-valued positive (counterclockwise) rotary component. + zn : np.ndarray + Time-domain complex-valued negative (clockwise) rotary component. + + Examples + -------- + + To obtain the rotary components of a real-valued signal: + >>> zp, zn = rotary_transform(u,v) + + Raises + ------ + ValueError + If the input arrays do not have the same shape. + + See Also + -------- + :func:`analytic_transform` + """ + # to implement: input one complex argument instead of two real arguments + # u and v arrays must have the same shape or list length. + if type(u) == list and type(v) == list: + if not len(u) == len(v): + raise ValueError("u and v must have the same length.") + else: + if not u.shape == v.shape: + raise ValueError("u and v must have the same shape.") + + muv = np.mean(u) + 1j * np.mean(v) + + if muv == xr.DataArray: + muv = muv.to_numpy() + + up = analytic_transform(u, boundary=boundary) + vp = analytic_transform(v, boundary=boundary) + + return 0.5 * (up + 1j * vp) + 0.5 * muv, 0.5 * (up - 1j * vp) + 0.5 * np.conj(muv) diff --git a/tests/signal_tests.py b/tests/signal_tests.py index 35936f0e..911c5e4b 100644 --- a/tests/signal_tests.py +++ b/tests/signal_tests.py @@ -1,5 +1,6 @@ from clouddrift.signal import ( analytic_transform, + rotary_transform, ) import numpy as np import unittest @@ -52,6 +53,12 @@ def test_analytic_transform_list(self): self.assertTrue(np.allclose(x - np.mean(x), np.real(z))) +def test_analytic_transform_array(self): + x = np.random.rand(99) + z = analytic_transform(x) + self.assertTrue(np.allclose(x - np.mean(x), np.real(z))) + + def test_analytic_transform_pandas(self): x = pd.Series(data=np.random.rand(99)) z = analytic_transform(x) @@ -62,3 +69,31 @@ def test_analytic_transform_xarray(self): x = xr.DataArray(data=np.random.rand(99)) z = analytic_transform(x) self.assertTrue(np.allclose(x - np.mean(x), np.real(z))) + + +def test_rotary_transform_array(self): + u = np.random.rand(99) + v = np.random.rand(99) + zp, zn = rotary_transform(u, v) + self.assertTrue(np.allclose(u + 1j * v, zp + np.conj(zn))) + + +def test_rotary_transform_list(self): + u = list(np.random.rand(99)) + v = list(np.random.rand(99)) + zp, zn = rotary_transform(u, v) + self.assertTrue(np.allclose(u + 1j * v, zp + np.conj(zn))) + + +def test_rotary_transform_pandas(self): + u = pd.Series(data=np.random.rand(99)) + v = pd.Series(data=np.random.rand(99)) + zp, zn = rotary_transform(u, v) + self.assertTrue(np.allclose(u + 1j * v, zp + np.conj(zn))) + + +def test_rotary_transform_xarray(self): + u = xr.DataArray(data=np.random.rand(99)) + v = xr.DataArray(data=np.random.rand(99)) + zp, zn = rotary_transform(u, v) + self.assertTrue(np.allclose(u + 1j * v, zp + np.conj(zn)))