code: def numpy_psd(x, f_sample=1.0):
      
      psd_of_x = (2.0/ (float(len(x)) * f_sample)) * numpy.abs(numpy.fft.rfft(x))**2
      f_axis = numpy.linspace(0, f_sample/2.0, len(psd_of_x)) # frequency axis
      return f_axis, psd_of_x
docs: calculate power spectral density of input signal x
          x = signal
          f_sample = sampling frequency in Hz. i.e. 1/fs is the time-interval
               in seconds between datapoints
          scale fft so that output corresponds to 1-sided PSD
          output has units of [X^2/Hz] where X is the unit of x