Python scipy.signal.find_peaks_cwt() Examples

The following are 6 code examples of scipy.signal.find_peaks_cwt(). You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module scipy.signal , or try the search function .
Example #1
Source File: spectools.py    From specidentify with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def detect_lines(x, y, kernal_size=3, centroid_kernal=default_kernal, 
                 center=False):
    """Detect lines goes through a 1-D spectra and detect peaks

    Parameters
    ----------
    x: ~numpy.ndarray
        Array describing the x-position 

    y: ~numpy.ndarray
        Array describing the counts in each x-position

    kernal_size: int
        Size for the detection kernal

    centroid_kernal: 
        Kernal to be used for centroiding

    center: boolean
        If True, centroid for detected peaks will be calculated

    Returns
    -------
    xp: ~numpy.ndarray
        Array of x-positions of peaks in the spectrum
    """
    # find all peaks
    xp = signal.find_peaks_cwt(y, np.array([kernal_size]))
    xp = np.array(xp)

    # set the output values
    if center:
        xdiff = int(0.5 * len(centroid_kernal) + 1)
        x_arr = np.arange(len(x))
        xp = xp * 1.0
        for i in range(len(xp)):
            xp[i] = mcentroid(x, y, kern=centroid_kernal, xdiff=xdiff, xc=x[xp[i]])

    return xp 
Example #2
Source File: peaksearch.py    From DynaPhoPy with MIT License 5 votes vote down vote up
def get_frequencies_from_correlation(correlation_vector,test_frequencies_range):

    frequencies = []
    for branch in range(correlation_vector.shape[1]):
        peakind = signal.find_peaks_cwt(correlation_vector[:,branch].real, np.arange(1,200) )

     #   plt.plot(test_frequencies_range,correlation_vector[:,branch].real)
     #   plt.plot([test_frequencies_range[i] for i in peakind],[correlation_vector[i,branch].real for i in peakind],'ro')
     #   plt.show()

        heights = [correlation_vector[i,branch] for i in peakind]
        max_height_index = heights.index(max(heights))
        frequencies.append(test_frequencies_range[peakind[max_height_index]])

    return np.array(frequencies) 
Example #3
Source File: stats.py    From enlopy with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def get_peaks(x, n):
    return find_peaks_cwt(x, widths=np.arange(1, n + 1), wavelet=ricker) 
Example #4
Source File: elementaldistribution.py    From dscribe with Apache License 2.0 5 votes vote down vote up
def test_single_continuous_property(self):
        # Tested on a water molecule
        system = ase.build.molecule("H2O")

        # Descriptor setup
        std = 0.1
        elements = ["H", "O"]
        peaks = [0.3, 2.0]
        values = dict(zip(elements, peaks))
        elemdist = ElementalDistribution(
            properties={
                "first_property": {
                    "type": "continuous",
                    "min": 0,
                    "max": 2.5,
                    "std": std,
                    "n": 50,
                    "values": values
                }
            }
        )

        # Features
        y = elemdist.create(system)
        y = y.todense().A1
        x = elemdist.get_axis("first_property")

        # Test that the peak positions match
        from scipy.signal import find_peaks_cwt
        peak_indices = find_peaks_cwt(y, [std])
        peak_loc = x[peak_indices]

        # Test that the peak locations match within some tolerance
        self.assertTrue(np.allclose(peaks, peak_loc, rtol=0, atol=0.05))

        # Plot for visual inspection
        # mpl.plot(x, y)
        # mpl.show() 
Example #5
Source File: encode_task_fraglen_stat_pe.py    From atac-seq-pipeline with MIT License 4 votes vote down vote up
def fragment_length_qc(data, prefix):
    results = []

    NFR_UPPER_LIMIT = 150
    MONO_NUC_LOWER_LIMIT = 150
    MONO_NUC_UPPER_LIMIT = 300

    # % of NFR vs res
    nfr_reads = data[data[:, 0] < NFR_UPPER_LIMIT][:, 1]
    percent_nfr = nfr_reads.sum() / data[:, 1].sum()
    results.append(
        QCGreaterThanEqualCheck('Fraction of reads in NFR', 0.4)(percent_nfr))

    # % of NFR vs mononucleosome
    mono_nuc_reads = data[
        (data[:, 0] > MONO_NUC_LOWER_LIMIT) &
        (data[:, 0] <= MONO_NUC_UPPER_LIMIT)][:, 1]

    percent_nfr_vs_mono_nuc = (
        nfr_reads.sum() /
        mono_nuc_reads.sum())
    results.append(
        QCGreaterThanEqualCheck('NFR / mono-nuc reads', 2.5)(
            percent_nfr_vs_mono_nuc))

    # peak locations
    pos_start_val = data[0, 0]  # this may be greater than 0
    peaks = find_peaks_cwt(data[:, 1], np.array([25]))
    nuc_range_metrics = [
        ('Presence of NFR peak', 20 - pos_start_val, 90 - pos_start_val),
        ('Presence of Mono-Nuc peak',
         120 - pos_start_val, 250 - pos_start_val),
        ('Presence of Di-Nuc peak',
         300 - pos_start_val, 500 - pos_start_val)]
    for range_metric in nuc_range_metrics:
        results.append(QCHasElementInRange(*range_metric)(peaks))

    out = prefix + '.nucleosomal.qc'
    with open(out, 'w') as fp:
        for elem in results:
            fp.write(
                '\t'.join(
                    [elem.metric, str(elem.qc_pass), elem.message]) + '\n')

    return out 
Example #6
Source File: encode_task_fraglen_stat_pe.py    From chip-seq-pipeline2 with MIT License 4 votes vote down vote up
def fragment_length_qc(data, prefix):
    results = []

    NFR_UPPER_LIMIT = 150
    MONO_NUC_LOWER_LIMIT = 150
    MONO_NUC_UPPER_LIMIT = 300

    # % of NFR vs res
    nfr_reads = data[data[:, 0] < NFR_UPPER_LIMIT][:, 1]
    percent_nfr = nfr_reads.sum() / data[:, 1].sum()
    results.append(
        QCGreaterThanEqualCheck('Fraction of reads in NFR', 0.4)(percent_nfr))

    # % of NFR vs mononucleosome
    mono_nuc_reads = data[
        (data[:, 0] > MONO_NUC_LOWER_LIMIT) &
        (data[:, 0] <= MONO_NUC_UPPER_LIMIT)][:, 1]

    percent_nfr_vs_mono_nuc = (
        nfr_reads.sum() /
        mono_nuc_reads.sum())
    results.append(
        QCGreaterThanEqualCheck('NFR / mono-nuc reads', 2.5)(
            percent_nfr_vs_mono_nuc))

    # peak locations
    pos_start_val = data[0, 0]  # this may be greater than 0
    peaks = find_peaks_cwt(data[:, 1], np.array([25]))
    nuc_range_metrics = [
        ('Presence of NFR peak', 20 - pos_start_val, 90 - pos_start_val),
        ('Presence of Mono-Nuc peak',
         120 - pos_start_val, 250 - pos_start_val),
        ('Presence of Di-Nuc peak',
         300 - pos_start_val, 500 - pos_start_val)]
    for range_metric in nuc_range_metrics:
        results.append(QCHasElementInRange(*range_metric)(peaks))

    out = prefix + '.nucleosomal.qc'
    with open(out, 'w') as fp:
        for elem in results:
            fp.write(
                '\t'.join(
                    [elem.metric, str(elem.qc_pass), elem.message]) + '\n')

    return out