Python scipy.signal.find_peaks() Examples

def detect_impulses(x, fs, k=10, tdist=1e-3):
    """Detect impulses in `x`

    The minimum height of impulses is defined by `a+k*b`
    where `a` is median of the envelope of `x` and `b` is median 
    absolute deviation (MAD) of the envelope of `x`.   

    :param x: real signal
    :param fs: sampling frequency in Hz
    :param k: multiple of MAD for the impulse minimum height (default: 10)
    :param tdist: minimum time difference between neighbouring impulses in sec (default: 1e-3)
    :returns: indices and heights of detected impulses

    >>> nsamp = 1000
    >>> ind_impulses = np.array([10, 115, 641, 888])
    >>> x = np.zeros((nsamp))
    >>> x[ind_impulses] = 1
    >>> x += np.random.normal(0, 0.1, nsamp)
    >>> ind_pks, h_pks = signal.detect_impulses(x, fs=100000, k=10, tdist=1e-3)
    """
    env = envelope(x)
    height = _np.median(env)+k*_np.median(_np.abs(env-_np.median(env)))
    distance = int(tdist*fs)
    ind_imp, properties = _sig.find_peaks(env, height=height, distance=distance)
    return ind_imp, properties["peak_heights"] 

def FeatureSpectralTonalPowerRatio(X, f_s, G_T=5e-4):

    X = X**2

    fSum = X.sum(axis=0)
    vtpr = np.zeros(fSum.shape)

    for n in range(0, X.shape[1]):
        if fSum[n] < G_T:
            continue

        # find local maxima above the threshold
        afPeaks = find_peaks(X[:, n], height=G_T)

        if not afPeaks[0].size:
            continue

        # calculate ratio
        vtpr[n] = X[afPeaks[0], n].sum() / fSum[n]

    return (vtpr) 

def preprocess(data, config):
    sr = config.sample_rate
    if sr == None:
      sr = 300
    data = np.nan_to_num(data) # removing NaNs and Infs
    from scipy.signal import resample
    data = resample(data, int(len(data) * 360 / sr) ) # resample to match the data sampling rate 360(mit), 300(cinc)
    from sklearn import preprocessing
    data = preprocessing.scale(data)
    from scipy.signal import find_peaks
    peaks, _ = find_peaks(data, distance=150)
    data = data.reshape(1,len(data))
    data = np.expand_dims(data, axis=2) # required by Keras
    return data, peaks

# predict 

def _calculate_results(self, s, interval, df, loader):
        if not df.empty:
            hist = pd.cut(df[Names.HEART_RATE], np.arange(30, 90), right=False).value_counts(sort=False)
            peaks, _ = find_peaks(hist)
            for peak in peaks:
                rest_hr = hist.index[peak].left
                measurements = hist.loc[rest_hr]
                if measurements > len(df) * 0.01:
                    log.debug(f'Rest HR is {rest_hr} with {measurements} values')
                    # conversion to int as value above is numpy int64
                    loader.add(Titles.REST_HR, Units.BPM, S.join(S.MIN, S.MSR), interval,
                               int(rest_hr), local_date_to_time(interval.start), StatisticJournalInteger,
                               'The rest heart rate')
                    return
                else:
                    log.debug(f'Skipping rest HR at {rest_hr} because too few measurements ({measurements}/{len(df)})')
        log.warning(f'Unable to calculate rest HR at {format_date(interval.start)}') 

def infer_pitch_v2(pitch, shortest=5, offset_interval=6):
    w_on = pitch[:,2]
    w_dura = pitch[:,1]

    peaks, properties = find_peaks(w_on, distance=shortest, width=5)
    if len(peaks) == 0:
        return []

    notes = []
    adjust = 5 if shortest==10 else 2
    for i in range(len(peaks)-1):
        notes.append({"start": peaks[i]-adjust, "end": peaks[i+1]-adjust, "stren": pitch[peaks[i], 2]})
    notes.append({"start": peaks[-1]-adjust, "end": len(w_on), "stren": pitch[peaks[-1], 2]})

    del_idx = []
    for idx, p in enumerate(peaks):
        upper = int(peaks[idx+1]) if idx<len(peaks)-2 else len(w_dura)
        for i in range(p, upper):
            if np.sum(w_dura[i:i+offset_interval]) == 0:
                if i - notes[idx]["start"] < shortest:
                    del_idx.append(idx)
                else:
                    notes[idx]["end"] = i
                break

    for ii, i in enumerate(del_idx):
        del notes[i-ii]

    return notes 

