Caclulate power spectral density with Welch's method

    Args
    ----
    x: ndarray
        sample array
    fs: float
        sampling frequency (1/dt)

    Returns
    -------
    f_welch: ndarray
        Discrete frequencies
    S_xx_welch: ndarray
        Estimated PSD at discrete frequencies `f_welch`
    P_welch: ndarray
        Signal power (integrated PSD)
    df_welch: ndarray
        Delta between discreet frequencies `f_welch`

def calc_PSD_welch(x, fs, nperseg):
    '''Caclulate power spectral density with Welch's method

    Args
    ----
    x: ndarray
        sample array
    fs: float
        sampling frequency (1/dt)

    Returns
    -------
    f_welch: ndarray
        Discrete frequencies
    S_xx_welch: ndarray
        Estimated PSD at discrete frequencies `f_welch`
    P_welch: ndarray
        Signal power (integrated PSD)
    df_welch: ndarray
        Delta between discreet frequencies `f_welch`
    '''
    import numpy
    import scipy.signal

    # Code source and description of FFT, DFT, etc.
    # http://stackoverflow.com/a/33251324/943773
    dt = 1/fs
    N = len(x)
    times = numpy.arange(N) / fs

    # Estimate PSD `S_xx_welch` at discrete frequencies `f_welch`
    f_welch, S_xx_welch = scipy.signal.welch(x, fs=fs, nperseg=nperseg)

    # Integrate PSD over spectral bandwidth
    # to obtain signal power `P_welch`
    df_welch = f_welch[1] - f_welch[0]
    P_welch = numpy.sum(S_xx_welch) * df_welch

    return f_welch, S_xx_welch, P_welch, df_welch