docs: 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