def test_k2_peaks_finite(self):
        """Tests the correct peak locations and intensities are found for the
        k=2 term in finite systems.
        """
        desc = MBTR(
            species=[1, 8],
            k2={
                "geometry": {"function": "distance"},
                "grid": {"min": -1, "max": 3, "sigma": 0.5, "n": 1000},
                "weighting": {"function": "unity"},
            },
            normalize_gaussians=False,
            periodic=False,
            flatten=True,
            sparse=False
        )
        features = desc.create(H2O)[0, :]
        pos = H2O.get_positions()
        x = desc.get_k2_axis()

        # Check the H-H peaks
        hh_feat = features[desc.get_location(("H", "H"))]
        hh_peak_indices = find_peaks(hh_feat, prominence=0.5)[0]
        hh_peak_locs = x[hh_peak_indices]
        hh_peak_ints = hh_feat[hh_peak_indices]
        self.assertTrue(np.allclose(hh_peak_locs, [np.linalg.norm(pos[0] - pos[2])], rtol=0, atol=1e-2))
        self.assertTrue(np.allclose(hh_peak_ints, [1], rtol=0, atol=1e-2))

        # Check the O-H peaks
        ho_feat = features[desc.get_location(("H", "O"))]
        ho_peak_indices = find_peaks(ho_feat, prominence=0.5)[0]
        ho_peak_locs = x[ho_peak_indices]
        ho_peak_ints = ho_feat[ho_peak_indices]
        self.assertTrue(np.allclose(ho_peak_locs, np.linalg.norm(pos[0] - pos[1]), rtol=0, atol=1e-2))
        self.assertTrue(np.allclose(ho_peak_ints, [2], rtol=0, atol=1e-2))

        # Check that everything else is zero
        features[desc.get_location(("H", "H"))] = 0
        features[desc.get_location(("H", "O"))] = 0
        self.assertEqual(features.sum(), 0) 

def test_k1_peaks_finite(self):
        """Tests the correct peak locations and intensities are found for the
        k=1 term.
        """
        desc = MBTR(
            species=[1, 8],
            k1={
                "geometry": {"function": "atomic_number"},
                "grid": {"min": 0, "max": 9, "sigma": 0.5, "n": 1000}
            },
            normalize_gaussians=False,
            periodic=False,
            flatten=True,
            sparse=False
        )
        features = desc.create(H2O)[0, :]
        x = desc.get_k1_axis()

        # Check the H peaks
        h_feat = features[desc.get_location(("H"))]
        h_peak_indices = find_peaks(h_feat, prominence=1)[0]
        h_peak_locs = x[h_peak_indices]
        h_peak_ints = h_feat[h_peak_indices]
        self.assertTrue(np.allclose(h_peak_locs, [1], rtol=0, atol=1e-2))
        self.assertTrue(np.allclose(h_peak_ints, [2], rtol=0, atol=1e-2))

        # Check the O peaks
        o_feat = features[desc.get_location(("O"))]
        o_peak_indices = find_peaks(o_feat, prominence=1)[0]
        o_peak_locs = x[o_peak_indices]
        o_peak_ints = o_feat[o_peak_indices]
        self.assertTrue(np.allclose(o_peak_locs, [8], rtol=0, atol=1e-2))
        self.assertTrue(np.allclose(o_peak_ints, [1], rtol=0, atol=1e-2))

        # Check that everything else is zero
        features[desc.get_location(("H"))] = 0
        features[desc.get_location(("O"))] = 0
        self.assertEqual(features.sum(), 0) 

