Python scipy.signal.argrelextrema() Examples

def find_peak_vextors_eagerly(price, offest=0):
    """
    （饥渴的）在 MACD 上坡的时候查找更多的极值点
    """
    xn = price

    # pass 0
    window_size, poly_order = 5, 1
    yy_sg = savgol_filter(xn, window_size, poly_order)

    # pass 1
    x_tp_min, x_tp_max = signal.argrelextrema(yy_sg, np.less)[0], signal.argrelextrema(yy_sg, np.greater)[0]
    n = int(len(price) / (len(x_tp_min) + len(x_tp_max))) * 2

    # peakutils 似乎一根筋只能查最大极值，通过曲线反相的方式查找极小点
    mirrors = (yy_sg * -1) + np.mean(price) * 2

    # pass 2 使用 peakutils 查找
    x_tp_max = peakutils.indexes(yy_sg, thres=0.01 / max(price), min_dist=n)
    x_tp_min = peakutils.indexes(mirrors, thres=0.01 / max(price), min_dist=n)

    return x_tp_min + offest, x_tp_max + offest 

def ma30_cross(data):
        MA5 = talib.MA(data.close, 5)
        MA30 = talib.MA(data.close, 30)
    
        MA30_CROSS_JX = CROSS(MA5, MA30)
        MA30_CROSS_JX_Integral = Timeline_Integral_with_cross_before(MA30_CROSS_JX)
        MA30_CROSS_SX = CROSS(MA30, MA5)
        MA30_CROSS_SX_Integral = Timeline_Integral_with_cross_before(MA30_CROSS_SX)
    
        MA30_CROSS = pd.DataFrame(columns=['MA30_CROSS', 'MA30_CROSS_JX', 'MA30_CROSS_SX', 'MA30_TP_CROSS_JX', 'MA30_TP_CROSS_SX'], index=data.index)
        MA30_CROSS.loc[MA30_CROSS_JX == 1, 'MA30_CROSS'] = 1
        MA30_CROSS.loc[MA30_CROSS_SX == 1, 'MA30_CROSS'] = -1
        MA30_CROSS['MA30_CROSS_JX'] = Timeline_Integral_with_cross_before(MA30_CROSS_JX)
        MA30_CROSS['MA30_CROSS_SX'] = Timeline_Integral_with_cross_before(MA30_CROSS_SX)
    
        # MA30 前29个是 NaN，处理会抛出 Warning，使用 [29:] 则不会计算 NaN，相应的 return_index+29
        MA30_tp_min, MA30_tp_max = signal.argrelextrema(MA30.values[29:], np.less)[0] + 29, signal.argrelextrema(MA30.values[29:], np.greater)[0] + 29
        MA30_TP_CROSS = pd.DataFrame(columns=['MA30_TP_CROSS_JX', 'MA30_TP_CROSS_SX'], index=data.index)
        MA30_TP_CROSS['MA30_TP_CROSS_SX'] = MA30_TP_CROSS['MA30_TP_CROSS_JX'] = 0
        MA30_TP_CROSS.iloc[MA30_tp_min, MA30_TP_CROSS.columns.get_loc('MA30_TP_CROSS_JX')] = 1
        MA30_TP_CROSS.iloc[MA30_tp_max, MA30_TP_CROSS.columns.get_loc('MA30_TP_CROSS_SX')] = 1
        MA30_CROSS['MA30_TP_CROSS_JX'] = Timeline_Integral_with_cross_before(MA30_TP_CROSS['MA30_TP_CROSS_JX'])
        MA30_CROSS['MA30_TP_CROSS_SX'] = Timeline_Integral_with_cross_before(MA30_TP_CROSS['MA30_TP_CROSS_SX'])
        return MA30_CROSS 

def find_ori_ter(c_skew, length):
    """
    find origin and terminus of replication based on 
    cumulative GC Skew
    """
    # find origin and terminus of replication based on 
    # cumulative gc skew min and max peaks
    c_skew_min = signal.argrelextrema(np.asarray(c_skew[1]), np.less, order = 1)[0].tolist()
    c_skew_max = signal.argrelextrema(np.asarray(c_skew[1]), np.greater, order = 1)[0].tolist()
    # return False if no peaks were detected
    if len(c_skew_min) == 0 or len(c_skew_min) == 0:
        return [False, False]
    else:
        c_skew_min = [[c_skew[0][i], c_skew[1][i]] for i in c_skew_min]
        c_skew_max = [[c_skew[0][i], c_skew[1][i]] for i in c_skew_max]
        ori, ter = check_peaks([c_skew_min, c_skew_max], length)
    return ori, ter 

