Python signal.detrend方法代码示例

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import detrend [as 别名]
def get_data(filename):
    data = xarray.open_dataset(
        pjoin(CELER_PATH, 'climate/surface', filename), decode_times=False)

    n_times = data[list(data.data_vars.keys())[0]].shape[0]

    X = np.array(data[list(data.data_vars.keys())[0]]).reshape(n_times, -1)

    # remove seasonality
    period = 12
    for m in range(period):
        # TODO using sklearn for preprocessing would be an improvement
        X[m::period] -= np.mean(X[m::period], axis=0)[None, :]
        X[m::period] /= np.std(X[m::period], axis=0)[None, :]
        if np.sum(np.isnan(X[m::period])) > 0:
            X[m::period] = np.where(np.isnan(X[m::period]), 0, X[m::period])

    # remove trend
    X = detrend(X, axis=0, type='linear')

    return X 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import detrend [as 别名]
def test_dft_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)
        npt.assert_almost_equal(ft.values, np.fft.fftn(da.values))

        ft = xrft.dft(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
        npt.assert_almost_equal(ft.values, np.fft.fftn(da_prime*window))

        da = xr.DataArray(np.random.rand(N,N), dims=['x','y'],
                         coords={'x':range(N,0,-1),'y':range(N,0,-1)}
                         )
        assert (xrft.power_spectrum(da, shift=False,
                                   density=True) >= 0.).all() 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import detrend [as 别名]
def test_dft_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)}
                         )
        with pytest.raises(ValueError):
            xrft.dft(da.chunk({'time':8}), dim=['y','x'], detrend='linear')
        ft = xrft.dft(da, shift=False)
        npt.assert_almost_equal(ft.values, np.fft.fftn(da.values))

        da_prime = xrft.detrendn(da[:,0].values, [0,1,2]) # cubic detrend over time, y, and x
        npt.assert_almost_equal(xrft.dft(da[:,0].drop('z'),
                                        dim=['time','y','x'],
                                        shift=False, detrend='linear'
                                        ).values,
                                np.fft.fftn(da_prime)) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import detrend [as 别名]
def test_average_all_segments(self):
        np.random.seed(1234)
        x = np.random.randn(1024)

        fs = 1.0
        window = 'hann'
        nperseg = 16
        noverlap = 8

        # Compare twosided, because onesided welch doubles non-DC terms to
        # account for power at negative frequencies. stft doesn't do this,
        # because it breaks invertibility.
        f, _, Z = stft(x, fs, window, nperseg, noverlap, padded=False,
                       return_onesided=False, boundary=None)
        fw, Pw = welch(x, fs, window, nperseg, noverlap, return_onesided=False,
                       scaling='spectrum', detrend=False)

        assert_allclose(f, fw)
        assert_allclose(np.mean(np.abs(Z)**2, axis=-1), Pw) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import detrend [as 别名]
def test_roundtrip_float32(self):
        np.random.seed(1234)

        settings = [('hann', 1024, 256, 128)]

        for window, N, nperseg, noverlap in settings:
            t = np.arange(N)
            x = 10*np.random.randn(t.size)
            x = x.astype(np.float32)

            _, _, zz = stft(x, nperseg=nperseg, noverlap=noverlap,
                            window=window, detrend=None, padded=False)

            tr, xr = istft(zz, nperseg=nperseg, noverlap=noverlap,
                           window=window)

            msg = '{0}, {1}'.format(window, noverlap)
            assert_allclose(t, t, err_msg=msg)
            assert_allclose(x, xr, err_msg=msg, rtol=1e-4)
            assert_(x.dtype == xr.dtype) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import detrend [as 别名]
def test_roundtrip_padded_signal(self):
        np.random.seed(1234)

        settings = [
                    ('boxcar', 101, 10, 0),
                    ('hann', 1000, 256, 128),
                    ]

        for window, N, nperseg, noverlap in settings:
            t = np.arange(N)
            x = 10*np.random.randn(t.size)

            _, _, zz = stft(x, nperseg=nperseg, noverlap=noverlap,
                            window=window, detrend=None, padded=True)

            tr, xr = istft(zz, noverlap=noverlap, window=window)

            msg = '{0}, {1}'.format(window, noverlap)
            # Account for possible zero-padding at the end
            assert_allclose(t, tr[:t.size], err_msg=msg)
            assert_allclose(x, xr[:x.size], err_msg=msg) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import detrend [as 别名]