def test_k2_peaks_finite(self):
        """Tests the correct peak locations and intensities are found for the
        k=2 term in finite systems.
        """
        desc = LMBTR(
            species=["H", "O"],
            k2={
                "geometry": {"function": "distance"},
                "grid": {"min": -1, "max": 3, "sigma": 0.5, "n": 1000},
                "weighting": {"function": "unity"},
            },
            normalize_gaussians=False,
            periodic=False,
            flatten=True,
            sparse=False
        )
        features = desc.create(H2O, [0])[0, :]
        pos = H2O.get_positions()
        x = desc.get_k2_axis()

        # Check the X-H peaks
        xh_feat = features[desc.get_location(("X", "H"))]
        xh_peak_indices = find_peaks(xh_feat, prominence=0.5)[0]
        xh_peak_locs = x[xh_peak_indices]
        xh_peak_ints = xh_feat[xh_peak_indices]
        self.assertTrue(np.allclose(xh_peak_locs, [np.linalg.norm(pos[0] - pos[2])], rtol=0, atol=1e-2))
        self.assertTrue(np.allclose(xh_peak_ints, [1], rtol=0, atol=1e-2))

        # Check the X-O peaks
        xo_feat = features[desc.get_location(("X", "O"))]
        xo_peak_indices = find_peaks(xo_feat, prominence=0.5)[0]
        xo_peak_locs = x[xo_peak_indices]
        xo_peak_ints = xo_feat[xo_peak_indices]
        self.assertTrue(np.allclose(xo_peak_locs, np.linalg.norm(pos[0] - pos[1]), rtol=0, atol=1e-2))
        self.assertTrue(np.allclose(xo_peak_ints, [1], rtol=0, atol=1e-2))

        # Check that everything else is zero
        features[desc.get_location(("X", "H"))] = 0
        features[desc.get_location(("X", "O"))] = 0
        self.assertEqual(features.sum(), 0) 

def add_noise(config):
    noises = dict()
    noises["trainset"] = list()
    noises["testset"] = list() 
    import csv
    try:
        testlabel = list(csv.reader(open('training2017/REFERENCE.csv')))
    except:
        cmd = "curl -O https://archive.physionet.org/challenge/2017/training2017.zip"
        os.system(cmd)
        os.system("unzip training2017.zip")
        testlabel = list(csv.reader(open('training2017/REFERENCE.csv')))
    for i, label in enumerate(testlabel):
      if label[1] == '~':
        filename = 'training2017/'+ label[0] + '.mat'
        from scipy.io import loadmat
        noise = loadmat(filename)
        noise = noise['val']
        _, size = noise.shape
        noise = noise.reshape(size,)
        noise = np.nan_to_num(noise) # removing NaNs and Infs
        from scipy.signal import resample
        noise= resample(noise, int(len(noise) * 360 / 300) ) # resample to match the data sampling rate 360(mit), 300(cinc)
        from sklearn import preprocessing
        noise = preprocessing.scale(noise)
        noise = noise/1000*6 # rough normalize, to be improved 
        from scipy.signal import find_peaks
        peaks, _ = find_peaks(noise, distance=150)
        choices = 10 # 256*10 from 9000
        picked_peaks = np.random.choice(peaks, choices, replace=False)
        for j, peak in enumerate(picked_peaks):
          if peak > config.input_size//2 and peak < len(noise) - config.input_size//2:
              start,end  = peak-config.input_size//2, peak+config.input_size//2
              if i > len(testlabel)/6:
                noises["trainset"].append(noise[start:end].tolist())
              else:
                noises["testset"].append(noise[start:end].tolist())
    return noises 

def get_picks_from_pdf(feature, phase, pick_set, height=0.5, distance=100):
    i = feature.phase.index(phase)
    peaks, properties = find_peaks(feature.pdf[-1, :, i], height=height, distance=distance)

    for p in peaks:
        if p:
            pick_time = UTCDateTime(feature.starttime) + p * feature.delta
            feature.pick_time.append(pick_time.isoformat())
            feature.pick_phase.append(feature.phase[i])
            feature.pick_set.append(pick_set) 

def find_peaks_simple(x):
    if x.ndim == 1:
        # Delegate to scipy's peak finder.
        return find_peaks(x)[0],
    else:
        # Use maximum filter for peak finding.
        y = maximum_filter(x, 3)
        return np.where(x == y) 

