def powerspec(data, tbin, Df=None, units=False, pointProcess=False):
freq, DATA = calculate_fft(data, tbin)<EOL>df = freq[<NUM_LIT:1>] - freq[<NUM_LIT:0>]<EOL>T = tbin * len(freq)<EOL>POW = np.abs(DATA) ** <NUM_LIT:2><EOL>if Df is not None:<EOL><INDENT>POW = [movav(x, Df, df) for x in POW]<EOL>cut = int(Df / df)<EOL>freq = freq[cut:]<EOL>POW = np.array([x[cut:] for x in POW])<EOL>POW = np.abs(POW)<EOL><DEDENT>assert(len(freq) == len(POW[<NUM_LIT:0>]))<EOL>if units is True:<EOL><INDENT>POW = mean(POW, units=units)<EOL>assert(len(freq) == len(POW))<EOL><DEDENT>if pointProcess:<EOL><INDENT>POW *= <NUM_LIT:1.> / T * <NUM_LIT>  <EOL><DEDENT>return freq, POW<EOL>
Calculate (smoothed) power spectra of all timeseries in data.
If units=True, power spectra are averaged across units.
Note that averaging is done on power spectra rather than data.

If pointProcess is True, power spectra are normalized by the length T of the
time series.


Parameters
----------
data : numpy.ndarray
    1st axis unit, 2nd axis time.
tbin : float
    Binsize in ms.
Df : float/None, 
    Window width of sliding rectangular filter (smoothing),
    None is no smoothing.
units : bool
    Average power spectrum.
pointProcess : bool
    If set to True, powerspectrum is normalized to signal length T.


Returns
-------
freq : tuple
    numpy.ndarray of frequencies.
POW : tuple
    if units=False: 
        2 dim numpy.ndarray; 1st axis unit, 2nd axis frequency
    if units=True:  
        1 dim numpy.ndarray; frequency series


Examples
--------    
>>> powerspec(np.array([analog_sig1, analog_sig2]), tbin, Df=Df)
Out[1]: (freq,POW)
>>> POW.shape
Out[2]: (2,len(analog_sig1))

>>> powerspec(np.array([analog_sig1, analog_sig2]), tbin, Df=Df, units=True)
Out[1]: (freq,POW)
>>> POW.shape
Out[2]: (len(analog_sig1),)