Python signal.detrend函数代码示例

def slidingWindow(P,inX=3,outX=32,inY=3,outY=64,maxM=50,norm=True):
	""" Enhance the constrast

		Cut off extreme values and demean the image
		Utilize scipy convolve2d to get the mean at a given pixel
		Remove local mean with inner exclusion region

		Args:
			P: 2-d numpy array image
			inX: inner exclusion region in the x-dimension
			outX: length of the window in the x-dimension
			inY: inner exclusion region in the y-dimension
			outY: length of the window in the y-dimension
			maxM: size of the output image in the y-dimension
			norm: boolean to cut off extreme values

		Returns:
			Q: 2-d numpy contrast enhanced
	"""
	Q = P.copy()
	m, n = Q.shape
	
	Q = exposure.equalize_hist(Q.astype('Float32'), nbins = 65)
	Q = detrend(Q.astype('Float32'), axis = 1)
	Q = detrend(Q.astype('Float32'), axis = 0)
	Q = wiener(Q.astype('Float32'), 4)

	return Q[:maxM,:]

def remt(st):
    '''
    rmean and retrend
    '''
    st.data = detrend(st.data, type='linear')
    st.data = detrend(st.data, type='constant')
    return st

 def generateCorrelations(self,doDetrend=True):
     # auto correlation corefficient of u
     if doDetrend:
         ux=signal.detrend(self.ux());
         uy=signal.detrend(self.uy());
         uz=signal.detrend(self.uz());
         umag=signal.detrend(self.Umag());
     else:
         ux=self.ux();
         uy=self.uy();
         uz=self.uz();
         umag=self.Umag();
     #ux=ux[-samples:-1]
     #uy=uy[-samples:-1]
     #uz=uz[-samples:-1]
     self.data['r11'],self.data['taur11'] = tt.xcorr_fft(ux, maxlags=None, norm='coeff')
     self.data['r22'],self.data['taur22'] = tt.xcorr_fft(uy, maxlags=None, norm='coeff')
     self.data['r33'],self.data['taur33'] = tt.xcorr_fft(uz, maxlags=None, norm='coeff')
     self.data['r12'],self.data['taur12'] = tt.xcorr_fft(ux,y=uy, maxlags=None, norm='coeff')
     self.data['r13'],self.data['taur13'] = tt.xcorr_fft(ux,y=uz, maxlags=None, norm='coeff')
     self.data['r23'],self.data['taur23'] = tt.xcorr_fft(uy,y=uz, maxlags=None, norm='coeff')
     self.data['rmag'],self.data['taurmag'] = tt.xcorr_fft(umag, maxlags=None, norm='coeff')
     # auto correlation of u
     self.data['R11'],self.data['tauR11'] = tt.xcorr_fft(ux, maxlags=None, norm='biased')
     self.data['R22'],self.data['tauR22'] = tt.xcorr_fft(uy, maxlags=None, norm='biased')
     self.data['R33'],self.data['tauR33'] = tt.xcorr_fft(uz, maxlags=None, norm='biased')

    def __init__(self, path):
        # convertir el archivo a objeto python
        self.doc2event(path)
        
        # @warning: se requieren convertir las matrices y los atributos que
        # tienen puntos a arrays y a atributos de estructuras respectivamenente

        #convertimos todos los sismos a campos de desplazamiento
        
        acelerometers_id = [76, 82, 118, 126, 146, 147]
        
        
        #transformar todo a campo de desplazamiento
        for s in self.seismograms:
            
            #guardar los datos crudos para analisis posteriores
            s.raw_data = s.data
            # si es acelerometro
            if s.site_id in acelerometers_id:
                s.data = sig.detrend(np.cumsum(sig.detrend(np.cumsum(s.data[
                               :,2:5], axis=0), axis=0), axis=0), axis=0)    
            # si es velocimetro
            else:
                s.data = sig.detrend(np.cumsum(s.data[:, 2:5], axis=0), axis=0)
            
            timevector = s.timevector
            s.data = pd.DataFrame(s.data, index = timevector)
            s.raw_data = pd.DataFrame(s.raw_data[:,2:5], index = timevector)
            
        pass