def PitchSpectralAcf(X, f_s):

    # initialize
    f_min = 300
    f = np.zeros(X.shape[1])

    # use spectral symmetry for robustness
    X[0, :] = np.max(X)
    X = np.concatenate((np.flipud(X), X), axis=0)

    # compute the ACF
    for n in range(0, X.shape[1]):

        if X[:, n].sum() < 1e-20:
            continue

        eta_min = int(round(f_min / f_s * (X.shape[0] - 2))) - 1

        afCorr = np.correlate(X[:, n], X[:, n], "full") / np.dot(X[:, n], X[:, n])
        afCorr = afCorr[np.arange(X.shape[0], afCorr.size)]

        # find the highest local maximum
        iPeaks = find_peaks(afCorr, height=0)
        if iPeaks[0].size:
            eta_min = np.max([eta_min, iPeaks[0][0] - 1])
        f[n] = np.argmax(afCorr[np.arange(eta_min, afCorr.size)]) + 1

        # find max index and convert to Hz (note: X has double length)
        f[n] = (f[n] + eta_min) / (X.shape[0] - 2) * f_s

    return (f) 

def panPeakDetect(detection, fs):    

    min_distance = int(0.25*fs)
    peaks, _ = signal.find_peaks(detection, distance=min_distance)      

    signal_peaks = []
    noise_peaks = []

    SPKI = 0.0
    NPKI = 0.0

    threshold_I1 = 0.0
    threshold_I2 = 0.0

    RR_missed = 0
    index = 0
    indexes = []

    missed_peaks = []
    for peak in peaks:

        if detection[peak] > threshold_I1:
               
            signal_peaks.append(peak)
            indexes.append(index)
            SPKI = 0.125*detection[signal_peaks[-1]] + 0.875*SPKI
            if RR_missed!=0:
                if signal_peaks[-1]-signal_peaks[-2]>RR_missed:
                    missed_section_peaks = peaks[indexes[-2]+1:indexes[-1]]
                    missed_section_peaks2 = []
                    for missed_peak in missed_section_peaks:
                        if missed_peak-signal_peaks[-2]>min_distance and signal_peaks[-1]-missed_peak>min_distance and detection[missed_peak]>threshold_I2:
                            missed_section_peaks2.append(missed_peak)

                    if len(missed_section_peaks2)>0:           
                        missed_peak = missed_section_peaks2[np.argmax(detection[missed_section_peaks2])]
                        missed_peaks.append(missed_peak)
                        signal_peaks.append(signal_peaks[-1])
                        signal_peaks[-2] = missed_peak   
        
        else:
            noise_peaks.append(peak)
            NPKI = 0.125*detection[noise_peaks[-1]] + 0.875*NPKI
        
        threshold_I1 = NPKI + 0.25*(SPKI-NPKI)
        threshold_I2 = 0.5*threshold_I1

        if len(signal_peaks)>8:
            RR = np.diff(signal_peaks[-9:])
            RR_ave = int(np.mean(RR))
            RR_missed = int(1.66*RR_ave)

        index = index+1

    signal_peaks.pop(0)

    return signal_peaks 

def get_standard_deviation_of_summed_textline_patch_along_width(self,img_patch,sigma_,multiplier=3.8 ):
        img_patch_sum_along_width=img_patch[:,:].sum(axis=1)

        img_patch_sum_along_width_updown=img_patch_sum_along_width[len(img_patch_sum_along_width)::-1]

        first_nonzero=(next((i for i, x in enumerate(img_patch_sum_along_width) if x), 0))
        last_nonzero=(next((i for i, x in enumerate(img_patch_sum_along_width_updown) if x), 0))

        last_nonzero=len(img_patch_sum_along_width)-last_nonzero


        y=img_patch_sum_along_width#[first_nonzero:last_nonzero]

        y_help=np.zeros(len(y)+20)

        y_help[10:len(y)+10]=y

        x=np.array( range(len(y)) )




        zneg_rev=-y_help+np.max(y_help)

        zneg=np.zeros(len(zneg_rev)+20)

        zneg[10:len(zneg_rev)+10]=zneg_rev

        z=gaussian_filter1d(y, sigma_)
        zneg= gaussian_filter1d(zneg, sigma_)


        peaks_neg, _ = find_peaks(zneg, height=0)
        peaks, _ = find_peaks(z, height=0)

        peaks_neg=peaks_neg-10-10
        
        interest_pos=z[peaks]
        
        interest_pos=interest_pos[interest_pos>10]
        
        interest_neg=z[peaks_neg]
        
        min_peaks_pos=np.mean(interest_pos)
        min_peaks_neg=0#np.min(interest_neg)
        
        dis_talaei=(min_peaks_pos-min_peaks_neg)/multiplier
        #print(interest_pos)
        grenze=min_peaks_pos-dis_talaei#np.mean(y[peaks_neg[0]:peaks_neg[len(peaks_neg)-1]])-np.std(y[peaks_neg[0]:peaks_neg[len(peaks_neg)-1]])/2.0

        interest_neg_fin=interest_neg[(interest_neg<grenze)]
        peaks_neg_fin=peaks_neg[(interest_neg<grenze)]
        interest_neg_fin=interest_neg[(interest_neg<grenze)]
        
        return interest_neg_fin,np.std(z) 

