#NEWLINE# def _spectral_helper( x , y , fs = 1.0 , window = 'hann' , nperseg = 256 , noverlap = None , nfft = None , detrend = 'constant' , return_onesided = True , scaling = 'spectrum' , axis = ( - 1 ) , mode = 'psd' ) : #NEWLINE# #INDENT# '\n    Calculate various forms of windowed FFTs for PSD, CSD, etc.\n\n    This is a helper function that implements the commonality between the\n    psd, csd, and spectrogram functions. It is not designed to be called\n    externally. The windows are not averaged over; the result from each window\n    is returned.\n\n    Parameters\n    ---------\n    x : array_like\n        Array or sequence containing the data to be analyzed.\n    y : array_like\n        Array or sequence containing the data to be analyzed. If this is\n        the same object in memoery as x (i.e. _spectral_helper(x, x, ...)),\n        the extra computations are spared.\n    fs : float, optional\n        Sampling frequency of the time series. Defaults to 1.0.\n    window : str or tuple or array_like, optional\n        Desired window to use. See `get_window` for a list of windows and\n        required parameters. If `window` is array_like it will be used\n        directly as the window and its length will be used for nperseg.\n        Defaults to \'hann\'.\n    nperseg : int, optional\n        Length of each segment.  Defaults to 256.\n    noverlap : int, optional\n        Number of points to overlap between segments. If None,\n        ``noverlap = nperseg // 2``.  Defaults to None.\n    nfft : int, optional\n        Length of the FFT used, if a zero padded FFT is desired.  If None,\n        the FFT length is `nperseg`. Defaults to None.\n    detrend : str or function or False, optional\n        Specifies how to detrend each segment. If `detrend` is a string,\n        it is passed as the ``type`` argument to `detrend`.  If it is a\n        function, it takes a segment and returns a detrended segment.\n        If `detrend` is False, no detrending is done.  Defaults to \'constant\'.\n    return_onesided : bool, optional\n        If True, return a one-sided spectrum for real data. If False return\n        a two-sided spectrum. Note that for complex data, a two-sided\n        spectrum is always returned.\n    scaling : { \'density\', \'spectrum\' }, optional\n        Selects between computing the cross spectral density (\'density\')\n        where `Pxy` has units of V**2/Hz and computing the cross spectrum\n        (\'spectrum\') where `Pxy` has units of V**2, if `x` and `y` are\n        measured in V and fs is measured in Hz.  Defaults to \'density\'\n    axis : int, optional\n        Axis along which the periodogram is computed; the default is over\n        the last axis (i.e. ``axis=-1``).\n    mode : str, optional\n        Defines what kind of return values are expected. Options are [\'psd\',\n        \'complex\', \'magnitude\', \'angle\', \'phase\'].\n\n    Returns\n    -------\n    freqs : ndarray\n        Array of sample frequencies.\n    t : ndarray\n        Array of times corresponding to each data segment\n    result : ndarray\n        Array of output data, contents dependant on *mode* kwarg.\n\n    References\n    ----------\n    .. [1] Stack Overflow, "Rolling window for 1D arrays in Numpy?",\n        http://stackoverflow.com/a/6811241\n    .. [2] Stack Overflow, "Using strides for an efficient moving average\n        filter", http://stackoverflow.com/a/4947453\n\n    Notes\n    -----\n    Adapted from matplotlib.mlab\n\n    .. versionadded:: 0.16.0\n    ' #NEWLINE# if ( mode not in [ 'psd' , 'complex' , 'magnitude' , 'angle' , 'phase' ] ) : #NEWLINE# #INDENT# raise ValueError ( ( "Unknown value for mode %s, must be one of: 'default', 'psd', 'complex', 'magnitude', 'angle', 'phase'" % mode ) ) #NEWLINE# #UNINDENT# same_data = ( y is x ) #NEWLINE# if ( ( not same_data ) and ( mode != 'psd' ) ) : #NEWLINE# #INDENT# raise ValueError ( "x and y must be equal if mode is not 'psd'" ) #NEWLINE# #UNINDENT# axis = int ( axis ) #NEWLINE# x = np . asarray ( x ) #NEWLINE# if ( not same_data ) : #NEWLINE# #INDENT# y = np . asarray ( y ) #NEWLINE# outdtype = np . result_type ( x , y , np . complex64 ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# outdtype = np . result_type ( x , np . complex64 ) #NEWLINE# #UNINDENT# if ( not same_data ) : #NEWLINE# #INDENT# xouter = list ( x . shape ) #NEWLINE# youter = list ( y . shape ) #NEWLINE# xouter . pop ( axis ) #NEWLINE# youter . pop ( axis ) #NEWLINE# try : #NEWLINE# #INDENT# outershape = np . broadcast ( np . empty ( xouter ) , np . empty ( youter ) ) . shape #NEWLINE# #UNINDENT# except ValueError : #NEWLINE# #INDENT# raise ValueError ( 'x and y cannot be broadcast together.' ) #NEWLINE# #UNINDENT# #UNINDENT# if same_data : #NEWLINE# #INDENT# if ( x . size == 0 ) : #NEWLINE# #INDENT# return ( np . empty ( x . shape ) , np . empty ( x . shape ) , np . empty ( x . shape ) ) #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( ( x . size == 0 ) or ( y . size == 0 ) ) : #NEWLINE# #INDENT# outshape = ( outershape + ( min ( [ x . shape [ axis ] , y . shape [ axis ] ] ) , ) ) #NEWLINE# emptyout = np . rollaxis ( np . empty ( outshape ) , ( - 1 ) , axis ) #NEWLINE# return ( emptyout , emptyout , emptyout ) #NEWLINE# #UNINDENT# #UNINDENT# if ( x . ndim > 1 ) : #NEWLINE# #INDENT# if ( axis != ( - 1 ) ) : #NEWLINE# #INDENT# x = np . rollaxis ( x , axis , len ( x . shape ) ) #NEWLINE# if ( ( not same_data ) and ( y . ndim > 1 ) ) : #NEWLINE# #INDENT# y = np . rollaxis ( y , axis , len ( y . shape ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# if ( not same_data ) : #NEWLINE# #INDENT# if ( x . shape [ ( - 1 ) ] != y . shape [ ( - 1 ) ] ) : #NEWLINE# #INDENT# if ( x . shape [ ( - 1 ) ] < y . shape [ ( - 1 ) ] ) : #NEWLINE# #INDENT# pad_shape = list ( x . shape ) #NEWLINE# pad_shape [ ( - 1 ) ] = ( y . shape [ ( - 1 ) ] - x . shape [ ( - 1 ) ] ) #NEWLINE# x = np . concatenate ( ( x , np . zeros ( pad_shape ) ) , ( - 1 ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# pad_shape = list ( y . shape ) #NEWLINE# pad_shape [ ( - 1 ) ] = ( x . shape [ ( - 1 ) ] - y . shape [ ( - 1 ) ] ) #NEWLINE# y = np . concatenate ( ( y , np . zeros ( pad_shape ) ) , ( - 1 ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# if ( x . shape [ ( - 1 ) ] < nperseg ) : #NEWLINE# #INDENT# warnings . warn ( 'nperseg = {0:d}, is greater than input length = {1:d}, using nperseg = {1:d}' . format ( nperseg , x . shape [ ( - 1 ) ] ) ) #NEWLINE# nperseg = x . shape [ ( - 1 ) ] #NEWLINE# #UNINDENT# nperseg = int ( nperseg ) #NEWLINE# if ( nperseg < 1 ) : #NEWLINE# #INDENT# raise ValueError ( 'nperseg must be a positive integer' ) #NEWLINE# #UNINDENT# if ( nfft is None ) : #NEWLINE# #INDENT# nfft = nperseg #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( nfft < nperseg ) : #NEWLINE# #INDENT# raise ValueError ( 'nfft must be greater than or equal to nperseg.' ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# nfft = int ( nfft ) #NEWLINE# #UNINDENT# #UNINDENT# if ( noverlap is None ) : #NEWLINE# #INDENT# noverlap = ( nperseg // 2 ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( noverlap >= nperseg ) : #NEWLINE# #INDENT# raise ValueError ( 'noverlap must be less than nperseg.' ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# noverlap = int ( noverlap ) #NEWLINE# #UNINDENT# #UNINDENT# if ( not detrend ) : #NEWLINE# #INDENT# def detrend_func( d ) : #NEWLINE# #INDENT# return d #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( not hasattr ( detrend , '__call__' ) ) : #NEWLINE# #INDENT# def detrend_func( d ) : #NEWLINE# #INDENT# return signaltools . detrend ( d , type = detrend , axis = ( - 1 ) ) #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# if ( axis != ( - 1 ) ) : #NEWLINE# #INDENT# def detrend_func( d ) : #NEWLINE# #INDENT# d = np . rollaxis ( d , ( - 1 ) , axis ) #NEWLINE# d = detrend ( d ) #NEWLINE# return np . rollaxis ( d , axis , len ( d . shape ) ) #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# detrend_func = detrend #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# if ( isinstance ( window , string_types ) or ( type ( window ) is tuple ) ) : #NEWLINE# #INDENT# win = get_window ( window , nperseg ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# win = np . asarray ( window ) #NEWLINE# if ( len ( win . shape ) != 1 ) : #NEWLINE# #INDENT# raise ValueError ( 'window must be 1-D' ) #NEWLINE# #UNINDENT# if ( win . shape [ 0 ] != nperseg ) : #NEWLINE# #INDENT# raise ValueError ( 'window must have length of nperseg' ) #NEWLINE# #UNINDENT# #UNINDENT# if ( np . result_type ( win , np . complex64 ) != outdtype ) : #NEWLINE# #INDENT# win = win . astype ( outdtype ) #NEWLINE# #UNINDENT# if ( mode == 'psd' ) : #NEWLINE# #INDENT# if ( scaling == 'density' ) : #NEWLINE# #INDENT# scale = ( 1.0 / ( fs * ( win * win ) . sum ( ) ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( scaling == 'spectrum' ) : #NEWLINE# #INDENT# scale = ( 1.0 / ( win . sum ( ) ** 2 ) ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# raise ValueError ( ( 'Unknown scaling: %r' % scaling ) ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# scale = 1 #NEWLINE# #UNINDENT# if ( return_onesided is True ) : #NEWLINE# #INDENT# if np . iscomplexobj ( x ) : #NEWLINE# #INDENT# sides = 'twosided' #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# sides = 'onesided' #NEWLINE# if ( not same_data ) : #NEWLINE# #INDENT# if np . iscomplexobj ( y ) : #NEWLINE# #INDENT# sides = 'twosided' #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# sides = 'twosided' #NEWLINE# #UNINDENT# if ( sides == 'twosided' ) : #NEWLINE# #INDENT# num_freqs = nfft #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( sides == 'onesided' ) : #NEWLINE# #INDENT# if ( nfft % 2 ) : #NEWLINE# #INDENT# num_freqs = ( ( nfft + 1 ) // 2 ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# num_freqs = ( ( nfft // 2 ) + 1 ) #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# result = _fft_helper ( x , win , detrend_func , nperseg , noverlap , nfft ) #NEWLINE# result = result [ ... , : num_freqs ] #NEWLINE# freqs = fftpack . fftfreq ( nfft , ( 1 / fs ) ) [ : num_freqs ] #NEWLINE# if ( not same_data ) : #NEWLINE# #INDENT# result_y = _fft_helper ( y , win , detrend_func , nperseg , noverlap , nfft ) #NEWLINE# result_y = result_y [ ... , : num_freqs ] #NEWLINE# result = ( np . conjugate ( result ) * result_y ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( mode == 'psd' ) : #NEWLINE# #INDENT# result = ( np . conjugate ( result ) * result ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( mode == 'magnitude' ) : #NEWLINE# #INDENT# result = np . absolute ( result ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( ( mode == 'angle' ) or ( mode == 'phase' ) ) : #NEWLINE# #INDENT# result = np . angle ( result ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# if ( mode == 'complex' ) : #NEWLINE# #INDENT# pass #NEWLINE# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# #UNINDENT# result *= scale #NEWLINE# if ( sides == 'onesided' ) : #NEWLINE# #INDENT# if ( nfft % 2 ) : #NEWLINE# #INDENT# result [ ... , 1 : ] *= 2 #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# result [ ... , 1 : ( - 1 ) ] *= 2 #NEWLINE# #UNINDENT# #UNINDENT# t = ( np . arange ( ( nperseg / 2 ) , ( ( x . shape [ ( - 1 ) ] - ( nperseg / 2 ) ) + 1 ) , ( nperseg - noverlap ) ) / float ( fs ) ) #NEWLINE# if ( ( sides != 'twosided' ) and ( not ( nfft % 2 ) ) ) : #NEWLINE# #INDENT# freqs [ ( - 1 ) ] *= ( - 1 ) #NEWLINE# #UNINDENT# if ( mode == 'phase' ) : #NEWLINE# #INDENT# result = np . unwrap ( result , axis = ( - 1 ) ) #NEWLINE# #UNINDENT# result = result . astype ( outdtype ) #NEWLINE# if ( same_data and ( mode != 'complex' ) ) : #NEWLINE# #INDENT# result = result . real #NEWLINE# #UNINDENT# if ( axis != ( - 1 ) ) : #NEWLINE# #INDENT# if ( axis < 0 ) : #NEWLINE# #INDENT# axis = ( ( len ( result . shape ) - 1 ) - axis ) #NEWLINE# #UNINDENT# result = np . rollaxis ( result , ( - 1 ) , axis ) #NEWLINE# #UNINDENT# else : #NEWLINE# #INDENT# result = np . rollaxis ( result , ( - 1 ) , ( - 2 ) ) #NEWLINE# #UNINDENT# return ( freqs , t , result )