Python numpy.bartlett函数代码示例

 def test_perform(self):
     x = tensor.lscalar()
     f = function([x], self.op(x))
     M = numpy.random.random_integers(3, 50, size=())
     assert numpy.allclose(f(M), numpy.bartlett(M))
     assert numpy.allclose(f(0), numpy.bartlett(0))
     assert numpy.allclose(f(-1), numpy.bartlett(-1))
     b = numpy.array([17], dtype='uint8')
     assert numpy.allclose(f(b[0]), numpy.bartlett(b[0]))

def smooth(x, windowLen, windowType="boxcar"):
    """Smooth a numpy array, returning the smoothed values
    
    Adapted from http://www.scipy.org/Cookbook/SignalSmooth"""
    if x.ndim != 1:
        raise ValueError("smooth only accepts 1 dimension arrays.")
 
    if x.size < windowLen:
        raise ValueError("Input vector needs to be bigger than window size.")

    if windowLen < 3:
        return x

    if windowType == "boxcar":
        w = np.ones(windowLen,'d')
    elif windowType == "hamming":
        w = np.hamming(windowLen)
    elif windowType == "hanning":
        w = np.hanning(windowLen)
    elif windowType == "bartlett":
        w = np.bartlett(windowLen)
    elif windowType == "blackman":
        w = np.blackman(windowLen)
    else:
        raise ValueError("windowType %s is not one of 'boxcar', 'hanning', 'hamming', 'bartlett', 'blackman'"
                         % windowType)

    s = np.r_[x[windowLen-1:0:-1],x,x[-1:-windowLen:-1]]

    y = np.convolve(w/w.sum(), s, mode='valid')

    return y[windowLen//2:-windowLen//2 + 1]

def power_spectrum(trace, time_window=1000, overlap_window=None, enforce_pow2=True, freq_range=None):
# -----------------------------------------------------------------------------
    '''
    Returns the power spectrum of trace. Parameters set to fit with what Paolo 
    does within his project.
    '''
    from matplotlib.mlab import psd
    from matplotlib.pylab import detrend_linear, detrend_mean, detrend_none
    from neurovivo.common import nextpow2
        
    dt = (trace._time[1]-trace._time[0]).rescale(pq.ms).item()
    Fs = 1000./dt
    NFFT=int(time_window/dt)
    if enforce_pow2:
        NFFT=nextpow2(NFFT)
        time_window = 1.*time_window*NFFT/int(time_window/dt)
    if overlap_window == None:
        overlap_window = 0.75 * time_window
    noverlap=int(overlap_window/dt)
    window = np.bartlett(NFFT)
    #print dt, Fs, NFFT, noverlap
    pows, freqs = psd(detrend_mean(trace._data), NFFT=NFFT, Fs=Fs, noverlap=noverlap, window=window, detrend=detrend_none)
    if freq_range == None:
        return pows, freqs
    else:
        freqs, pows = cmn.splice_two_vectors(freqs, pows, freq_range)
        return pows, freqs

def plotWindowFunc():
	plt.clf()
	plt.plot(np.arange(20), np.kaiser(20,3.5))
	plt.plot(np.arange(20), np.bartlett(20))
	plt.plot(np.arange(20), np.blackman(20))
	plt.plot(np.arange(20), np.hamming(20))
	plt.plot(np.arange(20), np.hanning(20))

 def apply_filter(self, array):
     if len(array.shape)!=1:
         return False # need 2D data
     nx = array.shape[0]
     # Apodization function
     if self.type == 1:
         ft = np.bartlett(nx)
     elif self.type == 2:
         ft = np.blackman(nx)
     elif self.type == 3:
         ft = np.hamming(nx)
     elif self.type == 4:
         ft = np.hanning(nx)
     elif self.type == 5:
         # Blackman-Harris window
         a0 = 0.35875
         a1 = 0.48829
         a2 = 0.14128
         a3 = 0.01168
         x = np.linspace(0,1,nx)
         ft = a0 - a1*np.cos(2*np.pi*x) + a2*np.cos(4*np.pi*x) + a3*np.cos(6*np.pi*x) 
     elif self.type == 6:
         # Lanczos window
         x = np.linspace(-1,1,nx)
         ft = np.sinc(x)
     elif self.type == 7:
         x = np.linspace(-1,1,nx)
         exec('x='+self.custom)
         ft = x
     else:
         ft = np.ones(nx)
     
     array.data = array.data*ft
     return True

def single_taper_spectrum(data, delta, taper_name=None):
    """
    Returns the spectrum and the corresponding frequencies for data with the
    given taper.
    """
    length = len(data)
    good_length = length // 2 + 1
    # Create the frequencies.
    # XXX: This might be some kind of hack
    freq = abs(np.fft.fftfreq(length, delta)[:good_length])
    # Create the tapers.
    if taper_name == 'bartlett':
        taper = np.bartlett(length)
    elif taper_name == 'blackman':
        taper = np.blackman(length)
    elif taper_name == 'boxcar':
        taper = np.ones(length)
    elif taper_name == 'hamming':
        taper = np.hamming(length)
    elif taper_name == 'hanning':
        taper = np.hanning(length)
    elif 'kaiser' in taper_name:
        taper = np.kaiser(length, beta=14)
    # Should never happen.
    else:
        msg = 'Something went wrong.'
        raise Exception(msg)
    # Detrend the data.
    data = detrend(data)
    # Apply the taper.
    data *= taper
    spec = abs(np.fft.rfft(data)) ** 2
    return spec, freq

def estimate_unbnd_conc_in_region(
        motif, score_cov, atacseq_cov, chipseq_rd_cov,
        frag_len, max_chemical_affinity_change):
    # trim the read coverage to account for the motif length
    trimmed_atacseq_cov = atacseq_cov[len(motif)+1:]
    chipseq_rd_cov = chipseq_rd_cov[len(motif)+1:]

    # normalzie the atacseq read coverage
    atacseq_weights = trimmed_atacseq_cov/trimmed_atacseq_cov.max()
    
    # build the smoothing window
    sm_window = np.ones(frag_len, dtype=float)/frag_len
    sm_window = np.bartlett(2*frag_len)
    sm_window = sm_window/sm_window.sum()

    def build_occ(log_tf_conc):
        raw_occ = logistic(log_tf_conc + score_cov/(R*T))
        occ = raw_occ*atacseq_weights
        smoothed_occ = np.convolve(sm_window, occ/occ.sum(), mode='same')

        return raw_occ, occ, smoothed_occ

    def calc_lhd(log_tf_conc):
        raw_occ, occ, smoothed_occ = build_occ(-log_tf_conc)
        #diff = (100*smoothed_occ - 100*rd_cov/rd_cov.sum())**2
        lhd = -(np.log(smoothed_occ + 1e-12)*chipseq_rd_cov).sum()
        #print log_tf_conc, diff.sum()
        return lhd

    res = brute(calc_lhd, ranges=(
        slice(0, max_chemical_affinity_change, 1.0),))[0]
    log_tf_conc = max(0, min(max_chemical_affinity_change, res))
                      
    return -log_tf_conc

def smooth(input_data, nth_octave = 6, window_type='hamming'):
    """ Smooth input data over 1/n octave """

    f_min = 30
    f_max = 20e3

    number_of_octaves = math.log(f_max / f_min, 2)

    # ideally, this should be computed from the display resolution
    number_of_points = 4048
    points_per_octave = number_of_points / number_of_octaves

    log_data = _distribute_over_log(input_data, f_min, f_max, 
                                                 number_of_points)

    window_length = points_per_octave / nth_octave

    if window_type == 'hamming':
        window = np.hamming(window_length)
    elif window_type == 'bartlett':
        window = np.bartlett(window_length)
    elif window_type == 'blackman':
        window = np.blackman(window_length)
    elif window_type == 'hanning':
        window = np.hanning(window_length)

    output = np.convolve(window / window.sum(), log_data, mode='same')
    return output

	def bartlett_window(self, show=0):
		apod = N.bartlett(self.data_points)
		for i in range(2):
			self.the_result.y[i] = self.the_result.y[i]*apod
		if show == 1:
			return self.the_result
		return self

def mdist_curves(mdist_real):
    """
    fit curve comparing median location error to the actual location error
    """
    CHUNKS = 10
    dists_ = np.array(list(mdist_real))
    dists = dists_[dists_[:,0].argsort()]

    dist_count = dists.shape[0]
    # cut off the start to make even chunks
    chunks = np.split(dists[(dist_count%CHUNKS):,:],CHUNKS)
    bins = utils.dist_bins(30)

    for index,chunk in enumerate(chunks):
        row = chunk[:,1]
        hist,b = np.histogram(row,bins)

        bin_mids = (b[:-1]+b[1:])/2
        bin_areas = np.pi*(b[1:]**2 - b[:-1]**2)
        scaled_hist = hist/(bin_areas*len(row))
        window = np.bartlett(5)
        smooth_hist=np.convolve(scaled_hist,window,mode='same')/sum(window)
        coeffs = np.polyfit(
                    np.log(bin_mids[1:121]),
                    np.log(smooth_hist[1:121]),
                    3)

        yield dict(
                coeffs = list(coeffs),
                cutoff = 0 if index==0 else chunk[0,0],
                local = scaled_hist[0],
                )

    def SelectWindowsFun(self, index):
        #global window
		x = {0: np.ones(self.CHUNK), 
				1:  np.hamming(self.CHUNK), 
				2: np.hanning(self.CHUNK), 
				3: np.bartlett(self.CHUNK), 
				4: np.blackman(self.CHUNK)}
		self.window = x[index]

def nonlinear(x):
    """
    Nonlinear energy operator for spike detection
    """
    xo = np.int32(x)
    y = [xo[n] ** 2 + xo[n - 1] * xo[n + 1] for n in range(1, len(x) - 1)]
    window = np.bartlett(12)
    return np.convolve(y, window)

def PowerSpectrum2(im, win=2, n1=1, n2=0):
    """2D spectrum estimation using the modified periodogram.
    This one includes a window function to decrease variance in \
    the estimate.
    
    :param x: input sequence
    :type x: numpy.array
    :param n1: starting index, x(n1)
    :type n1: int
    :param n2: ending index, x(n2)
    :type n2: int
    :param win: The window type \n
                1 = Rectangular \n
                2 = Hamming \n
                3 = Hanning \n
                4 = Bartlett \n
                5 = Blackman \n
    :type win: int
    
    :returns: spectrum estimate.
    :rtype: numpy.array
            
    .. note:: 
       If n1 and n2 are not specified the periodogram of the entire 
       sequence is computed.
    
    """
    
    if n2 == 0:
        n2 = len(im[:,1])
    
    N  = n2 - n1 + 1
    w  = np.ones((N))
    
    if (win == 2):
        w = np.hamming(N)
    elif (win == 3):
        w = np.hanning(N)
    elif (win == 4):
        w = np.bartlett(N)
    elif (win == 5): 
        w = np.blackman(N);
    
    
    
    xs, ys = im.shape
    if xs/ys != 1:
        raise ValueError('Dimensions must be equal')
    
    
    m = w[:] * w[:][np.newaxis,:]
    
    U  = np.linalg.norm(w)**2.0 / N**2.0
    
    fftim = np.abs(np.fft.fftshift(np.fft.fft2(((im) * m)))) / ( (N**2.0) * U)
    
    return fftim

    def fid_to_spec(self, window_f = ' ', beta = 2, t_start = 0., t_end = 100.0):
        '''
        FFT with zero-padding of the FID.
        Saves the result in self.spec
        
        Parameters
        ----------
        window_f: is the applied window function. The following window functions are
        applicable: hanning, hamming, blackman, bartlett, kaiser (beta (float) is only used
        by the kaiser window)
        t_start (float): specifies the beginning of the FID. The corresponding data points
        will be cut away.
        t_end (float): specifies the end of the FID. The corresponding data points
        will be cut away.
                
        Returns
        -------
        none
        
        Notes
        -----
        none
        '''
        if len(self.fid) != 0:

            fid = slice_fid(self.fid, t_start, t_end)
            window_f = window_f.lower()   # Choose the window function (default: rect / none)    
            
            if window_f == 'hanning':
                fid[:,1] = fid[:,1] * np.hanning(len(fid))
                
            elif window_f == 'hamming':
                fid[:,1] = fid[:,1] * np.hamming(len(fid))
                
            elif window_f == 'blackman':
                fid[:,1] = fid[:,1] * np.blackman(len(fid))
            
            elif window_f == 'bartlett':
                fid[:,1] = fid[:,1] * np.bartlett(len(fid))
                
            elif window_f == 'kaiser':
                fid[:,1] = fid[:,1] * np.kaiser(len(fid), beta)
            
            h = (int(np.sqrt(len(fid))) + 1) ** 2   # Zero padding to n**2 length to enhance computing speed
            spec = np.absolute(np.fft.rfft(fid[:,1], h)) ** 2 / h
            spec = spec[0:int(h/2)]
            freq = np.fft.fftfreq(h, np.abs(fid[2,0] - fid[1,0])) / 1.E6
            freq = freq[0:int(h/2)]  # Combine FFT and frequencies

            self.spec = np.column_stack([freq, spec])
            
            self.spec_params['npoints'] = len(spec)
            self.spec_params['max_f'] = np.max(freq)
            self.spec_params['min_f'] = np.min(freq)
            self.spec_params['delta_f'] = np.abs(freq[1]-freq[0])

    def applyWindow(self, samples, window='hanning'):

        if window == 'bartlett':
            return samples*np.bartlett(len(samples))
        elif window == 'blackman':
            return samples*np.hanning(len(samples))
        elif window == 'hamming':
            return samples*np.hamming(len(samples))
        elif window == 'kaiser':
            return samples*np.kaiser(len(samples))
        else:
            return samples*np.hanning(len(samples))