copower¶
- copower(dx, y1, y2, /, axis=0, **kwargs)[source]¶
Return the co-spectral decomposition and related quantities for two real-valued arrays along an arbitrary axis.
- Parameters
dx (float or array-like) – Dimension step size in physical units. Used to scale
fx
.y1 (array-like) – First input data.
y2 (array-like) – Second input data. Must have same shape as
y1
.axis (int, optional) – Location of the “time” axis.
wintype (str or (str, float), optional) – The window specification, passed to
window
. The resulting weights are used to window the data before carrying out spectral decompositions. See notes for details.nperseg (int, optional) – The time dimension window or segment length, passed to
window
. IfNone
, windowing is not carried out. See notes for details.detrend ({‘constant’, ‘linear’}, optional) – Passed as the
type
argument toscipy.signal.detrend
.'constant'
removes the mean and'linear'
removes the linear trend.
- Returns
fx (array-like) – Frequencies. Units are the inverse of the
dx
units.C (array-like) – Co-power spectrum.
Q (array-like) – Quadrature spectrum.
P1 (array-like) – Power spectrum for the first input data.
P2 (array-like) – Power spectrum for the second input data.
Coh (array-like, optional) – Coherence. Values should range from 0 to 1.
Phi (array-like, optional) – Phase difference in degrees.
Notes
The Fourier coefficients are scaled so that total variance is equal to one half the sum of the right-hand coefficients. This is more natural for the real-valued datasets typically used by physical scientists, and matches the convention from Professor Elizabeth Barnes’s objective analysis course notes. This differs from the numpy convention, which scales the coefficients so that total variance is equal to the sum of squares of all coefficients, or twice the right-hand coefficients.
Windowing is carried out by applying the
wintype
weights to successive time segments of lengthnperseg
(overlapping by one half the window length), taking spectral decompositions of each weighted segment, then taking the average of the result for all segments. Note that non-boxcar windowing reduces the total power amplitude and results in loss of information. It may often be preferable to follow the example of [RH91] and smooth in frequency space with a Gaussian filter after the decomposition has been carried out.The below example shows that the extent of power reduction resulting from non-boxcar windowing depends on the character of the signal.
>>> import numpy as np >>> import climopy as climo >>> state = np.random.RandomState(51423) >>> w = climo.window(200, 'hanning') >>> y1 = np.sin(np.arange(0, 8 * np.pi - 0.01, np.pi / 25)) # basic signal >>> y2 = state.rand(200) # complex signal >>> for y in (y1, y2): ... yvar = y.var() ... Y = (np.abs(np.fft.fft(y)[1:] / y.size) ** 2).sum() ... Yw = (np.abs(np.fft.fft(y * w)[1:] / y.size) ** 2).sum() ... print('Boxcar', Y / yvar) ... print('Hanning', Yw / yvar) Boxcar 1.0 Hanning 0.375 Boxcar 0.9999999999999999 Hanning 0.5728391743988162
References
- RH91
William J. Randel and Isaac M. Held. Phase Speed Spectra of Transient Eddy Fluxes and Critical Layer Absorption. Journal of the Atmospheric Sciences, 48(5):688–697, March 1991. doi:10.1175/1520-0469(1991)048<0688:PSSOTE>2.0.CO;2.
Examples
>>> import numpy as np >>> import xarray as xr >>> import climopy as climo >>> state = np.random.RandomState(51423) >>> ureg = climo.ureg >>> x = xr.DataArray( ... np.arange(1000, dtype=float) * ureg.day, dims=('time',), name='time' ... ) >>> y1 = xr.DataArray( ... state.rand(1000, 50) * ureg.K, dims=('time', 'space'), name='variable' ... ) >>> y2 = xr.DataArray( ... state.rand(1000, 50) * ureg.m / ureg.s, ... dims=('time', 'space'), name='variable', ... ) >>> f, C, Q, P1, P2, Coh, Phi = climo.copower(x, y1, y2, axis=0) >>> f.head() <xarray.DataArray 'time' (f: 5)> <Quantity([0.001 0.002 0.003 0.004 0.005], '1 / day')> Dimensions without coordinates: f Attributes: long_name: frequency >>> C.head() <xarray.DataArray 'variable' (f: 5, space: 5)> <Quantity([[-1.55726091e-04 -2.51099398e-04 -4.51984603e-05 3.38682324e-04 1.02699316e-04] [ 6.03921637e-05 1.84940350e-04 -6.98685563e-06 1.25120193e-04 -5.62898627e-04] [ 2.89041840e-05 1.62381854e-05 5.51445757e-05 -2.75903456e-04 -1.26185549e-05] [ 5.04572726e-05 -6.94459577e-05 2.43580083e-04 -4.08724137e-05 6.86952412e-05] [-2.18000152e-04 1.43110211e-05 -2.26535472e-04 -3.73436882e-05 -1.28860163e-04]], 'kelvin * meter / second')> Coordinates: * f (f) float64 0.001 0.002 0.003 0.004 0.005 Dimensions without coordinates: space >>> Q.head() <xarray.DataArray 'variable' (f: 5, space: 5)> <Quantity([[ 3.39059554e-05 -5.49872491e-05 3.83785647e-05 1.03330854e-04 -9.02664704e-06] [ 1.65124562e-05 -4.40473709e-05 -5.06947002e-07 1.96874623e-04 -1.18780187e-04] [-1.09859742e-04 2.31464562e-05 -4.98439846e-05 5.67827280e-05 3.39652860e-05] [ 1.22685839e-04 -6.57993024e-05 2.23153696e-04 -9.57458276e-06 1.21946552e-06] [-2.82846001e-05 4.62368401e-05 1.54712600e-05 -6.94771696e-05 -2.77645650e-05]], 'kelvin * meter / second')> Coordinates: * f (f) float64 0.001 0.002 0.003 0.004 0.005 Dimensions without coordinates: space