def process_merlin_label(bin_label_fname, text_lab_dir, phonedim=416, subphonedim=9):

    text_label = os.path.join(text_lab_dir, basename(bin_label_fname) + '.lab')
    assert os.path.isfile(text_label), 'No text file for %s '%(basename(bin_label_fname))
    
    labfrombin = get_speech(bin_label_fname, phonedim+subphonedim)
    
    ## fraction through phone (forwards)
    fraction_through_phone_forwards = labfrombin[:,-1]

    ## This is a suprisingly noisy signal which never seems to start at 0.0! Find minima:-
    (minima, ) = argrelextrema(fraction_through_phone_forwards, np.less)

    ## first frame is always a start: 
    minima = np.insert(minima, 0, 0)  

    ## check size against text file:
    labfromtext = merlin_state_label_to_phone(text_label)
    assert labfromtext.shape[0] == minima.shape[0]

    lab = labfrombin[minima,:-subphonedim] ## discard frame level feats, and take first frame of each phone

    return lab 

def peaks(spectra,frequency,number=3,thresh=0.01):
        """ Return the peaks from the Fourier transform
            Variables:
            number:     integer. number of peaks to print.
            thresh:     float. Threshhold intensity for printing.

            Returns: Energy (eV), Intensity (depends on type of spectra)
        """

        from scipy.signal import argrelextrema as pks
        # find all peak indices [idx], and remove those below thresh [jdx]
        idx = pks(np.abs(spectra),np.greater,order=3)
        jdx = np.where((np.abs(spectra[idx]) >= thresh))
        kdx = idx[0][jdx[0]] # indices of peaks matching criteria
        if number > len(kdx):
            number = len(kdx)
        print("First "+str(number)+" peaks (eV) found: ")
        for i in xrange(number):
            print("{0:.4f}".format(frequency[kdx][i]*27.2114),
                  "{0:.4f}".format(spectra[kdx][i])) 

def test_find_extrema(first_extrema):
    """Test ability to find peaks and troughs."""

    # Load signal
    sig = np.load(DATA_PATH + 'sim_stationary.npy')

    fs = 1000
    f_range = (6, 14)

    # find local maxima and minima using scipy
    maxima = argrelextrema(sig, np.greater)
    minima = argrelextrema(sig, np.less)

    # Find peaks and troughs using bycycle and make sure match scipy
    ps, ts = find_extrema(sig, fs, f_range, boundary=1, first_extrema=first_extrema)

    if first_extrema == 'trough':
        assert len(ps) == len(ts)
        assert ts[0] < ps[0]
        np.testing.assert_allclose(ps, maxima[0])
        np.testing.assert_allclose(ts[:len(ps)], minima[0][:len(ps)])
    elif first_extrema == 'peak':
        assert ps[0] < ts[0] 

def arg_peaks(data, extrema_chk = extrema_chk):
    """
    the function search peaks of spectrum and return positions of all peaks
    if extrema_chk == 1 uses numpy module
    if extrema_chk == 0 uses independent code (see below)
    """
    if extrema_chk == 1:
        return argrelextrema(data, np.greater)[0]
    else:
        diff_y = np.diff(data)
        extrm_y = np.diff(np.sign(diff_y))
        return np.where(extrm_y<0)[0]+1 

def find_extreme(self, obs):
        max_idx = argrelextrema(obs.open.values, np.greater, order=self.peak_order)[0]
        min_idx = argrelextrema(obs.open.values, np.less, order=self.peak_order)[0]
        extreme_idx = np.concatenate([max_idx, min_idx, [obs.shape[0] - 1]])
        extreme_idx.sort()
        return obs.open.iloc[extreme_idx] 

