Python signal.find_peaks_cwt函数代码示例

    def test_frequency(self):
        import scipy.signal as ss
        nsteps = 10000
        dt = 0.1

        t,w = self.regular_integrator.run(self.regular_w0, dt=dt, nsteps=nsteps)
        f,fft = fft_orbit(t, w)

        peak_ix = ss.find_peaks_cwt(fft[:,0], widths=np.linspace(dt*2, dt*100, 10))
        print(peak_ix)

        plt.clf()
        plt.axvline(self.regular_par[1]/(2*np.pi), linewidth=3., alpha=0.35, color='b')
        plt.axvline(1/(2*np.pi), linewidth=3., alpha=0.35, color='r')
        plt.semilogx(f[:,0], fft[:,0], marker=None)
        plt.savefig(os.path.join(plot_path,"pend_fft_regular.png"))

        # ----------------------------------------------------------------------
        t,w = self.chaotic_integrator.run(self.chaotic_w0, dt=dt, nsteps=nsteps)
        f,fft = fft_orbit(t, w)

        peak_ix = ss.find_peaks_cwt(fft[:,0], widths=np.linspace(dt*2, dt*100, 10))
        print(peak_ix)

        plt.clf()
        plt.axvline(self.chaotic_par[1]/(2*np.pi), linewidth=3., alpha=0.35, color='b')
        plt.axvline(1/(2*np.pi), linewidth=3., alpha=0.35, color='r')
        plt.semilogx(f[:,0], fft[:,0], marker=None)
        plt.savefig(os.path.join(plot_path,"pend_fft_chaotic.png"))

  def fft_analysis(self,position=None):
    '''
    Performes Fast fourier analysis on data sets and returns either five most significant frequencies or numpy arrays of 
    frequencies and power spectrum when specific position number is passed as parameter.
    '''

    if position is None:
      freqs_5=np.array([])
      
      for dp, pr in zip(self.data_points, self.prumery):
        time_step=(dp.data_time[-1]-dp.data_time[0])/dp.data.size
        ps=np.abs(np.fft.rfft(dp.data-pr))**2 #Power spectrum
      
        freqs = np.fft.rfftfreq(dp.data.size, time_step)
        peakind = signal.find_peaks_cwt(ps, np.arange(1,50))
        freqs_5 = np.append(freqs_5,freqs[peakind[:5]])
        
      return freqs[peakind[:5]] #Returns first five most significant frequencies
    
    else: #Returns arrays of frequencies and power spectrum and five most significant frequencies.
      for dp,pr in zip(self.data_points, self.prumery):
        if position == dp.position_no:
          time_step=(dp.data_time[-1]-dp.data_time[0])/dp.data.size
          ps=np.abs(np.fft.rfft(dp.data-pr))**2 #Power spectrum
          #ps=np.angle(np.fft.rfft(dp.data-pr)) #Phase spectrum
          #ps=20*np.log10(np.abs(np.fft.rfft(dp.data-pr))) #Amplitude spectrum in [dB]
        
          freqs = np.fft.rfftfreq(dp.data.size, time_step)
          peakind = signal.find_peaks_cwt(ps, np.arange(1,50))
          
          
          return freqs, ps, freqs[peakind[:5]]

def find_threshold(filename, base_offset=0, pos_pol=True, **kws):
    with trace_gen(filename, base_offset=base_offset, **kws) as gen:
        if pos_pol:
            min_data = [event[100:-100].max() for event in gen]
        else:
            min_data = [event[100:-100].min() for event in gen]
    hist = histogram(min_data, bins=2048)
    # second = 2 if pos_pol else -2
    # index = argrelmax(hist[0][hist[0] != 0], order=5)[0][second]
    second = 1 if pos_pol else -1
    try:
        index = find_peaks_cwt(hist[0][hist[0] != 0], arange(10, 40))[second]
    except IndexError:
        index = find_peaks_cwt(hist[0][hist[0] != 0], arange(10, 40))[0]
    return hist[1][hist[0] != 0][index] * 3 / 4, len(min_data)