def phases_from_complex(wts, continuous=False, do_detrend=False):
    """Calculates phases from 1d or 2d wavelet/hilbert arrays, dim0 is time"""
    if len(wts.shape) == 1:
        #1d
        phasen = n.arctan2(wts.imag,wts.real)
        if not (continuous or do_detrend):
            return phasen
        else:
            phasen = make_phases_continuous(phasen)
            if do_detrend:
                phasen = detrend(phasen,axis=0)
            return phasen
    elif len(wts.shape) == 2:
        #2d
        phasen = n.arctan2(wts.imag,wts.real)
        if not (continuous or do_detrend):
            return phasen
        else:
            phasen = make_phases_continuous(phasen)
            if do_detrend:
                phasen = detrend(phasen,axis=0)
            return phasen
        
    else:
        raise ValueError("Only 1d and 2d arrays supported")

    def preProcess(self,
                                    periodF0 = 0.06,
                                    deltaF_div_F0 = True,
                                    
                                    max_threshold = None,
                                    min_threshold = None,
                                    nan_to_zeros = True,
                                    
                                    detrend = False,
                                    
                                    #~ band_filter = None,
                                    
                                    gaussian_filter = None,
                                    
                                    f1 = None,
                                    f2 = None,
                                    
                                    **kargs):
        
        images = self.images
        if deltaF_div_F0:
            ind = self.t()<=self.t_start+periodF0
            m0 = mean(images[ind,:,:] , axis = 0)
            images = (images-m0)/m0*1000.
            
        if max_threshold is not None:
            #~ images[images>max_threshold] = max_threshold
            images[images>max_threshold] = nan
            

        if min_threshold is not None:
            #~ images[images<min_threshold] = min_threshold
            images[images<min_threshold] = nan
                
            
        if nan_to_zeros:
            images[isnan(images) ] = 0.

        if detrend and not nan_to_zeros:
            m = any(isnan(images) , axis = 0)
            images[isnan(images) ] = 0.
            images = signal.detrend( images , axis = 0)
            images[:,m] = nan
        elif detrend and nan_to_zeros:
            images = signal.detrend( images , axis = 0)
            
        if gaussian_filter is not None:
            images = ndimage.gaussian_filter( images , (0 , gaussian_filter , gaussian_filter))
            

        if f1 is not None or f2 is not None:
            from ..computing.filter import fft_passband_filter
            if f1 is None: f1=0.
            if f2 is None: f1=inf
            nq = self.sampling_rate/2.
            images = fft_passband_filter(images, f_low = f1/nq , f_high = f2/nq , axis = 0)
        
        return images