def _compute_peaks_or_valleys_of_first_derivative(s, do_peaks=True):
    """
    Takes a spectrogram and returns a 2D array of the form:

    0 0 0 1 0 0 1 0 0 0 1   <-- Frequency 0
    0 0 1 0 0 0 0 0 0 1 0   <-- Frequency 1
    0 0 0 0 0 0 1 0 1 0 0   <-- Frequency 2
    *** Time axis *******

    Where a 1 means that the value in that time bin in the spectrogram corresponds to
    a peak/valley in the first derivative.

    This function is used as part of the ASA algorithm and is not meant to be used publicly.
    """
    # Get the first derivative of each frequency in the time domain
    gradient = np.nan_to_num(np.apply_along_axis(np.gradient, 1, s), copy=False)

    # Calculate the value we will use for determinig whether something is an event or not
    threshold = np.squeeze(np.nanmean(gradient, axis=1) + np.nanstd(gradient, axis=1))

    # Look for relative extrema along the time dimension
    half_window = 4
    if do_peaks:
        indexes = [signal.argrelextrema(gradient[i, :], np.greater, order=half_window)[0] for i in range(gradient.shape[0])]
    else:
        indexes = [signal.argrelextrema(gradient[i, :], np.less, order=half_window)[0] for i in range(gradient.shape[0])]

    # indexes should now contain the indexes of possible extrema
    # But we need to filter out values that are not large enough, and we want the end result
    # to be a 1 or 0 mask corresponding to locations of extrema
    extrema = np.zeros(s.shape)
    for row_index, index_array in enumerate(indexes):
        # Each index_array is a list of indexes corresponding to all the extrema in a given row
        for col_index in index_array:
            if do_peaks and (gradient[row_index, col_index] > threshold[row_index]):
                extrema[row_index, col_index] = 1
            elif not do_peaks:
                # Note that we do not remove under-threshold values from the offsets - these will be taken care of later in the algo
                extrema[row_index, col_index] = 1
    return extrema, gradient 

def _compute_peaks_or_valleys_of_first_derivative(s, do_peaks=True):
    """
    Takes a spectrogram and returns a 2D array of the form:

    0 0 0 1 0 0 1 0 0 0 1   <-- Frequency 0
    0 0 1 0 0 0 0 0 0 1 0   <-- Frequency 1
    0 0 0 0 0 0 1 0 1 0 0   <-- Frequency 2
    *** Time axis *******

    Where a 1 means that the value in that time bin in the spectrogram corresponds to
    a peak/valley in the first derivative.

    This function is used as part of the ASA algorithm and is not meant to be used publicly.
    """
    # Get the first derivative of each frequency in the time domain
    gradient = np.nan_to_num(np.apply_along_axis(np.gradient, 1, s), copy=False)

    # Calculate the value we will use for determinig whether something is an event or not
    threshold = np.squeeze(np.nanmean(gradient, axis=1) + np.nanstd(gradient, axis=1))

    # Look for relative extrema along the time dimension
    half_window = 4
    if do_peaks:
        indexes = [signal.argrelextrema(gradient[i, :], np.greater, order=half_window)[0] for i in range(gradient.shape[0])]
    else:
        indexes = [signal.argrelextrema(gradient[i, :], np.less, order=half_window)[0] for i in range(gradient.shape[0])]

    # indexes should now contain the indexes of possible extrema
    # But we need to filter out values that are not large enough, and we want the end result
    # to be a 1 or 0 mask corresponding to locations of extrema
    extrema = np.zeros(s.shape)
    for row_index, index_array in enumerate(indexes):
        # Each index_array is a list of indexes corresponding to all the extrema in a given row
        for col_index in index_array:
            if do_peaks and (gradient[row_index, col_index] > threshold[row_index]):
                extrema[row_index, col_index] = 1
            elif not do_peaks:
                # Note that we do not remove under-threshold values from the offsets - these will be taken care of later in the algo
                extrema[row_index, col_index] = 1
    return extrema, gradient 

def minmax(candles: np.ndarray, order=3, sequential=False) -> EXTREMA:
    """
    minmax - Get extrema

    :param candles: np.ndarray
    :param order: int - default = 3
    :param sequential: bool - default=False

    :return: EXTREMA(min, max, last_min, last_max)
    """
    if not sequential and len(candles) > 240:
        candles = candles[-240:]

    low = candles[:, 4]
    high = candles[:, 3]

    minimaIdxs = argrelextrema(low, np.less, order=order, axis=0)
    maximaIdxs = argrelextrema(high, np.greater, order=order, axis=0)

    min = np.full_like(low, np.nan)
    max = np.full_like(high, np.nan)

    # set the extremas with the matching price
    min[minimaIdxs] = low[minimaIdxs]
    max[maximaIdxs] = high[maximaIdxs]

    # forward fill Nan values to get the last extrema
    last_min = np_ffill(min)
    last_max = np_ffill(max)

    if sequential:
        return EXTREMA(min, max, last_min, last_max)
    else:
        return EXTREMA(min[-1], max[-1], last_min[-1], last_max[-1]) 