def find_peaks(x, y, widthrange, rel_threshold=0.1):
    """Peak-finding in a 2d dataset.
    
    Parameters
    ----------
    x, y : array_like
        Input arrays.
    widthrange : tuple
        Lower and upper limit of peak widths to find.
    rel_threshold : float, optional
        Peaks with a height lower than this value times the height of the
        maximum in 'y' are ignored.

    Returns
    -------
    list
        Array indices of where the peaks were found in 'y'.

    See Also
    --------
    scipy.signal.find_peaks_cwt : Peak-finding using a continous wavelet
    transform technique.
    """
    dx = abs(x[1] - x[0])
    minwidth, maxwidth = widthrange
    widths = np.arange(floor(minwidth/dx), ceil(maxwidth/dx))
    peakpos = find_peaks_cwt(y, widths)
    maxy = max(y)
    return [pos for pos in peakpos if y[pos] >= rel_threshold*maxy]

def main(data):
    start = time.time()
    T = 1.0/SAMPLING_RATE # sampling interval
    Fs = 1.0 / T
    ir_data, bpm_data, avg_bpm_data = clean_data(data)
    hr_data = zero_mean(ir_data)

    cutoff_bpm = [50.0, 200.0]
    cutoff_hz = [x/60 for x in cutoff_bpm] # cutoff frequency in HZ
    cutoff = [x/(Fs/2) for x in cutoff_hz]
    [b, a] = signal.butter(2, cutoff, 'bandpass') # 2nd order butterworth filter

    is_valid = check_data_good(ir_data)
    if (is_valid == False):
        return "Heart Rate: Please Place Sensor on Feet"
    else:
        # filter out the noise from signal
        hr_filt = signal.lfilter(b, a, hr_data)

        pks = signal.find_peaks_cwt(hr_filt, np.arange(3, 10))
        num_pks = len(pks)
        beats_from_peaks = num_pks/2
        bpm_from_peaks = beats_from_peaks*60/TIME_SEC
        print("HR Found from Beats is = " + str(bpm_from_peaks) + " BPM")
        time_btw_peaks = sum(np.array(pks[1:num_pks]) - np.array(pks[0:-1]))/(num_pks - 1)
        bpm_from_peaks = SAMPLING_RATE*60/time_btw_peaks/2
        print("HR Found from Time is = " + str(bpm_from_peaks) + " BPM")
        end = time.time()
        print("total time = " + str(end - start))
        return bpm_from_peaks

def find_peak_data(data, n, spec=False):
    indices = find_peaks_cwt(data, arange(30, 80))
    values = [data[i] for i in indices]
    params = [val for i in range(n) for val in [values[i], indices[i], 30]]
    if spec:
        params.extend((0.2, 1))
    return params

    def execute_analysis(self):
        """
        The continuous wavelet transform peak identification algorithm from scipy.signal.find_peaks_cwt().
        In openmsi it is renamed "findpeaks_cwt" with one fewer underscore.
        """

        from scipy import signal as sig

        msidata = self['msidata']    #now in memory as hdf5 cube or np.ndarray
        widths = self['widths']      #should already be numpy ndarray
        min_snr = self['min_snr']

        shape_x = msidata.shape[0]
        shape_y = msidata.shape[1]

        mzindices = []

        for xi in xrange(0, shape_x):
            for yi in xrange(0, shape_y):
                print xi, yi

                # Load the spectrum
                m = msidata[xi, yi, :]

                # find indices of m where peaks are
                peak_indices = sig.find_peaks_cwt(m, widths=widths, wavelet=sig.ricker, min_snr=min_snr)
                mzindices.append(peak_indices)
        return mzindices

 def get_peaks(self, data, domain):
     thresh = 1.1  # TODO: Make this adjustable - issue#2
     peakrange = np.arange(1, 100)
     rawpeaks = np.array(sc_sg.find_peaks_cwt(data, peakrange))
     rawpeaks = rawpeaks[rawpeaks > rawpeaks.min()]  # Remove leftmost peak (false positive at input)
     threshpeaks = rawpeaks[data[rawpeaks] > thresh]
     return np.array([domain[threshpeaks], data[threshpeaks]]).T