def preprocess(matr, prepr, Fs, fc_min, fc_max, taper_fract):
    """
    :type matr: numpy.ndarray
    :param matr: time series of used stations (dim: [number of samples, number of stations])
    :type prepr: integer
    :param prepr: type of preprocessing. 0=None, 1=bandpass filter, 2=spectral whitening
    :type Fs: float
    :param Fs: sampling rate of data streams
    :type fc_min, fc_max: float
    :param fc_min, fc_max: corner frequencies used for preprocessing
    :type taper_fract: float
    :param taper_fract: percentage of frequency band which is tapered after spectral whitening

    :return: preprocessed data (dim: [number of samples, number of stations])
    """
    if prepr == 0:
        data = signal.detrend(matr, axis=0)

    elif prepr == 1:
        # generate frequency vector and butterworth filter
        b, a = signal.butter(4, np.array([fc_min, fc_max]) / Fs * 2, btype="bandpass")
        # filter data and normalize it by maximum energy
        data = signal.filtfilt(b, a, signal.detrend(matr, axis=0), axis=0)
        fact = np.sqrt(np.dot(np.ones((data.shape[0], 1)), np.sum(data**2, axis=0).reshape((1, data.shape[1]))))
        data = np.divide(data, fact)

    elif prepr == 2:
        nfft = nearest_powof2(matr.shape[0])
        Y = np.fft.fft(matr, n=nfft, axis=0)
        f = np.fft.fftfreq(nfft, 1./float(Fs))

        # whiten: discard all amplitude information within range fc
        Y_white = np.zeros(Y.shape)
        J = np.where((f > fc_min) & (f < fc_max))
        Y_white[J, :] = np.exp(1j * np.angle(Y[J, :]))

        # now taper within taper_fract
        deltaf = (fc_max - fc_min) * taper_fract
        Jdebut = np.where((f > fc_min) & (f < (fc_min + deltaf)))
        Jfin = np.where((f > (fc_max - deltaf)) & (f < fc_max))
        for ii in range(Y.shape[1]):
            if len(Jdebut[0]) > 1:
                Y_white[Jdebut, ii] = np.multiply(Y_white[Jdebut, ii],
                            np.sin(np.pi / 2 * np.arange(0, len(Jdebut[0])) / len(Jdebut[0]))**2)
            if len(Jfin[0]) > 1:
                Y_white[Jfin, ii] = np.multiply(Y_white[Jfin, ii],
                            np.cos(np.pi / 2 * np.arange(0, len(Jfin[0])) / len(Jfin[0]))**2)

        # perform inverse fft to obtain time signal
        # data = 2*np.real(np.fft.ifft(Y_white, n=nfft, axis=0))
        data = np.fft.ifft(Y_white, n=nfft, axis=0)
        # normalize it by maximum energy
        fact = np.sqrt(np.dot(np.ones((data.shape[0], 1)), np.sum(data**2, axis=0).reshape((1, data.shape[1]))))
        data = np.divide(data, fact)
    return data

    def generateStatistics(self,doDetrend=True):
        '''
        Generates statistics and populates member variable data.

        Arguments:
            doDetrend: detrend data bevor sigbal processing

        Populates the "data" python dict with with the following keys:
            rii:    [numpy.array of shape=(?)] Auto-correlation coefficent rii. For i=1,2,3
            taurii: [numpy.array of shape=(?)] Time lags for rii. For i=1,2,3
            Rii:    [numpy.array of shape=(?)] Auto-correlation Rii. For i=1,2,3
            tauRii: [numpy.array of shape=(?)] Time lags for Rii. For i=1,2,3

            uifrq:  [numpy.array of shape=(?)] u1 in frequency domain. For i=1,2,3
            uiamp:  [numpy.array of shape=(?)] amplitude of u1 in frequency domain. For i=1,2,3
            Seiifrq:[numpy.array of shape=(?)] Frequencies for energy spectrum Seii. For i=1,2,3
            Seii:   [numpy.array of shape=(?)] Energy spectrum Seii derived from Rii. For i=1,2,3
        '''
        # auto correlation corefficient of u
        if doDetrend:
            ux=signal.detrend(self.ux());
            uy=signal.detrend(self.uy());
            uz=signal.detrend(self.uz());
            umag=signal.detrend(self.Umag());
        else:
            ux=self.ux();
            uy=self.uy();
            uz=self.uz();
            umag=self.Umag();
        #ux=ux[-samples:-1]
        #uy=uy[-samples:-1]
        #uz=uz[-samples:-1]
        self.data['r11'],self.data['taur11'] = tt.xcorr_fft(ux, maxlags=None, norm='coeff')
        self.data['r22'],self.data['taur22'] = tt.xcorr_fft(uy, maxlags=None, norm='coeff')
        self.data['r33'],self.data['taur33'] = tt.xcorr_fft(uz, maxlags=None, norm='coeff')
        self.data['r12'],self.data['taur12'] = tt.xcorr_fft(ux,y=uy, maxlags=None, norm='coeff')
        self.data['r13'],self.data['taur13'] = tt.xcorr_fft(ux,y=uz, maxlags=None, norm='coeff')
        self.data['r23'],self.data['taur23'] = tt.xcorr_fft(uy,y=uz, maxlags=None, norm='coeff')
        self.data['rmag'],self.data['taurmag'] = tt.xcorr_fft(umag, maxlags=None, norm='coeff')
        # auto correlation of u
        self.data['R11'],self.data['tauR11'] = tt.xcorr_fft(ux, maxlags=None, norm='none')
        self.data['R22'],self.data['tauR22'] = tt.xcorr_fft(uy, maxlags=None, norm='none')
        self.data['R33'],self.data['tauR33'] = tt.xcorr_fft(uz, maxlags=None, norm='none')


        #u in frequency domain
        self.data['u1frq'],self.data['u1amp'] = tt.dofft(sig=ux,samplefrq=self.data['frq'])
        self.data['u2frq'],self.data['u2amp'] = tt.dofft(sig=uy,samplefrq=self.data['frq'])
        self.data['u3frq'],self.data['u3amp'] = tt.dofft(sig=uz,samplefrq=self.data['frq'])
        #Time energy sectrum Se11 (mean: Rii in frequency domain...)
        self.data['Se11frq'],self.data['Se11'] = tt.dofft(sig=self.data['R11'],samplefrq=self.data['frq'])
        self.data['Se22frq'],self.data['Se22'] = tt.dofft(sig=self.data['R22'],samplefrq=self.data['frq'])
        self.data['Se33frq'],self.data['Se33'] = tt.dofft(sig=self.data['R33'],samplefrq=self.data['frq'])