def find_peak_vextors(price, return_ref=False, offest=0):
    """
    采用巴特沃斯信号滤波器，自适应寻找最佳极值点，决定平均周期分段数量。
    使用 scipy.Gaussian 机器学习统计算法进行第二次分析
    If you meet a Warning message, To slove this need upgrade scipy=>1.2. 
    but QUANTAXIS incompatible scipy=>1.2

    Parameters
    ----------
    price : (N,) array_like
        传入需要查找极值点的价格-时间序列。
        The numerator coefficient vector of the filter.
    return_ref : bool or None, optional
        返回作为参照的平滑曲线，平滑曲线的目的是减少锯齿抖动，减少计算的极值点。
        Return the smoothed line for reference.
    offest : int or None, optional
        传递参数时可能会被 .dropna() 或者 price[t:0] 等切片手段移除 nparray 头部
        的 np.nan 元素，因为此函数返回的是向量节点的数组索引，为了跟原始参数对应，调用者
        可以指定一个补偿偏移量，在返回的最大最小值中每个索引都会追加这个偏移量。
        The number of elements index offest, for jump np.nan in price's head.

    Returns
    -------
    x_tp_min, x_tp_max : ndarray
        包含最大值和最少值索引的数组
        The min/max peakpoint's index in array.

    """
    xn = price

    # Create an order 3 lowpass butterworth filter.
    b, a = butter(3, 0.05)

    # Apply the filter to xn.  Use lfilter_zi to choose the initial condition
    # of the filter.
    zi = lfilter_zi(b, a)
    z, _ = lfilter(b, a, xn, zi=zi * xn[0])

    # Apply the filter again, to have a result filtered at an order
    # the same as filtfilt.
    z2, _ = lfilter(b, a, z, zi=zi * z[0])

    # Use filtfilt to apply the filter.  If you meet a Warning need upgrade to
    # scipy=>1.2 but QUANTAXIS incompatible scipy=>1.2
    y = filtfilt(b, a, xn)

    # pass 1
    x_tp_min, x_tp_max = signal.argrelextrema(y, np.less)[0], signal.argrelextrema(y, np.greater)[0]
    n = int(len(price) / (len(x_tp_min) + len(x_tp_max))) * 2

    # peakutils 似乎一根筋只能查最大极值，通过曲线反相的方式查找极小点
    mirrors = (price * -1) + np.mean(price) * 2

    # pass 2 使用 peakutils 查找
    x_tp_max = peakutils.indexes(price, thres=0.01 / max(price), min_dist=n)
    x_tp_min = peakutils.indexes(mirrors, thres=0.01 / max(price), min_dist=n)

    if (return_ref):
        return x_tp_min + offest, x_tp_max + offest, y
    else:
        return x_tp_min + offest, x_tp_max + offest 