def checkerboard_matrix_filtering(similarity_matrix, kernel_width, peak_range):
    """
    Moving the checkerboard matrix over the main diagonal of the similarity matrix one sample at a time.

    :param similarity_matrix:
    :param peak_range: the number of samples in which the peak detection algorithms finds a peak (TODO:clarify)
    :param kernel_width: the size of one quarter of the checkerboard matrix
    :return: peaks and convolution values
    """

    checkerboard_matrix = get_checkerboard_matrix(kernel_width)

    # The values calculated in this step are starting from the 'kernel_width' position and ending
    # at length - kernel_width
    d = []
    for i in range(0, similarity_matrix.shape[0] - 2 * kernel_width):
        base = similarity_matrix[i:i + kernel_width * 2, i:i + kernel_width * 2]
        d.append(np.sum(np.multiply(base, checkerboard_matrix)))

    # The missing values from 0 to kernel_width are calculated here
    top_left_d = []
    for i in range(0, kernel_width):
        base = similarity_matrix[0:i + kernel_width, 0:i + kernel_width]
        top_left_d.append(np.sum(np.multiply(base, checkerboard_matrix[kernel_width - i:, kernel_width - i:])))

    # The missing kernel_width values at the bottom right are set to 0
    convolution_values = top_left_d + d + [0 for i in range(0, kernel_width)]

    peaks = find_peaks_cwt(convolution_values, np.arange(1, peak_range))
    peaks = [0] + peaks + [len(convolution_values)-1]
    return peaks, convolution_values

def call_boundary_peaks(Pb):
    n = len(Pb)
    peaks = find_peaks_cwt(
        Pb, 
        np.array([0.5]), 
        wavelet=None, 
        max_distances=None, 
        gap_thresh=None, 
        min_length=None, 
        min_snr=1, 
        noise_perc=10)

    regions = []
    for i in peaks:
        top = Pb[i]
        x = [l for l in range(1, min(i, 5)) if Pb[i-l] > 0.4*top]
        y = [r for r in range(1, min(n-i, 5)) if Pb[i+r] > 0.4*top]
        start = (i - x[-1]) if len(x) else i
        end = (i + y[-1]) if len(y) else i
        regions.append([start-1, end+1, top])

    if len(regions):
        starts, ends, values = zip(*regions)
    else:
        starts, ends, values = [], [], []

    return np.array(starts), np.array(ends), np.array(values)

def dih_plotter3(dirname,savename,numplot):
    inlist = dih_tablereader(dirname)
    plotlist = inlist[0:numplot]
    colors = iter(cm.rainbow(np.linspace(0,1,len(plotlist)))) #creates color table
    for memberlist in plotlist:
        x = memberlist[0] #x coordinate data
        y = memberlist[1] #y coordinate data
        ysmooth = dih_boxcar(y)
        xshort = x[0:len(ysmooth)]
        peaklist =signal.find_peaks_cwt(ysmooth, np.arange(1,30))#continuous wavelet transformation
        plt.plot(xshort,ysmooth,color = next(colors))
        for num in peaklist:
            plt.plot(xshort[num],ysmooth[num],'gD')#places markers on peaks
        peak = max(ysmooth)
        peaklist2 = [i for i, j in enumerate(ysmooth) if j == peak]
        for num in peaklist2:
            plt.plot(xshort[num],ysmooth[num],'rD')

#finish up plot characteristics
    plt.title('Super Most Awesome Graph!')
    plt.ylabel('Flux')
    plt.xlabel('Time')       
    pylab.ylim([-5,5])
    pylab.xlim([0,6.3])
    plt.savefig(savename)#saves postscript file
    return plotlist

 def curPeakdetect(y, x):
     if ds["type"] == "cwt":
         _maxidx = find_peaks_cwt(y, ds["widths"])
         _max = zip(x[_maxidx], y[_maxidx])
     else:
         _max, _ = peakdetect(y, x, lookahead=ds["lookahead"], delta=ds["delta"])
     return _max