def reshape_4DSTEM_SumFrames(data):
    """
    Reshapes the lazy-imported stack of dimensions: (xxxxxx|Det_X, Det_Y) to the correct scan pattern
    shape: (x, y | Det_X, Det_Y) when the frame exposure times are not in headre bits.
    It utilises the over-exposed fly-back frame to identify the start of the lines in the first 20
    lines of frames,checks line length consistancy and finds the number of frames to skip at the
    beginning (this number is printed out as string output).

    Parameters
    ----------
    data : pyxem LazyElectronDiffraction2D
        imported mib file with diensions of: framenumbers|Det_X, Det_Y

    Returns
    -------
    data_reshaped : pyxem.signals.LazyElectronDiffraction2D
        reshaped data (x, y | Det_X, Det_Y)
    """
    # Assuming sacn_x, i.e. number of probe positions in a line is square root
    # of total number of frames
    scan_x = int(np.sqrt(data.axes_manager[0].size))
    # detector size in pixels
    det_size = data.axes_manager[1].size
    # crop the first ~20 lines
    data_crop = data.inav[0 : 20 * scan_x]
    data_crop_t = data_crop.T
    data_crop_t_sum = data_crop_t.sum()
    # summing over patterns
    intensity_array = data_crop_t_sum.data
    intensity_array = intensity_array.compute()
    peaks = find_peaks(intensity_array, distance=scan_x)
    # Difference between consecutive elements of the array
    lines = np.ediff1d(peaks[0])
    # Assuming the last element to be the line length
    line_len = lines[-1]
    if line_len != scan_x:
        raise ValueError("Fly_back does not correspond to correct line length!")
        #  Exits this routine and reshapes using scan size instead
    # Checking line lengths
    check = np.ravel(np.where(lines == line_len, True, False))
    # In case there is a False in there take the index of the last False
    if ~np.all(check):
        start_ind = np.where(check == False)[0][-1] + 2
        # Adding 2 - instead of 1 -  to be sure scan is OK
        skip_ind = peaks[0][start_ind]
    # In case they are all True take the index of the first True
    else:
        # number of frames to skip at the beginning
        skip_ind = peaks[0][0]
    # Number of lines
    n_lines = floor((data.data.shape[0] - skip_ind) / line_len)
    # with the skipped frames removed
    data_skip = data.inav[skip_ind : skip_ind + (n_lines * line_len)]
    data_skip.data = data_skip.data.reshape(n_lines, line_len, det_size, det_size)
    data_skip.axes_manager._axes.insert(0, data_skip.axes_manager[0].copy())
    data_skip.get_dimensions_from_data()
    print("Reshaping using the frame intensity sums of the first 20 lines")
    print("Number of frames skipped at the beginning: ", skip_ind)
    # Cropping the flyaback pixel at the start
    data_skip = data_skip.inav[1:]
    return data_skip 