def calc_1D_coherent_fraction(U, axisName, axis, p=0):
    """
    Calculates 1D degree of coherence (DoC). From its width in respect to the
    width of intensity distribution also infers the coherent fraction. Both
    widths are rms. The one of intensity is calculated over the whole axis, the
    one of DoC is calculated between the first local minima around the center
    provided that these minima are lower than 0.5.

    *U*: complex valued ndarray, shape(repeats, nx, ny)
        3D stack of field images. For a 1D cut along *axis*, the middle of the
        other dimension is sliced out.

    *axis*: str, one of 'x' or ('y' or 'z')
        Specifies the 1D axis of interest.

    *p*: float, distance to screen
        If non-zero, the calculated mutual intensity will be divided by *p*\ ².
        This is useful to get proper physical units of the returned intensity
        if the function is applied directly to the field stacks given by
        Undulator.multi_electron_stack() that is calculated in angular units.

    Returns a tuple of mutual intensity, 1D intensity, 1D DoC, rms width of
    intensity, rms width of DoC (between the local minima, see above), the
    position of the minima (only the positive side) and the coherent fraction.
    This tuple can be fed to the plotting function
    :func:`plot_1D_degree_of_coherence`.

    """
    repeats, binsx, binsz = U.shape
    if axisName == 'x':
        Uc = U[:, :, binsz//2]
    elif axisName in ['y', 'z']:
        Uc = U[:, binsx//2, :]
    else:
        raise ValueError("unknown axis")
    J = np.dot(Uc.T.conjugate(), Uc) / repeats
    if p > 0:  # p is distance
        J /= p**2
    II = np.abs(np.diag(J))  # intensity as the main diagonal
#    II = (np.abs(Uc)**2).sum(axis=0) / repeats  # intensity by definition
    J /= II**0.5 * II[:, np.newaxis]**0.5
    Jd = np.abs(np.diag(np.fliplr(J)))  # DoC as the cross-diagonal

    varI = (II * axis**2).sum() / II.sum()
    axisEx = 2 * axis

    lm = argrelextrema(Jd, np.less)[0]  # local min
    # for variance up to the 1st local minimum which is < 0.5:
    lm = lm[(axisEx[lm] > 0) & (Jd[lm] < 0.5)]
    if len(lm) > 0:
        cond = np.abs(axisEx) <= axisEx[lm[0]]
        limJd = axisEx[lm[0]]
    else:
        cond = slice(None)  # for unconstrained variance calculation
        limJd = None
    varJd = (Jd * axisEx**2)[cond].sum() / Jd[cond].sum()
    cohFr = (4*varI/varJd + 1)**(-0.5)
    return J, II, Jd, varI, varJd, limJd, cohFr 

def __init__(self, x, y, S=1.0, curve='concave', direction='increasing'):
        """
        Once instantiated, this class attempts to find the point of maximum
        curvature on a line. The knee is accessible via the `.knee` attribute.
        :param x: x values.
        :type x: list or array.
        :param y: y values.
        :type y: list or array.
        :param S: Sensitivity, original paper suggests default of 1.0
        :type S: float
        :param curve: If 'concave', algorithm will detect knees. If 'convex', it
            will detect elbows.
        :type curve: string
        :param direction: one of {"increasing", "decreasing"}
        :type direction: string
        """
        # Step 0: Raw Input
        self.x = x
        self.y = y
        self.curve = curve
        self.direction = direction
        self.N = len(self.x)
        self.S = S

        # Step 1: fit a smooth line
        uspline = interpolate.interp1d(self.x, self.y)
        self.Ds_x = np.linspace(np.min(self.x), np.max(self.x), self.N)
        self.Ds_y = uspline(self.Ds_x)

        # Step 2: normalize values
        self.xsn = self.__normalize(self.Ds_x)
        self.ysn = self.__normalize(self.Ds_y)

        # Step 3: Calculate difference curve
        self.xd = self.xsn
        if self.curve == 'convex' and direction == 'decreasing':
            self.yd = self.ysn + self.xsn
            self.yd = 1 - self.yd
        elif self.curve == 'concave' and direction == 'decreasing':
            self.yd = self.ysn + self.xsn
        elif self.curve == 'concave' and direction == 'increasing':
            self.yd = self.ysn - self.xsn
        if self.curve == 'convex' and direction == 'increasing':
            self.yd = abs(self.ysn - self.xsn)

        # Step 4: Identify local maxima/minima
        # local maxima
        self.xmx_idx = argrelextrema(self.yd, np.greater)[0]
        self.xmx = self.xd[self.xmx_idx]
        self.ymx = self.yd[self.xmx_idx]

        # local minima
        self.xmn_idx = argrelextrema(self.yd, np.less)[0]
        self.xmn = self.xd[self.xmn_idx]
        self.ymn = self.yd[self.xmn_idx]

        # Step 5: Calculate thresholds
        self.Tmx = self.__threshold(self.ymx)

        # Step 6: find knee
        self.knee, self.norm_knee, self.knee_x = self.find_knee() 

def __init__(self, x, y, S=1.0, curve='concave', direction='increasing'):
        """
        Once instantiated, this class attempts to find the point of maximum
        curvature on a line. The knee is accessible via the `.knee` attribute.
        :param x: x values.
        :type x: list or array.
        :param y: y values.
        :type y: list or array.
        :param S: Sensitivity, original paper suggests default of 1.0
        :type S: float
        :param curve: If 'concave', algorithm will detect knees. If 'convex', it
            will detect elbows.
        :type curve: string
        :param direction: one of {"increasing", "decreasing"}
        :type direction: string
        """
        # Step 0: Raw Input
        self.x = x
        self.y = y
        self.curve = curve
        self.direction = direction
        self.N = len(self.x)
        self.S = S

        # Step 1: fit a smooth line
        uspline = interpolate.interp1d(self.x, self.y)
        self.Ds_x = np.linspace(np.min(self.x), np.max(self.x), self.N)
        self.Ds_y = uspline(self.Ds_x)

        # Step 2: normalize values
        self.xsn = self.__normalize(self.Ds_x)
        self.ysn = self.__normalize(self.Ds_y)

        # Step 3: Calculate difference curve
        self.xd = self.xsn
        if self.curve == 'convex' and direction == 'decreasing':
            self.yd = self.ysn + self.xsn
            self.yd = 1 - self.yd
        elif self.curve == 'concave' and direction == 'decreasing':
            self.yd = self.ysn + self.xsn
        elif self.curve == 'concave' and direction == 'increasing':
            self.yd = self.ysn - self.xsn
        if self.curve == 'convex' and direction == 'increasing':
            self.yd = abs(self.ysn - self.xsn)

        # Step 4: Identify local maxima/minima
        # local maxima
        self.xmx_idx = argrelextrema(self.yd, np.greater)[0]
        self.xmx = self.xd[self.xmx_idx]
        self.ymx = self.yd[self.xmx_idx]

        # local minima
        self.xmn_idx = argrelextrema(self.yd, np.less)[0]
        self.xmn = self.xd[self.xmn_idx]
        self.ymn = self.yd[self.xmn_idx]

        # Step 5: Calculate thresholds
        self.Tmx = self.__threshold(self.ymx)

        # Step 6: find knee
        self.knee, self.norm_knee, self.knee_x = self.find_knee() 

def compute_th(PSD, barks, ATH, freqs):
    """ returns the global masking threshold
    """
    # Identification of tonal maskers
    # find the index of maskers that are the local maxima
    length = len(PSD)
    masker_index = signal.argrelextrema(PSD, np.greater)[0]
    
    
    # delete the boundary of maskers for smoothing
    if 0 in masker_index:
        masker_index = np.delete(0)
    if length - 1 in masker_index:
        masker_index = np.delete(length - 1)
    num_local_max = len(masker_index)

    # treat all the maskers as tonal (conservative way)
    # smooth the PSD 
    p_k = pow(10, PSD[masker_index]/10.)    
    p_k_prev = pow(10, PSD[masker_index - 1]/10.)
    p_k_post = pow(10, PSD[masker_index + 1]/10.)
    P_TM = 10 * np.log10(p_k_prev + p_k + p_k_post)
    
    # bark_psd: the first column bark, the second column: P_TM, the third column: the index of points
    _BARK = 0
    _PSD = 1
    _INDEX = 2
    bark_psd = np.zeros([num_local_max, 3])
    bark_psd[:, _BARK] = barks[masker_index]
    bark_psd[:, _PSD] = P_TM
    bark_psd[:, _INDEX] = masker_index
    
    # delete the masker that doesn't have the highest PSD within 0.5 Bark around its frequency 
    for i in range(num_local_max):
        next = i + 1
        if next >= bark_psd.shape[0]:
            break
            
        while bark_psd[next, _BARK] - bark_psd[i, _BARK]  < 0.5:
            # masker must be higher than quiet threshold
            if quiet(freqs[int(bark_psd[i, _INDEX])]) > bark_psd[i, _PSD]:
                bark_psd = np.delete(bark_psd, (i), axis=0)
            if next == bark_psd.shape[0]:
                break
                
            if bark_psd[i, _PSD] < bark_psd[next, _PSD]:
                bark_psd = np.delete(bark_psd, (i), axis=0)
            else:
                bark_psd = np.delete(bark_psd, (next), axis=0)
            if next == bark_psd.shape[0]:
                break        
    
    # compute the individual masking threshold
    delta_TM = 1 * (-6.025  -0.275 * bark_psd[:, 0])
    Ts = two_slops(bark_psd, delta_TM, barks) 
    Ts = np.array(Ts)
    
    # compute the global masking threshold
    theta_x = np.sum(pow(10, Ts/10.), axis=0) + pow(10, ATH/10.) 
 
    return theta_x