def calculatepeak(step, jump, dwell):
	totalpeak = np.array([])
	totalinterval = np.array([])
	totalfilterpeak = np.array([])
	totalfilterinterval = np.array([])
	for window in step:
		if window/jump % 2 == 0:
			window += jump
		ma = convertma(dwell, window, jump)
		interval = window / jump
		#Calculate the moving average of the result
		std = movstd(ma, interval)
		#Calculate the moving standard deviation of the result
		oripeak = signal.find_peaks_cwt(std, np.arange(interval * percentage, interval, interval*0.05))
		oristdamp = std[oripeak]
		#Obtain the standard deviation peak and the amplitude of the peak

		if len(oripeak) == 0:
			pass
		else:
			oripopendiffpeak = calpopendiffpeak(dwell, oripeak)
			#Calculate the difference in Popen from the raw data using the time of the peak
			oripopendiffstdamp = oristdamp / np.sqrt(1.0/12)
			#Calculate the difference in Popen based on the value of the standard deviation peak
			filterpeak, oripopendiffstdamp, oripopendiffpeak = filterpeaktime(dwell, std, oripeak, oripopendiffstdamp, oripopendiffpeak)

			totalpeak = np.append(totalpeak, oripeak)
			totalinterval = np.append(totalinterval, [interval] * len(oripeak))
			totalfilterpeak = np.append(totalfilterpeak, filterpeak)
			totalfilterinterval = np.append(totalfilterinterval, [interval] * len(filterpeak))
			for iprint in range(len(filterpeak)):
				print 'peak', filterpeak[iprint], 'Popen(avg)', oripopendiffpeak[iprint], 'Popen(std)', oripopendiffstdamp[iprint]

	np.savez(pathfilename+str(step[0])+'_'+str(step[-1]), totalpeak=totalpeak, totalinterval=totalinterval, totalfilterpeak=totalfilterpeak, totalfilterinterval=totalfilterinterval)

def determine_step_process(dwell, percentage, start, end, jump):
	attempt = np.array([])
	firstpeak = np.array([])
	lastpeak = np.array([])
	suminterval = np.array([])
	for window in np.arange(start, (end+jump*2), jump*2):
		#print 'testing interval:', interval,
		#window must be an odd number
		if window/jump % 2 == 0:
			window += jump

		ma = convertma(dwell, window, jump)
		interval = window / jump
		std = movstd(ma, interval)
		oripeak = signal.find_peaks_cwt(std, np.arange(interval * percentage, interval, interval*0.05))
		if oripeak != []:
			peaklength = np.append(oripeak, np.cumsum(dwell)[-1]/jump-1) - np.append(0, oripeak)
			if np.amin(peaklength) > interval:
				attempt = np.append(attempt, 1)
			else:
				attempt = np.append(attempt, 0)
			firstpeak = np.append(firstpeak, oripeak[0])
			lastpeak = np.append(lastpeak, oripeak[-1])
			suminterval = np.append(suminterval, interval)

	np.savez(pathfilename+'_'+str(start)+'_'+str(end), attempt = attempt, firstpeak = firstpeak, lastpeak = lastpeak, suminterval = suminterval)

def sort_spikes(sweep):
    filtered = np.array(sweep) - medfilt(sweep, 51)
    pks = find_peaks_cwt(filtered, np.arange(9, 20))
    offsets = np.zeros_like(pks)
    datamtrx = np.zeros(shape=(80, len(pks[3:])))

    for i, pk in enumerate(pks[3:]):
        offset = np.argmax(filtered[pk - 40 : pk + 40]) - 40
        datamtrx[:, i] = filtered[pk - 40 + offset : pk + 40 + offset]
        offsets[i + 3] = offset
        # plb.plot(filtered[pk-40+offset:pk+20+offset],color = 'k', alpha = 0.1)

    from scipy.linalg import svd

    U, s, Vt = svd(datamtrx, full_matrices=False)
    V = Vt.T

    ind = np.argsort(s)[::-1]
    U = U[:, ind]
    s = s[ind]
    V = V[:, ind]

    features = V
    es, idx = kmeans2(features[:, 0:4], 4, iter=50)
    colors = [([1, 0, 1], [1, 0, 0], [0, 0, 1], [0, 1, 1])[i] for i in idx]