def computeNoveltyFunction(cNoveltyName, afAudioData, f_s, afWindow=None, iBlockLength=4096, iHopLength=512):

    # compute window function for FFT
    if afWindow is None:
        afWindow = ToolComputeHann(iBlockLength)

    assert(afWindow.shape[0] == iBlockLength), "parameter error: invalid window dimension"

    #mypackage = __import__(".Novelty" + cNoveltyName, package="pyACA")
    hNoveltyFunc = getattr(pyACA, "Novelty" + cNoveltyName)

    # initialization
    fLengthLpInS = 0.3
    iLengthLp = np.max([2, math.ceil(fLengthLpInS * f_s / iHopLength)])

    # pre-processing
    afAudioData = ToolPreprocAudio(afAudioData, iBlockLength)

    # in the real world, we would do this block by block...
    [f, t, X] = spectrogram(afAudioData,
                            f_s,
                            afWindow,
                            iBlockLength,
                            iBlockLength - iHopLength,
                            iBlockLength,
                            False,
                            True,
                            'spectrum')

    #  scale the same as for matlab
    X = np.sqrt(X / 2)

    # novelty function
    d = hNoveltyFunc(X, f_s)

    # smooth novelty function
    b = np.ones(10) / 10
    d = filtfilt(b, 1, d)
    d[d < 0] = 0

    # compute threshold
    b = np.ones(iLengthLp) / iLengthLp
    G_T = .5 * np.mean(d[np.arange(1, d.shape[0])]) + filtfilt(b, 1, d)

    # find local maxima above the threshold
    iPeaks = find_peaks(d - G_T, height=0)

    return (d, t, iPeaks[0]) 

def test_k2_peaks_periodic(self):
        """Tests the correct peak locations and intensities are found for the
        k=2 term in periodic systems.
        """
        atoms = Atoms(
            cell=[
                [10, 0, 0],
                [10, 10, 0],
                [10, 0, 10],
            ],
            symbols=["H", "C"],
            scaled_positions=[
                [0.1, 0.5, 0.5],
                [0.9, 0.5, 0.5],
            ]
        )

        desc = LMBTR(
            species=["H", "C"],
            k2={
                "geometry": {"function": "distance"},
                "grid": {"min": 0, "max": 10, "sigma": 0.5, "n": 1000},
                "weighting": {"function": "exp", "scale": 0.8, "cutoff": 1e-3},
            },
            normalize_gaussians=False,
            periodic=True,
            flatten=True,
            sparse=False
        )
        features = desc.create(atoms, [0])[0, :]
        x = desc.get_k2_axis()

        # Calculate assumed locations and intensities.
        assumed_locs = np.array([2, 8])
        assumed_ints = np.exp(-0.8*np.array([2, 8]))

        # Check the X-C peaks
        xc_feat = features[desc.get_location(("X", "C"))]
        xc_peak_indices = find_peaks(xc_feat, prominence=0.001)[0]
        xc_peak_locs = x[xc_peak_indices]
        xc_peak_ints = xc_feat[xc_peak_indices]
        self.assertTrue(np.allclose(xc_peak_locs, assumed_locs, rtol=0, atol=1e-2))
        self.assertTrue(np.allclose(xc_peak_ints, assumed_ints, rtol=0, atol=1e-2))

        # Check that everything else is zero
        features[desc.get_location(("X", "C"))] = 0
        self.assertEqual(features.sum(), 0) 