def detrend(self, method='linear'):
        """ Remove linear trend from each voxel

        Args:
            type: ('linear','constant', optional) type of detrending

        Returns:
            out: (Brain_Data) detrended Brain_Data instance

        """

        if len(self.shape()) == 1:
            raise ValueError('Make sure there is more than one image in order '
                             'to detrend.')

        out = deepcopy(self)
        out.data = detrend(out.data, type=method, axis=0)
        return out 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import detrend [as 别名]
def fft(self):
        '''Compute the fast Fourier transform using the multitaper method.

        Returns
        -------
        fourier_coefficients : array, shape (n_time_windows, n_trials,
                                             n_tapers, n_fft_samples,
                                             n_signals)

        '''
        time_series = _add_axes(self.time_series)
        time_series = _sliding_window(
            time_series, window_size=self.n_time_samples_per_window,
            step_size=self.n_time_samples_per_step, axis=0)
        if self.detrend_type is not None:
            time_series = detrend(time_series, type=self.detrend_type)

        logger.info(self)

        return _multitaper_fft(
            self.tapers, time_series, self.n_fft_samples,
            self.sampling_frequency).swapaxes(2, -1) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import detrend [as 别名]
def detrend_wrap(detrend_func):
    """
    Wrapper function for `xrft.detrendn`.
    """
    def func(a, axes=None):
        if len(axes) > 3:
            raise ValueError("Detrending is only supported up to "
                            "3 dimensions.")
        if axes is None:
            axes = tuple(range(a.ndim))
        else:
            if len(set(axes)) < len(axes):
                raise ValueError("Duplicate axes are not allowed.")

        for each_axis in axes:
            if len(a.chunks[each_axis]) != 1:
                raise ValueError('The axis along the detrending is upon '
                                'cannot be chunked.')

        if len(axes) == 1:
            return dsar.map_blocks(sps.detrend, a, axis=axes[0],
                                   chunks=a.chunks, dtype=a.dtype
                                  )
        else:
            for each_axis in range(a.ndim):
                if each_axis not in axes:
                    if len(a.chunks[each_axis]) != a.shape[each_axis]:
                        raise ValueError("The axes other than ones to detrend "
                                        "over should have a chunk length of 1.")
            return dsar.map_blocks(detrend_func, a, axes,
                                   chunks=a.chunks, dtype=a.dtype
                                  )

    return func 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import detrend [as 别名]
def _apply_detrend(da, axis_num):
    """Wrapper function for applying detrending"""
    if da.chunks:
        func = detrend_wrap(detrendn)
        da = xr.DataArray(func(da.data, axes=axis_num),
                         dims=da.dims, coords=da.coords)
    else:
        if da.ndim == 1:
            da = xr.DataArray(sps.detrend(da),
                             dims=da.dims, coords=da.coords)
        else:
            da = detrendn(da, axes=axis_num)

    return da 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import detrend [as 别名]
