docs: def power_spectral_density(x, time_step, freq_range = None):
      """
      returns the *single sided* power spectral density of the time trace x which is sampled at intervals time_step
      Args:
          x (array):  timetrace
          time_step (float): sampling interval of x
          freq_range (array or tuple): frequency range in the form [f_min, f_max] to return only the spectrum within this range
  
      Returns:
  
      """
      N = len(x)
      P = 2 * np.abs(np.fft.rfft(x))**2 / N * time_step
      F = np.fft.rfftfreq(len(x), time_step)
  
      if freq_range is not None:
          brange = np.all([F >= freq_range[0], F <= freq_range[1]], axis=0)
          P = P[brange]
          F = F[brange]
  
      return F, P