def test_k3_peaks_finite(self):
        """Tests the correct peak locations and intensities are found for the
        k=3 term in finite systems.
        """
        desc = LMBTR(
            species=["H", "O"],
            k3={
                "geometry": {"function": "angle"},
                "grid": {"min": -10, "max": 180, "sigma": 5, "n": 2000},
                "weighting": {"function": "unity"},
            },
            normalize_gaussians=False,
            periodic=False,
            flatten=True,
            sparse=False
        )
        features = desc.create(H2O, [0])[0, :]
        x = desc.get_k3_axis()

        # Check the X-H-O peaks
        xho_assumed_locs = np.array([38])
        xho_assumed_ints = np.array([1])
        xho_feat = features[desc.get_location(("X", "H", "O"))]
        xho_peak_indices = find_peaks(xho_feat, prominence=0.5)[0]
        xho_peak_locs = x[xho_peak_indices]
        xho_peak_ints = xho_feat[xho_peak_indices]
        self.assertTrue(np.allclose(xho_peak_locs, xho_assumed_locs, rtol=0, atol=5e-2))
        self.assertTrue(np.allclose(xho_peak_ints, xho_assumed_ints, rtol=0, atol=5e-2))

        # Check the X-O-H peaks
        xoh_assumed_locs = np.array([104])
        xoh_assumed_ints = np.array([1])
        xoh_feat = features[desc.get_location(("X", "O", "H"))]
        xoh_peak_indices = find_peaks(xoh_feat, prominence=0.5)[0]
        xoh_peak_locs = x[xoh_peak_indices]
        xoh_peak_ints = xoh_feat[xoh_peak_indices]
        self.assertTrue(np.allclose(xoh_peak_locs, xoh_assumed_locs, rtol=0, atol=5e-2))
        self.assertTrue(np.allclose(xoh_peak_ints, xoh_assumed_ints, rtol=0, atol=5e-2))

        # Check the H-X-O peaks
        hxo_assumed_locs = np.array([38])
        hxo_assumed_ints = np.array([1])
        hxo_feat = features[desc.get_location(("H", "X", "O"))]
        hxo_peak_indices = find_peaks(hxo_feat, prominence=0.5)[0]
        hxo_peak_locs = x[hxo_peak_indices]
        hxo_peak_ints = hxo_feat[hxo_peak_indices]
        self.assertTrue(np.allclose(hxo_peak_locs, hxo_assumed_locs, rtol=0, atol=5e-2))
        self.assertTrue(np.allclose(hxo_peak_ints, hxo_assumed_ints, rtol=0, atol=5e-2))

        # Check that everything else is zero
        features[desc.get_location(("X", "H", "O"))] = 0
        features[desc.get_location(("X", "O", "H"))] = 0
        features[desc.get_location(("H", "X", "O"))] = 0
        self.assertEqual(features.sum(), 0) 

def PitchTimeAuditory(x, iBlockLength, iHopLength, f_s):

    # initialize
    iNumOfBlocks = math.ceil(x.size / iHopLength)
    f = np.zeros(iNumOfBlocks)
    f_max = 2000
    iNumBands = 20
    fLengthLpInS = 0.001

    iLengthLp = math.ceil(fLengthLpInS * f_s)

    # compute time stamps
    t = (np.arange(0, iNumOfBlocks) * iHopLength + (iBlockLength / 2)) / f_s

    # apply filterbank
    X = ToolGammatoneFb(x, f_s, iNumBands)

    # half wave rectification
    X[X < 0] = 0

    # smooth the results with a moving average filter
    b = np.ones(iLengthLp) / iLengthLp
    X = filtfilt(b, 1, X)

    for n in range(0, iNumOfBlocks):

        eta_min = int(round(f_s / f_max))
        afSumCorr = np.zeros(iBlockLength - 1)
        x_tmp = np.zeros(iBlockLength)

        i_start = n * iHopLength
        i_stop = np.min([x.size - 1, i_start + iBlockLength - 1])

        # compute ACF per band and summarize
        for k in range(0, iNumBands):
            # get current block
            if X[k, np.arange(i_start, i_stop + 1)].sum() < 1e-20:
                continue
            else:
                x_tmp[np.arange(0, i_stop - i_start + 1)] = X[k, np.arange(i_start, i_stop + 1)]

            afCorr = np.correlate(x_tmp, x_tmp, "full") / np.dot(x_tmp, x_tmp)

            # aggregate bands with simple sum before peak picking
            afSumCorr += afCorr[np.arange(iBlockLength, afCorr.size)]

        if afSumCorr.sum() < 1e-20:
            continue

        # find the highest local maximum
        iPeaks = find_peaks(afSumCorr, height=0)
        if iPeaks[0].size:
            eta_min = np.max([eta_min, iPeaks[0][0] - 1])
        f[n] = np.argmax(afSumCorr[np.arange(eta_min + 1, afSumCorr.size)]) + 1

        # convert to Hz
        f[n] = f_s / (f[n] + eta_min + 1)

    return (f, t) 