def test_dft_3d_dask(self, dask):
        """Test the discrete Fourier transform on 3D dask array data"""
        N=16
        da = xr.DataArray(np.random.rand(N,N,N), dims=['time','x','y'],
                          coords={'time':range(N),'x':range(N),
                                  'y':range(N)}
                         )
        if dask:
            da = da.chunk({'time': 1})
            daft = xrft.dft(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'])
            with pytest.raises(ValueError):
                xrft.dft(da, dim='x')

            da = da.chunk({'time':N})
            daft = xrft.dft(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')
                                   ) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import detrend [as 别名]
def test_dft_real_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, real='x', detrend='constant')
        # check that the frequency dimension was created properly
        assert ft.dims == ('freq_x',)
        # check that the coords are correct
        freq_x_expected = np.fft.rfftfreq(Nx, dx)
        npt.assert_allclose(ft['freq_x'], freq_x_expected)
        # check that a spacing variable was created
        assert ft['freq_x_spacing'] == freq_x_expected[1] - freq_x_expected[0]
        # make sure the function is lazy
        assert isinstance(ft.data, type(da.data))
        # check that the Fourier transform itself is correct
        data = (da - da.mean()).values
        ft_data_expected = np.fft.rfft(data)
        # because the zero frequency component is zero, there is a numerical
        # precision issue. Fixed by setting atol
        npt.assert_allclose(ft_data_expected, ft.values, atol=1e-14)

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

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import detrend [as 别名]
def test_cross_spectrum(self, dask):
        """Test the cross spectrum function"""
        N = 16
        dim = ['x','y']
        da = xr.DataArray(np.random.rand(2,N,N), dims=['time','x','y'],
                          coords={'time':np.array(['2019-04-18', '2019-04-19'],
                                                 dtype='datetime64'),
                                 'x':range(N), 'y':range(N)})
        da2 = xr.DataArray(np.random.rand(2,N,N), dims=['time','x','y'],
                          coords={'time':np.array(['2019-04-18', '2019-04-19'],
                                                 dtype='datetime64'),
                                  'x':range(N), 'y':range(N)})
        if 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)
        cs = xrft.cross_spectrum(da, da2, dim=dim, window=True, density=False,
                                detrend='constant')
        npt.assert_almost_equal(cs.values, np.real(daft*np.conj(daft2)))
        npt.assert_almost_equal(np.ma.masked_invalid(cs).mask.sum(), 0.)

        cs = xrft.cross_spectrum(da, da2, dim=dim, shift=True, window=True,
                                detrend='constant')
        test = (daft * np.conj(daft2)).real.values/N**4

        dk = np.diff(np.fft.fftfreq(N, 1.))[0]
        test /= dk**2
        npt.assert_almost_equal(cs.values, test)
        npt.assert_almost_equal(np.ma.masked_invalid(cs).mask.sum(), 0.) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import detrend [as 别名]
def test_spectrum_dim(self):
        N = 16
        da = xr.DataArray(np.random.rand(2,N,N), dims=['time','y','x'],
                         coords={'time':np.array(['2019-04-18', '2019-04-19'],
                                                dtype='datetime64'),
                                'y':range(N),'x':range(N)}
                         )

        ps = xrft.power_spectrum(da, dim='y', real='x', window=True,
                                density=False, detrend='constant')
        npt.assert_array_equal(ps.values,
                               xrft.power_spectrum(da, dim=['y'],
                                                  real='x', window=True,
                                                  density=False,
                                                  detrend='constant').values)

        da2 = xr.DataArray(np.random.rand(2,N,N), dims=['time','y','x'],
                          coords={'time':np.array(['2019-04-18', '2019-04-19'],
                                                 dtype='datetime64'),
                                  'y':range(N), 'x':range(N)})
        cs = xrft.cross_spectrum(da, da2, dim='y',
                                shift=True, window=True,
                                detrend='constant')
        npt.assert_array_equal(xrft.cross_spectrum(da, da2, dim=['y'],
                                                  shift=True, window=True,
                                                  detrend='constant').values,
                              cs.values
                              )
        assert ps.dims == ('time','freq_y','freq_x')
        assert cs.dims == ('time','freq_y','x') 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import detrend [as 别名]
def test_isotropic_ps_slope(N=512, dL=1., amp=1e1, s=-3.):
    """Test the spectral slope of isotropic power spectrum."""

    theta = xr.DataArray(synthetic_field(N, dL, amp, s),
                        dims=['y', 'x'],
                        coords={'y':range(N), 'x':range(N)})
    iso_ps = xrft.isotropic_power_spectrum(theta, detrend='constant',
                                         density=True)
    npt.assert_almost_equal(np.ma.masked_invalid(iso_ps[1:]).mask.sum(), 0.)
    y_fit, a, b = xrft.fit_loglog(iso_ps.freq_r.values[4:],
                                 iso_ps.values[4:])

    npt.assert_allclose(a, s, atol=.1)