def psd(t, y, pow2=False, window=None, rescale=False):
<EOL>f, Y = fft(t,y,pow2,window,rescale)<EOL>f = _n.abs(f[int(len(f)/<NUM_LIT:2>)::-<NUM_LIT:1>])<EOL>P = _n.abs(Y[int(len(Y)/<NUM_LIT:2>)::-<NUM_LIT:1>])**<NUM_LIT:2> / (f[<NUM_LIT:1>]-f[<NUM_LIT:0>])<EOL>if len(t)%<NUM_LIT:2> == <NUM_LIT:0>: P[<NUM_LIT:1>:len(P)-<NUM_LIT:1>] = P[<NUM_LIT:1>:len(P)-<NUM_LIT:1>]*<NUM_LIT:2><EOL>else:             P[<NUM_LIT:1>:]         = P[<NUM_LIT:1>:]*<NUM_LIT:2><EOL>return f, P<EOL>
Single-sided power spectral density, assuming real valued inputs. This goes 
through the numpy fourier transform process, assembling and returning
(frequencies, psd) given time and signal data y. 

Note it is defined such that sum(psd)*df, where df is the frequency
spacing, is the variance of the original signal for any range of frequencies.
This includes the DC and Nyquist components: 

    sqrt(psd[0]*df)  = average value of original time trace

    sqrt(psd[-1]*df) = amplitude of Nyquist component (for even # points)

Parameters
----------
t,y
    Time (t) and signal (y) arrays with which to perform the PSD. Note the t
    array is assumed to be evenly spaced.

pow2 = False
    Set this to true if you only want to keep the first 2^n data
    points (speeds up the FFT substantially)

window = None
    can be set to any of the windowing functions in numpy,
    e.g. window='hanning'. 

rescale = False
    If True, the FFT will be rescaled by the square root of the ratio of 
    variances before and after windowing, such that the integral 
    sum(PSD)*df is the variance of the *original* time-domain data.

returns frequencies, psd (y^2/Hz)