def test_k3_peaks_finite(self):
        """Tests that all the correct angles are present in finite systems.
        There should be n*(n-1)*(n-2)/2 unique angles where the division by two
        gets rid of duplicate angles.
        """
        desc = MBTR(
            species=["H", "O"],
            k3={
                "geometry": {"function": "angle"},
                "grid": {"min": -10, "max": 180, "sigma": 5, "n": 2000},
                "weighting": {"function": "unity"},
            },
            normalize_gaussians=False,
            periodic=False,
            flatten=True,
            sparse=False
        )
        features = desc.create(H2O)[0, :]
        x = desc.get_k3_axis()

        # Check the H-H-O peaks
        hho_assumed_locs = np.array([38])
        hho_assumed_ints = np.array([2])
        hho_feat = features[desc.get_location(("H", "H", "O"))]
        hho_peak_indices = find_peaks(hho_feat, prominence=0.5)[0]
        hho_peak_locs = x[hho_peak_indices]
        hho_peak_ints = hho_feat[hho_peak_indices]
        self.assertTrue(np.allclose(hho_peak_locs, hho_assumed_locs, rtol=0, atol=5e-2))
        self.assertTrue(np.allclose(hho_peak_ints, hho_assumed_ints, rtol=0, atol=5e-2))

        # Check the H-O-H peaks
        hoh_assumed_locs = np.array([104])
        hoh_assumed_ints = np.array([1])
        hoh_feat = features[desc.get_location(("H", "O", "H"))]
        hoh_peak_indices = find_peaks(hoh_feat, prominence=0.5)[0]
        hoh_peak_locs = x[hoh_peak_indices]
        hoh_peak_ints = hoh_feat[hoh_peak_indices]
        self.assertTrue(np.allclose(hoh_peak_locs, hoh_assumed_locs, rtol=0, atol=5e-2))
        self.assertTrue(np.allclose(hoh_peak_ints, hoh_assumed_ints, rtol=0, atol=5e-2))

        # Check that everything else is zero
        features[desc.get_location(("H", "H", "O"))] = 0
        features[desc.get_location(("H", "O", "H"))] = 0
        self.assertEqual(features.sum(), 0) 

def test_k3_peaks_periodic(self):
        """Tests that the final spectra does not change when translating atoms
        in a periodic cell. This is not trivially true unless the weight of
        angles is weighted according to the cell indices of the involved three
        atoms. Notice that the values of the geometry and weight functions are
        not equal before summing them up in the final graph.
        """
        scale = 0.85
        desc = MBTR(
            species=["H"],
            k3={
                "geometry": {"function": "angle"},
                "grid": {"min": 0, "max": 180, "sigma": 5, "n": 2000},
                "weighting": {"function": "exp", "scale": scale, "cutoff": 1e-3},
            },
            normalize_gaussians=False,
            periodic=True,
            flatten=True,
            sparse=False
        )

        atoms = Atoms(
            cell=[
                [10, 0, 0],
                [0, 10, 0],
                [0, 0, 10],
            ],
            symbols=3*["H"],
            scaled_positions=[
                [0.05, 0.40, 0.5],
                [0.05, 0.60, 0.5],
                [0.95, 0.5, 0.5],
            ],
            pbc=True
        )
        features = desc.create(atoms)[0, :]
        x = desc.get_k3_axis()

        # Calculate assumed locations and intensities.
        assumed_locs = np.array([45, 90])
        dist = 2+2*np.sqrt(2)  # The total distance around the three atoms
        weight = np.exp(-scale*dist)
        assumed_ints = np.array([4*weight, 2*weight])
        assumed_ints /= 2  # The periodic distances ar halved because they belong to different cells

        # Check the H-H-H peaks
        hhh_feat = features[desc.get_location(("H", "H", "H"))]
        hhh_peak_indices = find_peaks(hhh_feat, prominence=0.01)[0]
        hhh_peak_locs = x[hhh_peak_indices]
        hhh_peak_ints = hhh_feat[hhh_peak_indices]
        self.assertTrue(np.allclose(hhh_peak_locs, assumed_locs, rtol=0, atol=1e-1))
        self.assertTrue(np.allclose(hhh_peak_ints, assumed_ints, rtol=0, atol=1e-1))

        # Check that everything else is zero
        features[desc.get_location(("H", "H", "H"))] = 0
        self.assertEqual(features.sum(), 0)