def ts_ft(data, name):
	data /= np.max(np.abs(data), axis=0) 
	plt.subplot(311)
	plt.plot(ma.mean(ma.mean(data, 2), 1))
	plt.axis('tight')
	plt.subplot(312)
	plt.plot(signal.detrend(ma.mean(ma.mean(data, 2), 1), axis=0))
	plt.axis('tight')
	plt.subplot(313)
	plt.plot(np.log(np.fft.rfft(ma.mean(ma.mean(signal.detrend(data, axis=0), 2), 1)))[0:365])
	plt.axis('tight')
	plt.savefig(FIGDIR+'time_freq_'+name+'.png')
	plt.close('all')

def summarize_timeseries(functional_path, masks_path, summary):

    if type(summary) is not dict:
        summary = {'method': summary}

    masks_img = [nb.load(mask_path) for mask_path in masks_path]
    mask = np.sum(np.array([
        mask_img.get_data() for mask_img in masks_img
    ]), axis=0) > 0.0

    if mask.sum() == 0:
        raise Exception(
            "The provided mask does not contains voxels. "
            "Please check if mask is being eroded and if the segmentation worked correctly."
        )

    functional_img = nb.load(functional_path)
    masked_functional = functional_img.get_data()[mask]

    regressors = np.zeros(masked_functional.shape[-1])

    if summary['method'] == 'Mean':
        regressors = masked_functional.mean(0)

    if summary['method'] == 'NormMean':
        masked_functional /= np.linalg.norm(masked_functional, 2)
        regressors = np.nan_to_num(masked_functional).mean(0)

    if summary['method'] == 'DetrendNormMean':
        masked_functional = \
            signal.detrend(masked_functional, type='linear').T

        masked_functional /= np.linalg.norm(masked_functional, 2)
        regressors = np.nan_to_num(masked_functional).mean(0)

    if summary['method'] in ['DetrendPC', 'PC']:
        if summary['method'] == 'DetrendPC':
            Y = signal.detrend(masked_functional, type='linear').T
        else:
            Y = masked_functional.T

        Yc = Y - np.tile(Y.mean(0), (Y.shape[0], 1))
        Yc = np.nan_to_num(Yc / np.tile(np.array(Y.std(0)).reshape(1,Y.shape[1]), (Y.shape[0],1)))
        U, _, _ = np.linalg.svd(Yc)

        regressors = U[:, 0:summary['components']]

    output_file_path = os.path.join(os.getcwd(), 'summary_regressors.1D')
    np.savetxt(output_file_path, regressors, fmt='%.18f')

    return output_file_path

