--- a
+++ b/ppfunctions_1.py
@@ -0,0 +1,742 @@
+# -*- coding: utf-8 -*-
+"""
+Personal Processing Functions 1 (ppfunctions_1)
+Created on Mon Apr  9 11:48:37 2018
+@author: Kevin MAchado
+
+Module file containing Python definitions and statements
+"""
+
+# Libraries
+import numpy as np
+from scipy.signal import kaiserord, lfilter, firwin, butter
+from scipy.fftpack import fft
+from scipy import signal as sg
+#import peakutils                                # Librery to help in peak detection
+
+# Functions 
+# Normal energy = np.sum(pcgFFT1**2)
+# Normalized average Shannon Energy = sum((l**2)*np.log(l**2))/l.shape[0]
+# Third Order Shannon Energy = sum((l**3)*np.log(l**3))/l.shape[0]
+# -----------------------------------------------------------------------------
+#                           Statistic Variables
+# -----------------------------------------------------------------------------
+def _variance(x):
+    miu = 0.0
+    vari = 0.0
+    for i in range(x.shape[0]):
+        miu =+ x[i]
+    miu = miu/len(x)
+        
+    for i in range(x.shape[0]):
+        vari = vari + (x[i]-miu)**2
+    vari = vari/len(x)
+    return vari
+
+def _StandarD(x):
+    miu = 0.0
+    vari = 0.0
+    for i in range(x.shape[0]):
+        miu += x[i]
+    miu = miu/len(x)
+    for i in range(x.shape[0]):
+        vari = vari + (x[i]-miu)**2
+    vari = vari/len(x)
+    _Std = vari **(.5)
+    return _Std
+
+def _CV(x):
+    # Coefficient of Variation
+    var = _variance(x)
+    std = _StandarD(x)
+    
+    return 100*(var/std)
+
+# -----------------------------------------------------------------------------
+# PDS
+# -----------------------------------------------------------------------------
+def vec_nor(x):
+    """
+    Normalize the amplitude of a vector from -1 to 1
+    """
+    nVec=np.zeros(len(x));		   # Initializate derivate vector
+    nVec = np.divide(x, max(x))
+    nVec=nVec-np.mean(nVec);
+    nVec=np.divide(nVec,np.max(nVec));
+        
+    return nVec
+
+def running_sum(x):
+    """
+    Running Sum Algorithm of an input signal is y[n]= x[n] + y[n-1] 
+    """
+    y = np.zeros(len(x))
+    for i in range(len(x)):
+        y[i] = x[i] + y[i-1]
+        
+    return vec_nor(y)
+
+def derivate_1 (x):
+    """
+    Derivate of an input signal as y[n]= x[n] - x[n-1] 
+    """
+    y=np.zeros(len(x));				# Initializate derivate vector
+    for i in range(len(x)):
+        y[i]=x[i]-x[i-1];		
+    return vec_nor(y)
+        
+def derivate (x):
+    """
+    Derivate of an input signal as y[n]= x[n+1]- x[n-1] 
+    """
+    lenght=x.shape[0]				# Get the length of the vector
+    y=np.zeros(lenght);				# Initializate derivate vector
+    for i in range(lenght-1):
+        y[i]=x[i-1]-x[i];		
+    return y
+
+def derivate_positive (x):
+    """
+    Derivate of an input signal as y[n]= x[n+1]- x[n-1] 
+    for all values where the signal is positive
+    """
+    lenght=x.shape[0]				# Get the length of the vector
+    y=np.zeros(lenght);				# Initializate derivate vector
+    for i in range(lenght-1):
+        if x[i]>0:
+            y[i]=x[i-1]-x[i];		
+    return y
+# -----------------------------------------------------------------------------
+# Energy
+# -----------------------------------------------------------------------------
+def Energy_value (x):
+    """
+    Energy of an input signal  
+    """
+    y = np.sum(x**2)
+    return y
+
+def shannonE_value (x):
+    """
+    Shannon energy of an input signal  
+    """
+    y = sum((x**2)*np.log(x**2))/x.shape[0]
+    return y
+
+def shannonE_vector (x):
+    """
+    Shannon energy of an input signal  
+    """
+    mu = -(x**2)*np.log(x**2)/x.shape[0]
+    y = -(((x**2)*np.log(x**2)) - mu)/np.std(x)
+    return y
+
+def shannonE_vector_1 (x):
+    """
+    Shannon energy of an input signal  
+    """
+    N = x.shape[0]
+#    Se = -(1/N) * 
+    mu = -(x**2)*np.log(x**2)/x.shape[0]
+    y = -(((x**2)*np.log(x**2)) - mu)/np.std(x)
+    return y
+
+def E_VS_LF (pcgFFT1, vTfft1, on):
+    """
+    Energy of PCG Vibratory Spectrum in low Frequencies (E_VS_LF)
+    (frequency components, frequency value vector, on = on percentage or not)
+    According with [1] The total vibratory spectrum can be divided into 7 bands.
+    This is a modification of this 7 bands
+    1. 0-5Hz, 2. 5-25Hz; 3. 25-60Hz; 4. 60-120Hz; 5. 120-400Hz
+
+The PCG signal producess vibrations in the spectrum between 0-2k Hz. 
+[1] Abbas, Abbas K. (Abbas Khudair), Bassam, Rasha and Morgan & Claypool Publishers Phonocardiography signal processing. Morgan & Claypool Publishers, San Rafael, Calif, 2009.
+    """
+    c1 = (np.abs(vTfft1-5)).argmin()
+    c2 = (np.abs(vTfft1-25)).argmin()
+    c3 = (np.abs(vTfft1-120)).argmin()
+    c4 = (np.abs(vTfft1-240)).argmin()
+    c5 = (np.abs(vTfft1-500)).argmin()
+    
+    # All vector energy
+    xAll = Energy_value(pcgFFT1)
+
+    # Procesando de 0.01-5 Hz
+    pcgFFT_F1 = pcgFFT1[0:c1]
+    x1 = Energy_value(pcgFFT_F1)
+    
+    # Procesando de 5-25 Hz
+    pcgFFT_F2 = pcgFFT1[c1:c2]
+    x2 = Energy_value(pcgFFT_F2)
+    
+    # Procesando de 25-120 Hz
+    pcgFFT_F3 = pcgFFT1[c2:c3]
+    x3 = Energy_value(pcgFFT_F3)
+    
+    # Procesando de 120-240 Hz
+    pcgFFT_F4 = pcgFFT1[c3:c4]
+    x4 = Energy_value(pcgFFT_F4)
+    
+    # Procesando de 240-500 Hz
+    pcgFFT_F5 = pcgFFT1[c4:c5]
+    x5 = Energy_value(pcgFFT_F5)
+    
+    x = np.array([xAll, x1, x2, x3, x4, x5])
+    
+    if (on == 'percentage'):
+        x = 100*(x/x[0])
+
+    return x
+
+def E_VS (pcgFFT1, vTfft1, on):
+    """
+    Energy of PCG Vibratory Spectrum
+    (frequency components, frequency value vector, on = on percentage or not)
+    According with [1] The total vibratory spectrum can be divided into 7 bands:
+    1. 0-5Hz, 2. 5-25Hz; 3. 25-120Hz; 4. 120-240Hz; 5. 240-500Hz; 6. 500-1000Hz; 7. 1000-2000Hz
+
+The PCG signal producess vibrations in the spectrum between 0-2k Hz. 
+[1] Abbas, Abbas K. (Abbas Khudair), Bassam, Rasha and Morgan & Claypool Publishers Phonocardiography signal processing. Morgan & Claypool Publishers, San Rafael, Calif, 2009.
+    """
+    c1 = (np.abs(vTfft1-5)).argmin()
+    c2 = (np.abs(vTfft1-25)).argmin()
+    c3 = (np.abs(vTfft1-120)).argmin()
+    c4 = (np.abs(vTfft1-240)).argmin()
+    c5 = (np.abs(vTfft1-500)).argmin()
+    c6 = (np.abs(vTfft1-1000)).argmin()
+    c7 = (np.abs(vTfft1-2000)).argmin()
+    
+    # All vector energy
+    xAll = Energy_value(pcgFFT1)
+
+    # Procesando de 0.01-5 Hz
+    pcgFFT_F1 = pcgFFT1[0:c1]
+    x1 = Energy_value(pcgFFT_F1)
+    
+    # Procesando de 5-25 Hz
+    pcgFFT_F2 = pcgFFT1[c1:c2]
+    x2 = Energy_value(pcgFFT_F2)
+    
+    # Procesando de 25-120 Hz
+    pcgFFT_F3 = pcgFFT1[c2:c3]
+    x3 = Energy_value(pcgFFT_F3)
+    
+    # Procesando de 120-240 Hz
+    pcgFFT_F4 = pcgFFT1[c3:c4]
+    x4 = Energy_value(pcgFFT_F4)
+    
+    # Procesando de 240-500 Hz
+    pcgFFT_F5 = pcgFFT1[c4:c5]
+    x5 = Energy_value(pcgFFT_F5)
+    
+    # Procesando de 500-1000 Hz
+    pcgFFT_F6 = pcgFFT1[c5:c6]
+    x6 = Energy_value(pcgFFT_F6)
+    
+    # Procesando de 1000-2000 Hz
+    pcgFFT_F7 = pcgFFT1[c6:c7]
+    x7 = Energy_value(pcgFFT_F7)
+    
+    x = np.array([xAll, x1, x2, x3, x4, x5, x6, x7])
+    
+    if (on == 'percentage'):
+        x = 100*(x/x[0])
+
+    return x
+#-------------------------------------------
+def features_1(data_P1, fs):
+    
+    # Defining the Vibratory Frequency Bands
+    bVec = [0.01, 120, 240, 350, 425, 500, 999]
+    # Initializing Vectors
+    band_matrix = np.zeros((len(bVec),len(data_P1))) # Band Matrix
+    power_matrix = np.zeros((len(bVec),int(1+(len(data_P1)/2)))) # Band Matrix
+    freqs_matrix = np.zeros((len(bVec),int(1+(len(data_P1)/2)))) # Band Matrix
+    SePCG = np.zeros(len(bVec)-1)                            # Shannon Energy in each Band
+    SePWR = np.zeros(len(bVec)-1)                            # Shannon Energy in each Band
+      
+    
+    for i in range(len(bVec)-1):
+        band_matrix[i,:] = butter_bp_fil(data_P1, bVec[i], bVec[i+1], fs)       
+        freqs_matrix[i,:], power_matrix[i,:] = sg.periodogram(band_matrix[i,:], fs, scaling = 'density')
+        SePCG[i] = Energy_value(band_matrix[i,:])
+        SePWR[i] = Energy_value(power_matrix[i,:])
+     
+    SePWR = 1*np.log10(SePWR)
+    SePCG = 1*np.log10(SePCG)
+    
+    return SePWR, SePCG
+
+# -----------------------------------------------------------------------------
+# Filter Processes
+# -----------------------------------------------------------------------------
+    
+def recursive_moving_average_F(X, Fs, M):
+    """
+    The recursive Moving Average Filter its an algorithm to implement the typical
+    moving average filter more faster. The algorithm is written as:
+    y[i] = y[i-1] + x[i+p] - x[i-q]
+    p = (M - 1)/2,   q = p + 1
+    Ref: Udemy Course "Digital Signal Processing (DSP) From Ground Up™ in Python", 
+    available in: https://www.udemy.com/python-dsp/
+    ---------------------------------------------------------------------------
+    X: input signal
+    Fs: Sampling Frequency
+    M: number of points in M range of time the moving average
+    
+    """
+    p = int(((M*Fs)-1)/2)
+    q = p + 1
+    Y = np.zeros(len(X))
+    for i in range(len(X)):
+        Y[i] = Y[i-1] + X[i-p] - X[i-q]
+        
+    return vec_nor(Y)
+
+def butter_bp_coe(lowcut, highcut, fs, order=1):
+    """
+    Butterworth passband filter coefficients b and a
+    Ref: 
+    [1] https://timsainb.github.io/spectrograms-mfccs-and-inversion-in-python.html
+    [2] https://gist.github.com/kastnerkyle/179d6e9a88202ab0a2fe
+    """
+    nyq = 0.5 * fs
+    low = lowcut / nyq
+    high = highcut / nyq
+    b, a = butter(order, [low, high], btype='band')
+    return b, a
+
+def butter_bp_fil(data, lowcut, highcut, fs, order=1):
+    """
+    Butterworth passband filter
+    Ref: 
+    [1] https://timsainb.github.io/spectrograms-mfccs-and-inversion-in-python.html
+    [2] https://gist.github.com/kastnerkyle/179d6e9a88202ab0a2fe
+    """
+    b, a = butter_bp_coe(lowcut, highcut, fs, order=order)
+    y = lfilter(b, a, data)
+    return vec_nor(y)
+
+
+def Fpass(X,lp):
+    """
+    Fpass is the function to pass the coefficients of a filter trough a signal'
+    """
+    llp=np.size(lp)	  	        # Get the length of the lowpass vector		
+
+    x=np.convolve(X,lp);		        # Disrete convolution 
+    x=x[int(llp/2):-int(llp/2)];
+    x=x-(np.mean(x));
+    x=x/np.max(x);
+    
+    y=vec_nor(x);				# Vector Normalizing
+        
+    return y
+
+def FpassBand(X,hp,lp):
+    """
+    FpassBand is the function that develop a pass band filter of the signal 'x' through the
+    discrete convolution of this 'x' first with the coeficients of a High Pass Filter 'hp' and then
+    with the discrete convolution of this result with a Low Pass Filter 'lp'
+    """
+    llp=np.shape(lp)	  	        # Get the length of the lowpass vector		
+    llp=llp[0];				       # Get the value of the length
+    lhp=np.shape(hp)			    # Get the length of the highpass vector		
+    lhp=lhp[0];				       # Get the value of the length	
+
+    x=np.convolve(X,lp);		        # Disrete convolution 
+    x=x[int(llp/2):-int(llp/2)];
+    x=x-(np.mean(x));
+    x=x/np.max(x);
+	
+    y=np.convolve(x,hp);			# Disrete onvolution
+    y=y[int(lhp/2):-int(lhp/2)];
+    y=y-np.mean(y);
+    y=y/np.max(y);
+
+    x=np.convolve(y,lp);		        # Disrete convolution 
+    x=x[int(llp/2):-int(llp/2)];
+    x=x-(np.mean(x));
+    x=x/np.max(x);
+	
+    y=np.convolve(x,hp);			# Disrete onvolution
+    y=y[int(lhp/2):-int(lhp/2)];
+    y=y-np.mean(y);
+    y=y/np.max(y);
+        
+    y=vec_nor(y);				# Vector Normalizing
+        
+    return y
+
+def bandpass_filter(data, lowcut, highcut, signal_freq, filter_order):
+    """
+    Method responsible for creating and applying Butterworth filter.
+    :param deque data: raw data
+    :param float lowcut: filter lowcut frequency value
+    :param float highcut: filter highcut frequency value
+    :param int signal_freq: signal frequency in samples per second (Hz)
+    :param int filter_order: filter order
+    :return array: filtered data
+    """
+    nyquist_freq = 0.5 * signal_freq
+    low = lowcut / nyquist_freq
+    high = highcut / nyquist_freq
+    b, a = butter(filter_order, [low, high], btype="band")
+    y = lfilter(b, a, data)
+
+    return y
+def FpassBand_1(X,Fs, f1, f2):
+    """
+    Ref: http://scipy-cookbook.readthedocs.io/items/FIRFilter.html
+    http://lagrange.univ-lyon1.fr/docs/scipy/0.17.1/generated/scipy.signal.firwin.html
+    FpassBand_1 is a function to develop a passband filter using 'firwin'
+    and 'lfilter' functions from the "Scipy" library
+    """
+
+    # The Nyquist rate of the signal.
+    nyq_rate = Fs / 2.0
+    # The desired width of the transition from pass to stop,
+    # relative to the Nyquist rate.  We'll design the filter
+    # with a 5 Hz transition width.
+    width = 5.0/nyq_rate
+    # The desired attenuation in the stop band, in dB.
+    ripple_db = 60.0
+    # Compute the order and Kaiser parameter for the FIR filter.
+    N, beta = kaiserord(ripple_db, width)
+    # 
+    taps = firwin(100, [f1, f2], pass_zero=False)
+#    taps = firwin(N, L_cutoff_hz/nyq_rate, window=('kaiser', beta))
+#    taps_2 = firwin(N, H_cutoff_hz/nyq_rate, pass_zero=True)
+    # Use lfilter to filter x with the FIR filter.
+    X_pb= lfilter(taps, 1.0, X)
+   # X_pb= lfilter(taps_2, 1.0, X_l)
+    
+    return X_pb[N-1:]
+
+def FpassBand_2(X, f1, f2, Fs):
+    
+    Y = FhighPass(X, Fs, f1)
+    
+    Y = FlowPass(X, Fs, f2)
+    
+    return Y
+
+def FhighPass(X, Fs, H_cutoff_hz):
+    """
+    Ref: http://scipy-cookbook.readthedocs.io/items/FIRFilter.html
+    http://lagrange.univ-lyon1.fr/docs/scipy/0.17.1/generated/scipy.signal.firwin.html
+    FhighPass is a function to develop a highpass filter using 'firwin'
+    and 'lfilter' functions from the "Scipy" library
+    """
+    # The Nyquist rate of the signal.
+    nyq_rate = Fs / 2.0
+    # The desired width of the transition from pass to stop,
+    # relative to the Nyquist rate.  We'll design the filter
+    # with a 5 Hz transition width.
+    width = 5.0/nyq_rate
+    # The desired attenuation in the stop band, in dB.
+    ripple_db = 60.0
+    # Compute the order and Kaiser parameter for the FIR filter.
+    N, beta = kaiserord(ripple_db, width)
+    # Use firwin with a Kaiser window to create a lowpass FIR filter.
+    taps_2 = firwin(N, H_cutoff_hz/nyq_rate, pass_zero=False)
+    # Use lfilter to filter x with the FIR filter.
+    X_h= lfilter(taps_2, 1.0, X)
+    
+    return X_h[N-1:]
+    
+def FlowPass(X, Fs, L_cutoff_hz):
+    """
+    Ref: http://scipy-cookbook.readthedocs.io/items/FIRFilter.html
+    http://lagrange.univ-lyon1.fr/docs/scipy/0.17.1/generated/scipy.signal.firwin.html
+    FlowPass is a function to develop a lowpass filter using 'firwin'
+    and 'lfilter' functions from the "Scipy" library
+    """
+    # The Nyquist rate of the signal.
+    nyq_rate = Fs / 2.0
+    # The desired width of the transition from pass to stop,
+    # relative to the Nyquist rate.  We'll design the filter
+    # with a 5 Hz transition width.
+    width = 5.0/nyq_rate
+    # The desired attenuation in the stop band, in dB.
+    ripple_db = 60.0
+    # Compute the order and Kaiser parameter for the FIR filter.
+    N, beta = kaiserord(ripple_db, width)
+    # Use firwin with a Kaiser window to create a lowpass FIR filter.
+    taps = firwin(N, L_cutoff_hz/nyq_rate, window=('kaiser', beta))
+    # Use lfilter to filter x with the FIR filter.
+    X_l= lfilter(taps, 1.0, X)
+    
+    return X_l[N-1:]
+
+# -----------------------------------------------------------------------------
+           # Segmentation By Running Sum Algorithm & Filters
+# -----------------------------------------------------------------------------
+def seg_RSA1(x, fs):
+    # # Appliying Running Sum Algorithm of PCG filtered from 0.01Hz to 1kHz
+    F_x = running_sum(vec_nor(butter_bp_fil(x, 0.01, 50, fs)))
+    # Smoothing the signal by filtering from 0.01Hz to 5Hz
+    F_x = butter_bp_fil(F_x,0.01, 2, fs)
+    # Appliying 1st derivative to indentify slope sign changes
+    xx = derivate_1(F_x)
+    #
+    xx = butter_bp_fil(xx, 0.01,2, fs)
+    # Transforming positives to 1 & negatives to -1
+    xxS = np.sign(xx)
+    
+    return xxS, F_x
+# -----------------------------------------------------------------------------
+                   # Segmentation By Derivatives
+# -----------------------------------------------------------------------------                                
+def seg_Der1(x, fs):
+    # Segmenting in an specific frequency band
+    pcgPeaks, peaks, allpeaks = PDP(x, fs)
+    # -----------------------------------------------------------------------------
+    # Time processing
+    dT = 0.4          # [1] mean S1 duration 122ms, mean S2 duration 92ms
+    timeV = []
+    timeV.append(0)
+    pointV = []
+    pointV.append(0)
+    segmV = np.zeros(len(x))
+    
+    for i in range(len(pcgPeaks)-1):
+        if pcgPeaks[i]>0.5:
+            timeV.append(i/fs)       # Gives the time when a peak get found
+            pointV.append(i)          # Gives the time when a peak get found
+            if (pointV[-1]/fs)-(pointV[-2]/fs)> dT:
+                segmV[pointV[-2]:pointV[-1]] = 0.4  # Marc a diastolic segment
+            else:
+                segmV[pointV[-2]:pointV[-1]] = 0.6  # Marc a systolic segment
+    
+    return segmV, pcgPeaks
+
+def seg_Der2(x, fs):
+    # Appliying Running Sum Algorithm of PCG filtered from 0.01Hz to 1kHz
+    F_x = running_sum(vec_nor(butter_bp_fil(x, 0.01, 999, fs)))
+    # Smoothing the signal by filtering from 0.01Hz to 5Hz
+    F_x = butter_bp_fil(F_x,0.01, 5, fs)
+    # Time to be represented in samples
+    time_samples = 0.5
+    # Number of samples to move over the signal
+    mC = int(time_samples * fs)                       
+    peaks = np.zeros(len(F_x))
+    p = sg.find_peaks(F_x, distance=mC)
+    # Defining peaks as +1
+    for i in range (len(p[0][:])):
+        peaks[p[0][i]] = 1
+    
+    return peaks                     
+# -----------------------------------------------------------------------------
+                                # Peak Detection
+# -----------------------------------------------------------------------------
+def PDP(Xf, samplerate):
+    """
+    Peak Detection Process
+    """
+    timeCut = samplerate*0.25                      # Time to count another pulse
+    vCorte = 0.6                                   # Amplitude threshold
+    
+    Xf = vec_nor(Xf)                               # Normalize signal
+    dX = derivate_positive(Xf);				      # Derivate of the signal
+    dX = vec_nor(dX);			                  # Vector Normalizing
+    
+    size=np.shape(Xf)				                 # Rank or dimension of the array
+    fil=size[0];					                     # Number of rows
+ 
+    positive=np.zeros((1,fil+1));                   # Initializating Vector 
+    positive=positive[0];                           # Getting the Vector
+
+    points=np.zeros((1,fil));                       # Initializating the Peak Points Vector
+    points=points[0];                               # Getting the point vector
+
+    points1=np.zeros((1,fil));                      # Initializating the Peak Points Vector
+    points1=points1[0];                             # Getting the point vector
+       
+    '''
+    FIRST! having the positives values of the slope as 1
+    And the negative values of the slope as 0
+    '''
+    for i in range(0,fil):
+        if Xf[i] > 0:
+            if dX[i]>0:
+                positive[i] = Xf[i];
+            else:
+                positive[i] = 0;
+    '''
+    SECOND! a peak will be found when the ith value is equal to 1 &&
+    the ith+1 is equal to 0
+    '''
+    for i in range(0,fil):
+        if (positive[i]==Xf[i] and positive[i+1]==0):
+            points[i] = Xf[i];
+        else:
+            points[i] = 0;
+    '''
+    THIRD! Define a minimun Peak Height
+    '''
+    p=0;
+    for i in range(0,fil):
+        if (Xf[i] > vCorte and p==0):
+            p = i
+            points1[i] = Xf[i]
+        else:
+            points1[i] = 0
+            if (p+timeCut < i):
+                p = 0
+                    
+    return points1, points, positive[0:(len(positive)-1)]
+# -----------------------------------------------------------------------------
+# Peak Detection 2
+# -----------------------------------------------------------------------------
+def PDP_2(Xf, samplerate):
+    """
+    Peak Detection Process
+    """
+    timeCut = samplerate*0.25                      # Time to count another pulse
+    vCorte = 0.6                                   # Amplitude threshold
+    
+    #Xf = vec_nor(Xf)                               # Normalize signal
+    dX = np.diff(Xf);				      # Derivate of the signal
+    dX = vec_nor(dX);			                  # Vector Normalizing
+    
+    size=np.shape(Xf)				                 # Rank or dimension of the array
+    fil=size[0];					                     # Number of rows
+ 
+    positive=np.zeros((1,fil+1));                   # Initializating Vector 
+    positive=positive[0];                           # Getting the Vector
+
+    points=np.zeros((1,fil));                       # Initializating the Peak Points Vector
+    points=points[0];                               # Getting the point vector
+
+    points1=np.zeros((1,fil));                      # Initializating the Peak Points Vector
+    points1=points1[0];                             # Getting the point vector
+       
+    '''
+    FIRST! having the positives values of the slope as 1
+    And the negative values of the slope as 0
+    '''
+    for i in range(0,fil):
+        if Xf[i] > 0:
+            if dX[i]>0:
+                positive[i] = Xf[i];
+            else:
+                positive[i] = 0;
+    '''
+    SECOND! a peak will be found when the ith value is equal to 1 &&
+    the ith+1 is equal to 0
+    '''
+    for i in range(0,fil):
+        if (positive[i]==Xf[i] and positive[i+1]==0):
+            points[i] = Xf[i];
+        else:
+            points[i] = 0;
+    '''
+    THIRD! Define a minimun Peak Height
+    '''
+    p=0;
+    for i in range(0,fil):
+        if (Xf[i] > vCorte and p==0):
+            p = i
+            points1[i] = Xf[i]
+        else:
+            points1[i] = 0
+            if (p+timeCut < i):
+                p = 0
+                    
+    return points1, points, positive[0:(len(positive)-1)]
+# -----------------------------------------------------------------------------
+# Peak Detection 3 - Janko Slavic
+# -----------------------------------------------------------------------------
+def findpeaks(data, spacing=1, limit=None):
+    """
+    Janko Slavic peak detection algorithm and implementation.
+    https://github.com/jankoslavic/py-tools/tree/master/findpeaks
+    Finds peaks in `data` which are of `spacing` width and >=`limit`.
+    :param ndarray data: data
+    :param float spacing: minimum spacing to the next peak (should be 1 or more)
+    :param float limit: peaks should have value greater or equal
+    :return array: detected peaks indexes array
+    """
+    len = data.size
+    x = np.zeros(len + 2 * spacing)
+    x[:spacing] = data[0] - 1.e-6
+    x[-spacing:] = data[-1] - 1.e-6
+    x[spacing:spacing + len] = data
+    peak_candidate = np.zeros(len)
+    peak_candidate[:] = True
+    for s in range(spacing):
+        start = spacing - s - 1
+        h_b = x[start: start + len]  # before
+        start = spacing
+        h_c = x[start: start + len]  # central
+        start = spacing + s + 1
+        h_a = x[start: start + len]  # after
+        peak_candidate = np.logical_and(peak_candidate, np.logical_and(h_c > h_b, h_c > h_a))
+
+    ind = np.argwhere(peak_candidate)
+    ind = ind.reshape(ind.size)
+    if limit is not None:
+        ind = ind[data[ind] > limit]
+    return ind
+
+# -----------------------------------------------------------------------------
+                              # Fast Fourier Transform
+# -----------------------------------------------------------------------------
+def fft_k(data, samplerate, showFrequency):
+    '''
+    Fast Fourier Transform using fftpack from scipy
+    Ref: https://docs.scipy.org/doc/scipy/reference/tutorial/fftpack.html
+    '''
+    # FFT Full Vector 'k' coefficients
+    pcgFFT = fft(data)
+    # FFT positives values                                                    
+    short_pcgFFT = 2.0/np.size(data) * np.abs(pcgFFT[0:(np.size(data)//2)])
+    #short_pcgFFT = 2.0/np.size(data) * np.abs(pcgFFT[(np.size(data)//2):None])  vector selected from the middle to the end
+    # Vector of frequencies (X-axes)
+    vTfft = np.linspace(0.0, 1.0/(2.0*(1/samplerate)), np.size(data)//2)  
+    # find the value closest to a value   
+    idx = (np.abs(vTfft-showFrequency)).argmin()             
+    
+    return short_pcgFFT[0:idx], vTfft[0:idx]
+
+def fft_k_N(data, samplerate, showFrequency):
+    '''
+    Normalized Fast Fourier Transform
+    Ref: https://docs.scipy.org/doc/scipy/reference/tutorial/fftpack.html
+    '''
+    # FFT Full Vector 'k' coefficients
+    pcgFFT = fft(data)
+    # FFT positives values from the middle to the end (to evoid the interference at beginning)
+    short_pcgFFT = 2.0/np.size(data) * np.abs(pcgFFT[0:(np.size(data)//2)])
+    #short_pcgFFT = 2.0/np.size(data) * np.abs(pcgFFT[(np.size(data)//2):None])  vector selected from the middle to the end
+    # Vector of frequencies (X-axes)
+    vTfft = np.linspace(0.0, 1.0/(2.0*(1/samplerate)), np.size(data)//2)  
+    # find the value closest to a value   
+    idx = (np.abs(vTfft-showFrequency)).argmin()            
+    
+    return vec_nor(short_pcgFFT[0:idx]), vTfft[0:idx]
+# -----------------------------------------------------------------------------
+                              # PCG Audio Pre-Processing
+# -----------------------------------------------------------------------------
+def pre_pro_audio_PCG(x, fs):
+# Ensure having a Mono sound
+    if len(x.shape)>1:
+    # select the left size
+        x = x[:,0]
+    # Resampling Audio PCG to 2k Hz
+    Frs = 2000
+    Nrs = int(Frs*(len(x)/fs))
+    if fs > Frs:
+        x = sg.resample(x, Nrs)
+    
+    return vec_nor(x), Frs
+# -----------------------------------------------------------------------------
+#             Pre-Process: Signal basic information (SBI)
+# -----------------------------------------------------------------------------
+def pre_pro_basicInfo_PCG(x, fs):
+    # find the time duration of the sound
+    t_sound = 1/fs*len(x)
+    # make a vector time for the sound
+    vec_time = np.linspace(0, t_sound, len(x))
+    return t_sound, vec_time
\ No newline at end of file