def standardize(X,stdtype='column'):
	"""Standardizes a two dimensional input matrix X, by either row or column.  Resulting matrix will have 
	row or column mean 0 and row or column std. dev. equal to 1.0."""
	if len(X.shape) > 2:
		print 'ERROR: standardize() not defines for matrices that are not two-dimenional'
		return
	if stdtype == 'column':
		F = signal.detrend(X,type='constant',axis=0)
		F = F/std(F,axis=0)
	else:
	    F = signal.detrend(X.T,type='constant',axis=0)
	    F = F/std(F,axis=0)
	    F = F.T
	return F

    def generateAutoCorrelations(self,doDetrend=True):
        # auto correlation corefficient of u
        if doDetrend:
            ux=signal.detrend(self.ux());
            uy=signal.detrend(self.uy());
            uz=signal.detrend(self.uz());
        else:
            ux=self.ux();
            uy=self.uy();
            uz=self.uz();

        self.data['r11'],self.data['taur11'] = tt.xcorr_fft(ux, maxlags=None, norm='coeff')
        self.data['r22'],self.data['taur22'] = tt.xcorr_fft(uy, maxlags=None, norm='coeff')
        self.data['r33'],self.data['taur33'] = tt.xcorr_fft(uz, maxlags=None, norm='coeff')

    def applyMeth(self, x, fMeth):
        """Apply the methods
        x : array signal
        fMeth : list of methods
        -> 3D array of the transform signal
        """
        npts, ntrial = x.shape
        nFce = len(fMeth)
        xf = n.zeros((nFce, npts, ntrial))

        # Detrend the signal :
        if self.dtrd:
            x = detrend(x, axis=0)

        # Apply methods :
        for k in range(0, nFce):  # For each frequency in the tuple
            xf[k, ...] = fMeth[k](x)

        # Correction for the wavelet (due to the wavelet width):
        if (self.method == 'wavelet') and (self.wltCorr is not None):
            w = 3*self.wltWidth
            xf[:, 0:w, :] = xf[:, w+1:2*w+1, :]
            xf[:, npts-w:npts, :] = xf[:, npts-2*w-1:npts-w-1, :]

        return xf

def _iter_contrasts(n_subjects, factor_levels, effect_picks):
    """ Aux Function: Setup contrasts """
    sc, sy, = [], []

    # prepare computation of Kronecker products
    for n_levels in factor_levels:
        # for each factor append
        # 1) column vector of length == number of levels,
        # 2) square matrix with diagonal == number of levels

        # main + interaction effects for contrasts
        sc.append([np.ones([n_levels, 1]),
                   detrend(np.eye(n_levels), type='constant')])
        # main + interaction effects for component means
        sy.append([np.ones([n_levels, 1]) / n_levels, np.eye(n_levels)])
        # XXX component means not returned at the moment

    for (c1, c2, c3) in defaults_twoway_rm['iter_contrasts'][effect_picks]:
        # c1 selects the first factors' level in the column vector
        # c3 selects the actual factor
        # c2 selects either its column vector or diag matrix
        c_ = np.kron(sc[0][c1], sc[c3][c2])
        # for 3 way anova accumulation of c_ across factors required
        df1 = matrix_rank(c_)
        df2 = df1 * (n_subjects - 1)
        yield c_, df1, df2

def enframe(x, win, inc):
    """
    Splits the vector up into (overlapping) frames beginning at increments
    of inc. Each frame is multiplied by the window win().
    The length of the frames is given by the length of the window win().
    The centre of frame I is x((I-1)*inc+(length(win)+1)/2) for I=1,2,...

    :param x: signal to split in frames
    :param win: window multiplied to each frame, length determines frame length
    :param inc: increment to shift frames, in samples
    :return f: output matrix, each frame occupies one row
    :return length, no_win: length of each frame in samples, number of frames
    """
    nx = len(x)
    nwin = len(win)
    if (nwin == 1):
        length = win
    else:
        # length = next_pow_2(nwin)
        length = nwin
    nf = int(fix((nx - length + inc) // inc))
    # f = np.zeros((nf, length))
    indf = inc * np.arange(nf)
    inds = np.arange(length) + 1
    f = x[(np.transpose(np.vstack([indf] * length)) +
           np.vstack([inds] * nf)) - 1]
    if (nwin > 1):
        w = np.transpose(win)
        f = f * np.vstack([w] * nf)
    f = signal.detrend(f, type='constant')
    no_win, _ = f.shape
    return f, length, no_win