signature:
def get_converted_psd ( self , sides ) :
    """This function returns the PSD in the * * sides * * format
 : param str sides : the PSD format in [ ' onesided ' , ' twosided ' , ' centerdc ' ]
 : return : the expected PSD .
 . . doctest : :
 from spectrum import *
 p = pcovar ( marple _ data , 15)
 centerdc _ psd = p . get _ converted _ psd ( ' centerdc ' )
 . . note : : this function does not change the object , in particular , it
 does not change the : attr : ` psd ` attribute . If you want to change
 the psd on the fly , change the attribute : attr : ` sides ` ."""

implementation:
if sides == self . sides : # nothing to be done is sides = : attr : ` sides
    return self . __psd
if self . datatype == 'complex' :
    assert sides != 'onesided' , "complex datatype so sides cannot be onesided."
if self . sides == 'onesided' :
    logging . debug ( 'Current sides is onesided' )
    if sides == 'twosided' :
        logging . debug ( '--->Converting to twosided' )
        # here we divide everything by 2 to get the twosided version
        # N = self . NFFT
        newpsd = numpy . concatenate ( ( self . psd [ 0 : - 1 ] / 2. , list ( reversed ( self . psd [ 0 : - 1 ] / 2. ) ) ) )
        # so we need to multiply by 2 the 0 and FS / 2 frequencies
        newpsd [ - 1 ] = self . psd [ - 1 ]
        newpsd [ 0 ] *= 2.
    elif sides == 'centerdc' : # FIXME . this assumes data is even so PSD is stored as
    # P0 X1 X2 X3 P1
        logging . debug ( '--->Converting to centerdc' )
        P0 = self . psd [ 0 ]
        P1 = self . psd [ - 1 ]
        newpsd = numpy . concatenate ( ( self . psd [ - 1 : 0 : - 1 ] / 2. , self . psd [ 0 : - 1 ] / 2. ) )
        # so we need to multiply by 2 the 0 and F2/2 frequencies
        # newpsd [ - 1 ] = P0 / 2
        newpsd [ 0 ] = P1
elif self . sides == 'twosided' :
    logging . debug ( 'Current sides is twosided' )
    if sides == 'onesided' : # we assume that data is stored as X0 , X1 , X2 , X3 , XN
    # that is original data is even .
        logging . debug ( 'Converting to onesided assuming ori data is even' )
        midN = ( len ( self . psd ) - 2 ) / 2
        newpsd = numpy . array ( self . psd [ 0 : int ( midN ) + 2 ] * 2 )
        newpsd [ 0 ] /= 2
        newpsd [ - 1 ] = self . psd [ - 1 ]
    elif sides == 'centerdc' :
        newpsd = stools . twosided_2_centerdc ( self . psd )
elif self . sides == 'centerdc' : # same as twosided to onesided
    logging . debug ( 'Current sides is centerdc' )
    if sides == 'onesided' :
        logging . debug ( '--->Converting to onesided' )
        midN = int ( len ( self . psd ) / 2 )
        P1 = self . psd [ 0 ]
        newpsd = numpy . append ( self . psd [ midN : ] * 2 , P1 )
    elif sides == 'twosided' :
        newpsd = stools . centerdc_2_twosided ( self . psd )
else :
    raise ValueError ( "sides must be set to 'onesided', 'twosided' or 'centerdc'" )
return newpsd


signature:
def plot ( self , filename = None , norm = False , ylim = None , sides = None , ** kargs ) :
    """a simple plotting routine to plot the PSD versus frequency .
 : param str filename : save the figure into a file
 : param norm : False by default . If True , the PSD is normalised .
 : param ylim : readjust the y range .
 : param sides : if not provided , : attr : ` sides ` is used . See : attr : ` sides `
 for details .
 : param kargs : any optional argument accepted by : func : ` pylab . plot ` .
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import *
 p = Periodogram ( marple _ data )
 p . plot ( norm = True , marker = ' o ' )"""

implementation:
import pylab
from pylab import ylim as plt_ylim
# First , check that psd attribute is up - to - date
# just to get the PSD to be recomputed if needed
_ = self . psd
# check that the input sides parameter is correct if provided
if sides is not None :
    if sides not in self . _sides_choices :
        raise errors . SpectrumChoiceError ( sides , self . _sides_choices )
# if sides is provided but identical to the current psd , nothing to do .
# if sides not provided , let us use self . sides
if sides is None or sides == self . sides :
    frequencies = self . frequencies ( )
    psd = self . psd
    sides = self . sides
elif sides is not None : # if sides argument is different from the attribute , we need to
# create a new PSD / Freq ; indeed we do not want to change the
# attribute itself
# if data is complex , one - sided is wrong in any case .
    if self . datatype == 'complex' :
        if sides == 'onesided' :
            raise ValueError ( "sides cannot be one-sided with complex data" )
    logging . debug ( "sides is different from the one provided. Converting PSD" )
    frequencies = self . frequencies ( sides = sides )
    psd = self . get_converted_psd ( sides )
if len ( psd ) != len ( frequencies ) :
    raise ValueError ( "PSD length is %s and freq length is %s" % ( len ( psd ) , len ( frequencies ) ) )
if 'ax' in list ( kargs . keys ( ) ) :
    save_ax = pylab . gca ( )
    pylab . sca ( kargs [ 'ax' ] )
    rollback = True
    del kargs [ 'ax' ]
else :
    rollback = False
if norm :
    pylab . plot ( frequencies , 10 * stools . log10 ( psd / max ( psd ) ) , ** kargs )
else :
    pylab . plot ( frequencies , 10 * stools . log10 ( psd ) , ** kargs )
pylab . xlabel ( 'Frequency' )
pylab . ylabel ( 'Power (dB)' )
pylab . grid ( True )
if ylim :
    plt_ylim ( ylim )
if sides == 'onesided' :
    pylab . xlim ( 0 , self . sampling / 2. )
elif sides == 'twosided' :
    pylab . xlim ( 0 , self . sampling )
elif sides == 'centerdc' :
    pylab . xlim ( - self . sampling / 2. , self . sampling / 2. )
if filename :
    pylab . savefig ( filename )
if rollback :
    pylab . sca ( save_ax )
del psd , frequencies


signature:
def power ( self ) :
    r"""Return the power contained in the PSD
 if scale _ by _ freq is False , the power is :
 . . math : : P = N \ sum _ { k = 1 } ^ { N } P _ { xx } ( k )
 else , it is
 . . math : : P = \ sum _ { k = 1 } ^ { N } P _ { xx } ( k ) \ frac { df } { 2 \ pi }
 . . todo : : check these equations"""

implementation:
if self . scale_by_freq == False :
    return sum ( self . psd ) * len ( self . psd )
else :
    return sum ( self . psd ) * self . df / ( 2. * numpy . pi )


signature:
def periodogram ( self ) :
    """An alias to : class : ` ~ spectrum . periodogram . Periodogram `
 The parameters are extracted from the attributes . Relevant attributes
 ares : attr : ` window ` , attr : ` sampling ` , attr : ` NFFT ` , attr : ` scale _ by _ freq ` ,
 : attr : ` detrend ` .
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import datasets
 from spectrum import FourierSpectrum
 s = FourierSpectrum ( datasets . data _ cosine ( ) , sampling = 1024 , NFFT = 512)
 s . periodogram ( )
 s . plot ( )"""

implementation:
from . periodogram import speriodogram
psd = speriodogram ( self . data , window = self . window , sampling = self . sampling , NFFT = self . NFFT , scale_by_freq = self . scale_by_freq , detrend = self . detrend )
self . psd = psd


signature:
def ipy_notebook_skeleton ( ) :
    """Returns a dictionary with the elements of a Jupyter notebook"""

implementation:
py_version = sys . version_info
notebook_skeleton = { "cells" : [ ] , "metadata" : { "kernelspec" : { "display_name" : "Python " + str ( py_version [ 0 ] ) , "language" : "python" , "name" : "python" + str ( py_version [ 0 ] ) } , "language_info" : { "codemirror_mode" : { "name" : "ipython" , "version" : py_version [ 0 ] } , "file_extension" : ".py" , "mimetype" : "text/x-python" , "name" : "python" , "nbconvert_exporter" : "python" , "pygments_lexer" : "ipython" + str ( py_version [ 0 ] ) , "version" : '{0}.{1}.{2}' . format ( * sys . version_info [ : 3 ] ) } } , "nbformat" : 4 , "nbformat_minor" : 0 }
return notebook_skeleton


signature:
def rst2md ( text ) :
    """Converts the RST text from the examples docstrigs and comments
 into markdown text for the IPython notebooks"""

implementation:
top_heading = re . compile ( r'^=+$\s^([\w\s-]+)^=+$' , flags = re . M )
text = re . sub ( top_heading , r'# \1' , text )
math_eq = re . compile ( r'^\.\. math::((?:.+)?(?:\n+^  .+)*)' , flags = re . M )
text = re . sub ( math_eq , lambda match : r'$${0}$$' . format ( match . group ( 1 ) . strip ( ) ) , text )
inline_math = re . compile ( r':math:`(.+)`' )
text = re . sub ( inline_math , r'$\1$' , text )
return text


signature:
def add_markdown_cell ( self , text ) :
    """Add a markdown cell to the notebook
 Parameters
 code : str
 Cell content"""

implementation:
markdown_cell = { "cell_type" : "markdown" , "metadata" : { } , "source" : [ rst2md ( text ) ] }
self . work_notebook [ "cells" ] . append ( markdown_cell )


signature:
def save_file ( self ) :
    """Saves the notebook to a file"""

implementation:
with open ( self . write_file , 'w' ) as out_nb :
    json . dump ( self . work_notebook , out_nb , indent = 2 )


signature:
def arma2psd ( A = None , B = None , rho = 1. , T = 1. , NFFT = 4096 , sides = 'default' , norm = False ) :
    r"""Computes power spectral density given ARMA values .
 This function computes the power spectral density values
 given the ARMA parameters of an ARMA model . It assumes that
 the driving sequence is a white noise process of zero mean and
 variance : math : ` \ rho _ w ` . The sampling frequency and noise variance are
 used to scale the PSD output , which length is set by the user with the
 ` NFFT ` parameter .
 : param array A : Array of AR parameters ( complex or real )
 : param array B : Array of MA parameters ( complex or real )
 : param float rho : White noise variance to scale the returned PSD
 : param float T : Sample interval in seconds to scale the returned PSD
 : param int NFFT : Final size of the PSD
 : param str sides : Default PSD is two - sided , but sides can be set to centerdc .
 . . warning : : By convention , the AR or MA arrays does not contain the
 A0 = 1 value .
 If : attr : ` B ` is None , the model is a pure AR model . If : attr : ` A ` is None ,
 the model is a pure MA model .
 : return : two - sided PSD
 . . rubric : : Details :
 AR case : the power spectral density is :
 . . math : : P _ { ARMA } ( f ) = T \ rho _ w \ left | \ frac { B ( f ) } { A ( f ) } \ right | ^ 2
 where :
 . . math : : A ( f ) = 1 + \ sum _ { k = 1 } ^ q b ( k ) e ^ { - j2 \ pi fkT }
 . . math : : B ( f ) = 1 + \ sum _ { k = 1 } ^ p a ( k ) e ^ { - j2 \ pi fkT }
 . . rubric : : * * Example : * *
 . . plot : :
 : width : 80%
 : include - source :
 import spectrum . arma
 from pylab import plot , log10 , legend
 plot ( 10 * log10 ( spectrum . arma . arma2psd ( [ 1,0.5 ] , [ 0.5,0.5 ] ) ) , label = ' ARMA ( 2,2 ) ' )
 plot ( 10 * log10 ( spectrum . arma . arma2psd ( [ 1,0.5 ] , None ) ) , label = ' AR ( 2 ) ' )
 plot ( 10 * log10 ( spectrum . arma . arma2psd ( None , [ 0.5,0.5 ] ) ) , label = ' MA ( 2 ) ' )
 legend ( )
 : References : [ Marple ] _"""

implementation:
if NFFT is None :
    NFFT = 4096
if A is None and B is None :
    raise ValueError ( "Either AR or MA model must be provided" )
psd = np . zeros ( NFFT , dtype = complex )
if A is not None :
    ip = len ( A )
    den = np . zeros ( NFFT , dtype = complex )
    den [ 0 ] = 1. + 0j
    for k in range ( 0 , ip ) :
        den [ k + 1 ] = A [ k ]
    denf = fft ( den , NFFT )
if B is not None :
    iq = len ( B )
    num = np . zeros ( NFFT , dtype = complex )
    num [ 0 ] = 1. + 0j
    for k in range ( 0 , iq ) :
        num [ k + 1 ] = B [ k ]
    numf = fft ( num , NFFT )
# Changed in version 0.6.9 ( divided by T instead of multiply )
if A is not None and B is not None :
    psd = rho / T * abs ( numf ) ** 2. / abs ( denf ) ** 2.
elif A is not None :
    psd = rho / T / abs ( denf ) ** 2.
elif B is not None :
    psd = rho / T * abs ( numf ) ** 2.
psd = np . real ( psd )
# The PSD is a twosided PSD .
# to obtain the centerdc
if sides != 'default' :
    from . import tools
    assert sides in [ 'centerdc' ]
    if sides == 'centerdc' :
        psd = tools . twosided_2_centerdc ( psd )
if norm == True :
    psd /= max ( psd )
return psd


signature:
def arma_estimate ( X , P , Q , lag ) :
    """Autoregressive and moving average estimators .
 This function provides an estimate of the autoregressive
 parameters , the moving average parameters , and the driving
 white noise variance of an ARMA ( P , Q ) for a complex or real data sequence .
 The parameters are estimated using three steps :
 * Estimate the AR parameters from the original data based on a least
 squares modified Yule - Walker technique ,
 * Produce a residual time sequence by filtering the original data
 with a filter based on the AR parameters ,
 * Estimate the MA parameters from the residual time sequence .
 : param array X : Array of data samples ( length N )
 : param int P : Desired number of AR parameters
 : param int Q : Desired number of MA parameters
 : param int lag : Maximum lag to use for autocorrelation estimates
 : return :
 * A - Array of complex P AR parameter estimates
 * B - Array of complex Q MA parameter estimates
 * RHO - White noise variance estimate
 . . note : :
 * lag must be > = Q ( MA order )
 * * dependencies * * :
 * : meth : ` spectrum . correlation . CORRELATION `
 * : meth : ` spectrum . covar . arcovar `
 * : meth : ` spectrum . arma . ma `
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import arma _ estimate , arma2psd , marple _ data
 import pylab
 a , b , rho = arma _ estimate ( marple _ data , 15 , 15 , 30)
 psd = arma2psd ( A = a , B = b , rho = rho , sides = ' centerdc ' , norm = True )
 pylab . plot ( 10 * pylab . log10 ( psd ) )
 pylab . ylim ( [ - 50,0 ] )
 : reference : [ Marple ] _"""

implementation:
R = CORRELATION ( X , maxlags = lag , norm = 'unbiased' )
R0 = R [ 0 ]
# C Estimate the AR parameters ( no error weighting is used ) .
# C Number of equation errors is M - Q .
MPQ = lag - Q + P
N = len ( X )
Y = np . zeros ( N - P , dtype = complex )
for K in range ( 0 , MPQ ) :
    KPQ = K + Q - P + 1
    if KPQ < 0 :
        Y [ K ] = R [ - KPQ ] . conjugate ( )
    if KPQ == 0 :
        Y [ K ] = R0
    if KPQ > 0 :
        Y [ K ] = R [ KPQ ]
# The resize is very important for the normalissation .
Y . resize ( lag )
if P <= 4 :
    res = arcovar_marple ( Y . copy ( ) , P )
    # ! Eq . ( 10.12)
    ar_params = res [ 0 ]
else :
    res = arcovar ( Y . copy ( ) , P )
    # ! Eq . ( 10.12)
    ar_params = res [ 0 ]
# the . copy is used to prevent a reference somewhere . this is a bug
# to be tracked down .
Y . resize ( N - P )
# C Filter the original time series
for k in range ( P , N ) :
    SUM = X [ k ]
    # SUM + = sum ( [ ar _ params [ j ] * X [ k - j - 1 ] for j in range ( 0 , P ) ] )
    for j in range ( 0 , P ) :
        SUM = SUM + ar_params [ j ] * X [ k - j - 1 ]
        # ! Eq . ( 10.17)
    Y [ k - P ] = SUM
# Estimate the MA parameters ( a " long " AR of order at least 2 * IQ
# C is suggested )
# Y . resize ( N - P )
ma_params , rho = ma ( Y , Q , 2 * Q )
# ! Eq . ( 10.3)
return ar_params , ma_params , rho


signature:
def ma ( X , Q , M ) :
    """Moving average estimator .
 This program provides an estimate of the moving average parameters
 and driving noise variance for a data sequence based on a
 long AR model and a least squares fit .
 : param array X : The input data array
 : param int Q : Desired MA model order ( must be > 0 and < M )
 : param int M : Order of " long " AR model ( suggest at least 2 * Q )
 : return :
 * MA - Array of Q complex MA parameter estimates
 * RHO - Real scalar of white noise variance estimate
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import arma2psd , ma , marple _ data
 import pylab
 # Estimate 15 Ma parameters
 b , rho = ma ( marple _ data , 15 , 30)
 # Create the PSD from those MA parameters
 psd = arma2psd ( B = b , rho = rho , sides = ' centerdc ' )
 # and finally plot the PSD
 pylab . plot ( pylab . linspace ( - 0.5 , 0.5 , 4096 ) , 10 * pylab . log10 ( psd / max ( psd ) ) )
 pylab . axis ( [ - 0.5 , 0.5 , - 30 , 0 ] )
 : reference : [ Marple ] _"""

implementation:
if Q <= 0 or Q >= M :
    raise ValueError ( 'Q(MA) must be in ]0,lag[' )
# C Fit a high - order AR to the data
a , rho , _c = yulewalker . aryule ( X , M , 'biased' )
# ! Eq . ( 10.5)
# add an element unity to the AR parameter array
a = np . insert ( a , 0 , 1 )
# C Find MA parameters from autocorrelations by Yule - Walker method
ma_params , _p , _c = yulewalker . aryule ( a , Q , 'biased' )
# ! Eq . ( 10.7)
return ma_params , rho


signature:
def CORRELOGRAMPSD ( X , Y = None , lag = - 1 , window = 'hamming' , norm = 'unbiased' , NFFT = 4096 , window_params = { } , correlation_method = 'xcorr' ) :
    """PSD estimate using correlogram method .
 : param array X : complex or real data samples X ( 1 ) to X ( N )
 : param array Y : complex data samples Y ( 1 ) to Y ( N ) . If provided , computes
 the cross PSD , otherwise the PSD is returned
 : param int lag : highest lag index to compute . Must be less than N
 : param str window _ name : see : mod : ` window ` for list of valid names
 : param str norm : one of the valid normalisation of : func : ` xcorr ` ( biased ,
 unbiased , coeff , None )
 : param int NFFT : total length of the final data sets ( padded with zero
 if needed ; default is 4096)
 : param str correlation _ method : either ` xcorr ` or ` CORRELATION ` .
 CORRELATION should be removed in the future .
 : return :
 * Array of real ( cross ) power spectral density estimate values . This is
 a two sided array with negative values following the positive ones
 whatever is the input data ( real or complex ) .
 . . rubric : : Description :
 The exact power spectral density is the Fourier transform of the
 autocorrelation sequence :
 . . math : : P _ { xx } ( f ) = T \ sum _ { m = - \ infty } ^ { \ infty } r _ { xx } [ m ] exp ^ { - j2 \ pi fmT }
 The correlogram method of PSD estimation substitutes a finite sequence of
 autocorrelation estimates : math : ` \ hat { r } _ { xx } ` in place of : math : ` r _ { xx } ` .
 This estimation can be computed with : func : ` xcorr ` or : func : ` CORRELATION ` by
 chosing a proprer lag ` L ` . The estimated PSD is then
 . . math : : \ hat { P } _ { xx } ( f ) = T \ sum _ { m = - L } ^ { L } \ hat { r } _ { xx } [ m ] exp ^ { - j2 \ pi fmT }
 The lag index must be less than the number of data samples ` N ` . Ideally , it
 should be around ` L / 10 ` [ Marple ] _ so as to avoid greater statistical
 variance associated with higher lags .
 To reduce the leakage of the implicit rectangular window and therefore to
 reduce the bias in the estimate , a tapering window is normally used and lead
 to the so - called Blackman and Tukey correlogram :
 . . math : : \ hat { P } _ { BT } ( f ) = T \ sum _ { m = - L } ^ { L } w [ m ] \ hat { r } _ { xx } [ m ] exp ^ { - j2 \ pi fmT }
 The correlogram for the cross power spectral estimate is
 . . math : : \ hat { P } _ { xx } ( f ) = T \ sum _ { m = - L } ^ { L } \ hat { r } _ { xx } [ m ] exp ^ { - j2 \ pi fmT }
 which is computed if : attr : ` Y ` is not provide . In such case ,
 : math : ` r _ { yx } = r _ { xy } ` so we compute the correlation only once .
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import CORRELOGRAMPSD , marple _ data
 from spectrum . tools import cshift
 from pylab import log10 , axis , grid , plot , linspace
 psd = CORRELOGRAMPSD ( marple _ data , marple _ data , lag = 15)
 f = linspace ( - 0.5 , 0.5 , len ( psd ) )
 psd = cshift ( psd , len ( psd ) / 2)
 plot ( f , 10 * log10 ( psd / max ( psd ) ) )
 axis ( [ - 0.5,0.5 , - 50,0 ] )
 grid ( True )
 . . seealso : : : func : ` create _ window ` , : func : ` CORRELATION ` , : func : ` xcorr ` ,
 : class : ` pcorrelogram ` ."""

implementation:
N = len ( X )
assert lag < N , 'lag must be < size of input data'
assert correlation_method in [ 'CORRELATION' , 'xcorr' ]
if Y is None :
    Y = numpy . array ( X )
    crosscorrelation = False
else :
    crosscorrelation = True
if NFFT is None :
    NFFT = N
psd = numpy . zeros ( NFFT , dtype = complex )
# Window should be centered around zero . Moreover , we want only the
# positive values . So , we need to use 2 * lag + 1 window and keep values on
# the right side .
w = Window ( 2. * lag + 1 , window , ** window_params )
w = w . data [ lag + 1 : ]
# compute the cross correlation
if correlation_method == 'CORRELATION' :
    rxy = CORRELATION ( X , Y , maxlags = lag , norm = norm )
elif correlation_method == 'xcorr' :
    rxy , _l = xcorr ( X , Y , maxlags = lag , norm = norm )
    rxy = rxy [ lag : ]
# keep track of the first elt .
psd [ 0 ] = rxy [ 0 ]
# create the first part of the PSD
psd [ 1 : lag + 1 ] = rxy [ 1 : ] * w
# create the second part .
# First , we need to compute the auto or cross correlation ryx
if crosscorrelation is True : # compute the cross correlation
    if correlation_method == 'CORRELATION' :
        ryx = CORRELATION ( Y , X , maxlags = lag , norm = norm )
    elif correlation_method == 'xcorr' :
        ryx , _l = xcorr ( Y , X , maxlags = lag , norm = norm )
        ryx = ryx [ lag : ]
    # print len ( ryx ) , len ( psd [ - 1 : NPSD - lag - 1 : - 1 ] )
    psd [ - 1 : NFFT - lag - 1 : - 1 ] = ryx [ 1 : ] . conjugate ( ) * w
else : # autocorrelation no additional correlation call required
    psd [ - 1 : NFFT - lag - 1 : - 1 ] = rxy [ 1 : ] . conjugate ( ) * w
psd = numpy . real ( fft ( psd ) )
return psd


signature:
def _get_data ( url ) :
    """Helper function to get data over http or from a local file"""

implementation:
if url . startswith ( 'http://' ) : # Try Python 2 , use Python 3 on exception
    try :
        resp = urllib . urlopen ( url )
        encoding = resp . headers . dict . get ( 'content-encoding' , 'plain' )
    except AttributeError :
        resp = urllib . request . urlopen ( url )
        encoding = resp . headers . get ( 'content-encoding' , 'plain' )
    data = resp . read ( )
    if encoding == 'plain' :
        pass
    elif encoding == 'gzip' :
        data = StringIO ( data )
        data = gzip . GzipFile ( fileobj = data ) . read ( )
    else :
        raise RuntimeError ( 'unknown encoding' )
else :
    with open ( url , 'r' ) as fid :
        data = fid . read ( )
return data


signature:
def _select_block ( str_in , start_tag , end_tag ) :
    """Select first block delimited by start _ tag and end _ tag"""

implementation:
start_pos = str_in . find ( start_tag )
if start_pos < 0 :
    raise ValueError ( 'start_tag not found' )
depth = 0
for pos in range ( start_pos , len ( str_in ) ) :
    if str_in [ pos ] == start_tag :
        depth += 1
    elif str_in [ pos ] == end_tag :
        depth -= 1
    if depth == 0 :
        break
sel = str_in [ start_pos + 1 : pos ]
return sel


signature:
def _parse_dict_recursive ( dict_str ) :
    """Parse a dictionary from the search index"""

implementation:
dict_out = dict ( )
pos_last = 0
pos = dict_str . find ( ':' )
while pos >= 0 :
    key = dict_str [ pos_last : pos ]
    if dict_str [ pos + 1 ] == '[' : # value is a list
        pos_tmp = dict_str . find ( ']' , pos + 1 )
        if pos_tmp < 0 :
            raise RuntimeError ( 'error when parsing dict' )
        value = dict_str [ pos + 2 : pos_tmp ] . split ( ',' )
        # try to convert elements to int
        for i in range ( len ( value ) ) :
            try :
                value [ i ] = int ( value [ i ] )
            except ValueError :
                pass
    elif dict_str [ pos + 1 ] == '{' : # value is another dictionary
        subdict_str = _select_block ( dict_str [ pos : ] , '{' , '}' )
        value = _parse_dict_recursive ( subdict_str )
        pos_tmp = pos + len ( subdict_str )
    else :
        raise ValueError ( 'error when parsing dict: unknown elem' )
    key = key . strip ( '"' )
    if len ( key ) > 0 :
        dict_out [ key ] = value
    pos_last = dict_str . find ( ',' , pos_tmp )
    if pos_last < 0 :
        break
    pos_last += 1
    pos = dict_str . find ( ':' , pos_last )
return dict_out


signature:
def parse_sphinx_searchindex ( searchindex ) :
    """Parse a Sphinx search index
 Parameters
 searchindex : str
 The Sphinx search index ( contents of searchindex . js )
 Returns
 filenames : list of str
 The file names parsed from the search index .
 objects : dict
 The objects parsed from the search index ."""

implementation:
# Make sure searchindex uses UTF - 8 encoding
if hasattr ( searchindex , 'decode' ) :
    searchindex = searchindex . decode ( 'UTF-8' )
# parse objects
query = 'objects:'
pos = searchindex . find ( query )
if pos < 0 :
    raise ValueError ( '"objects:" not found in search index' )
sel = _select_block ( searchindex [ pos : ] , '{' , '}' )
objects = _parse_dict_recursive ( sel )
# parse filenames
query = 'filenames:'
pos = searchindex . find ( query )
if pos < 0 :
    raise ValueError ( '"filenames:" not found in search index' )
filenames = searchindex [ pos + len ( query ) + 1 : ]
filenames = filenames [ : filenames . find ( ']' ) ]
filenames = [ f . strip ( '"' ) for f in filenames . split ( ',' ) ]
return filenames , objects


signature:
def embed_code_links ( app , exception ) :
    """Embed hyperlinks to documentation into example code"""

implementation:
if exception is not None :
    return
# No need to waste time embedding hyperlinks when not running the examples
# XXX : also at the time of writing this fixes make html - noplot
# for some reason I don ' t fully understand
if not app . builder . config . plot_gallery :
    return
# XXX : Whitelist of builders for which it makes sense to embed
# hyperlinks inside the example html . Note that the link embedding
# require searchindex . js to exist for the links to the local doc
# and there does not seem to be a good way of knowing which
# builders creates a searchindex . js .
if app . builder . name not in [ 'html' , 'readthedocs' ] :
    return
print ( 'Embedding documentation hyperlinks in examples..' )
gallery_conf = app . config . sphinx_gallery_conf
gallery_dirs = gallery_conf [ 'gallery_dirs' ]
if not isinstance ( gallery_dirs , list ) :
    gallery_dirs = [ gallery_dirs ]
for gallery_dir in gallery_dirs :
    _embed_code_links ( app , gallery_conf , gallery_dir )


signature:
def _get_link ( self , cobj ) :
    """Get a valid link , False if not found"""

implementation:
fname_idx = None
full_name = cobj [ 'module_short' ] + '.' + cobj [ 'name' ]
if full_name in self . _searchindex [ 'objects' ] :
    value = self . _searchindex [ 'objects' ] [ full_name ]
    if isinstance ( value , dict ) :
        value = value [ next ( iter ( value . keys ( ) ) ) ]
    fname_idx = value [ 0 ]
elif cobj [ 'module_short' ] in self . _searchindex [ 'objects' ] :
    value = self . _searchindex [ 'objects' ] [ cobj [ 'module_short' ] ]
    if cobj [ 'name' ] in value . keys ( ) :
        fname_idx = value [ cobj [ 'name' ] ] [ 0 ]
if fname_idx is not None :
    fname = self . _searchindex [ 'filenames' ] [ fname_idx ] + '.html'
    if self . _is_windows :
        fname = fname . replace ( '/' , '\\' )
        link = os . path . join ( self . doc_url , fname )
    else :
        link = posixpath . join ( self . doc_url , fname )
    if hasattr ( link , 'decode' ) :
        link = link . decode ( 'utf-8' , 'replace' )
    if link in self . _page_cache :
        html = self . _page_cache [ link ]
    else :
        html = get_data ( link , self . gallery_dir )
        self . _page_cache [ link ] = html
    # test if cobj appears in page
    comb_names = [ cobj [ 'module_short' ] + '.' + cobj [ 'name' ] ]
    if self . extra_modules_test is not None :
        for mod in self . extra_modules_test :
            comb_names . append ( mod + '.' + cobj [ 'name' ] )
    url = False
    if hasattr ( html , 'decode' ) : # Decode bytes under Python 3
        html = html . decode ( 'utf-8' , 'replace' )
    for comb_name in comb_names :
        if hasattr ( comb_name , 'decode' ) : # Decode bytes under Python 3
            comb_name = comb_name . decode ( 'utf-8' , 'replace' )
        if comb_name in html :
            url = link + u'#' + comb_name
    link = url
else :
    link = False
return link


signature:
def create_all_psd ( ) :
    f = pylab . linspace ( 0 , 1 , 4096 )
    pylab . figure ( figsize = ( 12 , 8 ) )
    # MA model
    p = spectrum . pma ( xx , 64 , 128 ) ;
    p ( ) ;
    p . plot ( )
    """# ARMA 15 order
 a , b , rho = spectrum . arma _ estimate ( data , 15,15 , 30)
 psd = spectrum . arma2psd ( A = a , B = b , rho = rho )
 newpsd = tools . cshift ( psd , len ( psd ) / / 2 ) # switch positive and negative freq
 pylab . plot ( f , 10 * pylab . log10 ( newpsd / max ( newpsd ) ) , label = ' ARMA 15,15 ' )"""

implementation:
# YULE WALKER
p = spectrum . pyule ( xx , 7 , NFFT = 4096 , scale_by_freq = False ) ;
p . plot ( )
# equivalent to
# plot ( [ x for x in p . frequencies ( ) ] , 10 * log10 ( p . psd ) ) ; grid ( True )
# burg method
p = spectrum . pburg ( xx , 7 , scale_by_freq = False ) ;
p . plot ( )
# pcovar
p = spectrum . pcovar ( xx , 7 , scale_by_freq = False ) ;
p . plot ( )
# pmodcovar
p = spectrum . pmodcovar ( xx , 7 , scale_by_freq = False ) ;
p . plot ( )
# correlogram
p = spectrum . pcorrelogram ( xx , lag = 60 , NFFT = 512 , scale_by_freq = False ) ;
p . plot ( )
# minvar
p = spectrum . pminvar ( xx , 7 , NFFT = 256 , scale_by_freq = False ) ;
p . plot ( )
# pmusic
p = spectrum . pmusic ( xx , 10 , 4 , scale_by_freq = False ) ;
p . plot ( )
# pmusic
p = spectrum . pev ( xx , 10 , 4 , scale_by_freq = False ) ;
p . plot ( )
# periodogram
p = spectrum . Periodogram ( xx , scale_by_freq = False ) ;
p . plot ( )
legend ( [ "MA 32" , "pyule 7" , "pburg 7" , "pcovar" , "pmodcovar" , "correlogram" , "minvar" , "pmusic" , "pev" , "periodgram" ] )
pylab . ylim ( [ - 80 , 80 ] )


signature:
def tf2zp ( b , a ) :
    """Convert transfer function filter parameters to zero - pole - gain form
 Find the zeros , poles , and gains of this continuous - time system :
 . . warning : : b and a must have the same length .
 from spectrum import tf2zp
 b = [ 2,3,0]
 a = [ 1 , 0.4 , 1]
 [ z , p , k ] = tf2zp ( b , a ) % Obtain zero - pole - gain form
 1.5
 -0.2000 + 0.9798i
 -0.2000 - 0.9798i
 : param b : numerator
 : param a : denominator
 : param fill : If True , check that the length of a and b are the same . If not , create a copy of the shortest element and append zeros to it .
 : return : z ( zeros ) , p ( poles ) , g ( gain )
 Convert transfer function f ( x ) = sum ( b * x ^ n ) / sum ( a * x ^ n ) to
 zero - pole - gain form f ( x ) = g * prod ( 1 - z * x ) / prod ( 1 - p * x )
 . . todo : : See if tf2ss followed by ss2zp gives better results . These
 are available from the control system toolbox . Note that
 the control systems toolbox doesn ' t bother , but instead uses
 . . seealso : : scipy . signal . tf2zpk , which gives the same results but uses a different
 algorithm ( z ^ - 1 instead of z ) ."""

implementation:
from numpy import roots
assert len ( b ) == len ( a ) , "length of the vectors a and b must be identical. fill with zeros if needed."
g = b [ 0 ] / a [ 0 ]
z = roots ( b )
p = roots ( a )
return z , p , g


signature:
def eqtflength ( b , a ) :
    """Given two list or arrays , pad with zeros the shortest array
 : param b : list or array
 : param a : list or array
 . . doctest : :
 > > > from spectrum . transfer import eqtflength
 > > > a = [ 1,2]
 > > > b = [ 1,2,3,4]
 > > > a , b , = eqtflength ( a , b )"""

implementation:
d = abs ( len ( b ) - len ( a ) )
if d != 0 :
    if len ( a ) > len ( b ) :
        try :
            b . extend ( [ 0. ] * d )
        except :
            b = np . append ( b , [ 0 ] * d )
    elif len ( b ) > len ( a ) :
        try :
            a . extend ( [ 0. ] * d )
        except :
            a = np . append ( a , [ 0 ] * d )
    return b , a
else :
    return b , a


signature:
def tf2zpk ( b , a ) :
    """Return zero , pole , gain ( z , p , k ) representation from a numerator ,
 denominator representation of a linear filter .
 Convert zero - pole - gain filter parameters to transfer function form
 : param ndarray b : numerator polynomial .
 : param ndarray a : numerator and denominator polynomials .
 : return :
 * z : ndarray Zeros of the transfer function .
 * p : ndarray Poles of the transfer function .
 * k : float System gain .
 If some values of b are too close to 0 , they are removed . In that case , a
 BadCoefficients warning is emitted .
 . . doctest : :
 > > > import scipy . signal
 > > > from spectrum . transfer import tf2zpk
 > > > [ b , a ] = scipy . signal . butter ( 3 . , . 4)
 > > > z , p , k = tf2zpk ( b , a )
 . . seealso : : : func : ` zpk2tf `
 . . note : : wrapper of scipy function tf2zpk"""

implementation:
import scipy . signal
z , p , k = scipy . signal . tf2zpk ( b , a )
return z , p , k


signature:
def ss2zpk ( a , b , c , d , input = 0 ) :
    """State - space representation to zero - pole - gain representation .
 : param A : ndarray State - space representation of linear system .
 : param B : ndarray State - space representation of linear system .
 : param C : ndarray State - space representation of linear system .
 : param D : ndarray State - space representation of linear system .
 : param int input : optional For multiple - input systems , the input to use .
 : return :
 * z , p : sequence Zeros and poles .
 * k : float System gain .
 . . note : : wrapper of scipy function ss2zpk"""

implementation:
import scipy . signal
z , p , k = scipy . signal . ss2zpk ( a , b , c , d , input = input )
return z , p , k


signature:
def zpk2tf ( z , p , k ) :
    r"""Return polynomial transfer function representation from zeros and poles
 : param ndarray z : Zeros of the transfer function .
 : param ndarray p : Poles of the transfer function .
 : param float k : System gain .
 : return :
 b : ndarray Numerator polynomial .
 a : ndarray Numerator and denominator polynomials .
 : func : ` zpk2tf ` forms transfer function polynomials from the zeros , poles , and gains
 of a system in factored form .
 zpk2tf ( z , p , k ) finds a rational transfer function
 . . math : : \ frac { B ( s ) } { A ( s ) } = \ frac { b _ 1 s ^ { n - 1 } + \ dots b _ { n - 1 } s + b _ n } { a _ 1 s ^ { m - 1 } + \ dots a _ { m - 1 } s + a _ m }
 given a system in factored transfer function form
 . . math : : H ( s ) = \ frac { Z ( s ) } { P ( s ) } = k \ frac { ( s - z _ 1 ) ( s - z _ 2 ) \ dots ( s - z _ m ) } { ( s - p _ 1 ) ( s - p _ 2 ) \ dots ( s - p _ n ) }
 with p being the pole locations , and z the zero locations , with as many .
 The gains for each numerator transfer function are in vector k .
 The zeros and poles must be real or come in complex conjugate pairs .
 The polynomial denominator coefficients are returned in row vector a and
 the polynomial numerator coefficients are returned in matrix b , which has
 as many rows as there are columns of z .
 Inf values can be used as place holders in z if some columns have fewer zeros than others .
 . . note : : wrapper of scipy function zpk2tf"""

implementation:
import scipy . signal
b , a = scipy . signal . zpk2tf ( z , p , k )
return b , a


signature:
def zpk2ss ( z , p , k ) :
    """Zero - pole - gain representation to state - space representation
 : param sequence z , p : Zeros and poles .
 : param float k : System gain .
 : return :
 * A , B , C , D : ndarray State - space matrices .
 . . note : : wrapper of scipy function zpk2ss"""

implementation:
import scipy . signal
return scipy . signal . zpk2ss ( z , p , k )


signature:
def create_window ( N , name = None , ** kargs ) :
    r"""Returns the N - point window given a valid name
 : param int N : window size
 : param str name : window name ( default is * rectangular * ) . Valid names
 are stored in : func : ` ~ spectrum . window . window _ names ` .
 : param kargs : optional arguments are :
 * * beta * : argument of the : func : ` window _ kaiser ` function ( default is 8.6)
 * * attenuation * : argument of the : func : ` window _ chebwin ` function ( default is 50dB )
 * * alpha * : argument of the
 1 . : func : ` window _ gaussian ` function ( default is 2.5)
 2 . : func : ` window _ blackman ` function ( default is 0.16)
 3 . : func : ` window _ poisson ` function ( default is 2)
 4 . : func : ` window _ cauchy ` function ( default is 3)
 * * mode * : argument : func : ` window _ flattop ` function ( default is * symmetric * , can be * periodic * )
 * * r * : argument of the : func : ` window _ tukey ` function ( default is 0.5 ) .
 The following windows have been simply wrapped from existing librairies like
 NumPy :
 * * * Rectangular * * : : func : ` window _ rectangle ` ,
 * * * Bartlett * * or Triangular : see : func : ` window _ bartlett ` ,
 * * * Hanning * * or Hann : see : func : ` window _ hann ` ,
 * * * Hamming * * : see : func : ` window _ hamming ` ,
 * * * Kaiser * * : see : func : ` window _ kaiser ` ,
 * * * chebwin * * : see : func : ` window _ chebwin ` .
 The following windows have been implemented from scratch :
 * * * Blackman * * : See : func : ` window _ blackman `
 * * * Bartlett - Hann * * : see : func : ` window _ bartlett _ hann `
 * * * cosine or sine * * : see : func : ` window _ cosine `
 * * * gaussian * * : see : func : ` window _ gaussian `
 * * * Bohman * * : see : func : ` window _ bohman `
 * * * Lanczos or sinc * * : see : func : ` window _ lanczos `
 * * * Blackman Harris * * : see : func : ` window _ blackman _ harris `
 * * * Blackman Nuttall * * : see : func : ` window _ blackman _ nuttall `
 * * * Nuttall * * : see : func : ` window _ nuttall `
 * * * Tukey * * : see : func : ` window _ tukey `
 * * * Parzen * * : see : func : ` window _ parzen `
 * * * Flattop * * : see : func : ` window _ flattop `
 * * * Riesz * * : see : func : ` window _ riesz `
 * * * Riemann * * : see : func : ` window _ riemann `
 * * * Poisson * * : see : func : ` window _ poisson `
 * * * Poisson - Hanning * * : see : func : ` window _ poisson _ hanning `
 . . todo : : on request taylor , potter , Bessel , expo ,
 rife - vincent , Kaiser - Bessel derived ( KBD )
 . . plot : :
 : width : 80%
 : include - source :
 from pylab import plot , legend
 from spectrum import create _ window
 data = create _ window ( 51 , ' hamming ' )
 plot ( data , label = ' hamming ' )
 data = create _ window ( 51 , ' kaiser ' )
 plot ( data , label = ' kaiser ' )
 legend ( )
 . . plot : :
 : width : 80%
 : include - source :
 from pylab import plot , log10 , linspace , fft , clip
 from spectrum import create _ window , fftshift
 A = fft ( create _ window ( 51 , ' hamming ' ) , 2048 ) / 25.5
 mag = abs ( fftshift ( A ) )
 freq = linspace ( - 0.5,0.5 , len ( A ) )
 response = 20 * log10 ( mag )
 mindB = - 60
 response = clip ( response , mindB , 100)
 plot ( freq , response )
 . . seealso : : : func : ` window _ visu ` , : func : ` Window ` , : mod : ` spectrum . dpss `"""

implementation:
if name is None :
    name = 'rectangle'
name = name . lower ( )
assert name in list ( window_names . keys ( ) ) , """window name %s not implemented or incorrect. Try to use one of %s""" % ( name , window_names )
# create the function name
f = eval ( window_names [ name ] )
windows_with_parameters = { 'kaiser' : { 'beta' : eval ( window_names [ 'kaiser' ] ) . __defaults__ [ 0 ] } , 'blackman' : { 'alpha' : eval ( window_names [ 'blackman' ] ) . __defaults__ [ 0 ] } , 'cauchy' : { 'alpha' : eval ( window_names [ 'cauchy' ] ) . __defaults__ [ 0 ] } , 'flattop' : { 'mode' : eval ( window_names [ 'flattop' ] ) . __defaults__ [ 0 ] } , 'gaussian' : { 'alpha' : eval ( window_names [ 'gaussian' ] ) . __defaults__ [ 0 ] } , 'chebwin' : { 'attenuation' : eval ( window_names [ 'chebwin' ] ) . __defaults__ [ 0 ] } , 'tukey' : { 'r' : eval ( window_names [ 'tukey' ] ) . __defaults__ [ 0 ] } , 'poisson' : { 'alpha' : eval ( window_names [ 'poisson' ] ) . __defaults__ [ 0 ] } , 'poisson_hanning' : { 'alpha' : eval ( window_names [ 'poisson_hanning' ] ) . __defaults__ [ 0 ] } , 'taylor' : { 'nbar' : eval ( window_names [ 'taylor' ] ) . __defaults__ [ 0 ] , 'sll' : eval ( window_names [ 'taylor' ] ) . __defaults__ [ 0 ] } , }
if name not in list ( windows_with_parameters . keys ( ) ) :
    if len ( kargs ) == 0 : # no parameters , so we directly call the function
        w = f ( N )
    else :
        raise ValueError ( """
            Parameters do not match any of the window. The window provided
            do not expect any parameters. Try to remove the parameters""" )
elif name in list ( windows_with_parameters . keys ( ) ) : # user optional parameters are provided , scan them :
    dargs = { }
    for arg in list ( kargs . keys ( ) ) : # check that the parameters are valid , and obtain the default value
        try :
            default = windows_with_parameters [ name ] [ arg ]
        except :
            raise ValueError ( """
                    Invalid optional argument (%s) for %s window.
                    Valid optional arguments are (%s)""" % ( arg , name , list ( windows_with_parameters [ name ] . keys ( ) ) ) )
        # add the user parameter to the list of parameters
        dargs [ arg ] = kargs . get ( arg , default )
    # call the proper function with the optional arguments
    w = f ( N , ** dargs )
return w


signature:
def enbw ( data ) :
    r"""Computes the equivalent noise bandwidth
 . . math : : ENBW = N \ frac { \ sum _ { n = 1 } ^ { N } w _ n ^ 2 } { \ left ( \ sum _ { n = 1 } ^ { N } w _ n \ right ) ^ 2}
 . . doctest : :
 > > > from spectrum import create _ window , enbw
 > > > w = create _ window ( 64 , ' rectangular ' )
 > > > enbw ( w )
 1.0
 The following table contains the ENBW values for some of the
 implemented windows in this module ( with N = 16384 ) . They have been
 double checked against litterature ( Source : [ Harris ] _ , [ Marple ] _ ) .
 If not present , it means that it has not been checked .
 name ENBW litterature
 rectangular 1 . 1.
 triangle 1.3334 1.33
 Hann 1.5001 1.5
 Hamming 1.3629 1.36
 blackman 1.7268 1.73
 kaiser 1.7
 blackmanharris , 4 2.004 2.
 riesz 1.2000 1.2
 riemann 1.32 1.3
 parzen 1.917 1.92
 tukey 0.25 1.102 1.1
 bohman 1.7858 1.79
 poisson 2 1.3130 1.3
 hanningpoisson 0.5 1.609 1.61
 cauchy 1.489 1.48
 lanczos 1.3"""

implementation:
N = len ( data )
return N * np . sum ( data ** 2 ) / np . sum ( data ) ** 2


signature:
def _kaiser ( n , beta ) :
    """Independant Kaiser window
 For the definition of the Kaiser window , see A . V . Oppenheim & R . W . Schafer , " Discrete - Time Signal Processing " .
 The continuous version of width n centered about x = 0 is :
 . . note : : 2 times slower than scipy . kaiser"""

implementation:
from scipy . special import iv as besselI
m = n - 1
k = arange ( 0 , m )
k = 2. * beta / m * sqrt ( k * ( m - k ) )
w = besselI ( 0 , k ) / besselI ( 0 , beta )
return w


signature:
def window_visu ( N = 51 , name = 'hamming' , ** kargs ) :
    """A Window visualisation tool
 : param N : length of the window
 : param name : name of the window
 : param NFFT : padding used by the FFT
 : param mindB : the minimum frequency power in dB
 : param maxdB : the maximum frequency power in dB
 : param kargs : optional arguments passed to : func : ` create _ window `
 This function plot the window shape and its equivalent in the Fourier domain .
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' kaiser ' , beta = 8 . )"""

implementation:
# get the default parameters
mindB = kargs . pop ( 'mindB' , - 100 )
maxdB = kargs . pop ( 'maxdB' , None )
norm = kargs . pop ( 'norm' , True )
# create a window object
w = Window ( N , name , ** kargs )
# plot the time and frequency windows
w . plot_time_freq ( mindB = mindB , maxdB = maxdB , norm = norm )


signature:
def window_kaiser ( N , beta = 8.6 , method = 'numpy' ) :
    r"""Kaiser window
 : param N : window length
 : param beta : kaiser parameter ( default is 8.6)
 To obtain a Kaiser window that designs an FIR filter with
 sidelobe attenuation of : math : ` \ alpha ` dB , use the following : math : ` \ beta ` where
 : math : ` \ beta = \ pi \ alpha ` .
 . . math : :
 w _ n = \ frac { I _ 0 \ left ( \ pi \ alpha \ sqrt { 1 - \ left ( \ frac { 2n } { M } - 1 \ right ) ^ 2 } \ right ) } { I _ 0 ( \ pi \ alpha ) }
 where
 * : math : ` I _ 0 ` is the zeroth order Modified Bessel function of the first kind .
 * : math : ` \ alpha ` is a real number that determines the shape of the
 window . It determines the trade - off between main - lobe width and side
 lobe level .
 * the length of the sequence is N = M + 1.
 The Kaiser window can approximate many other windows by varying
 the : math : ` \ beta ` parameter :
 beta Window shape
 0 Rectangular
 5 Similar to a Hamming
 6 Similar to a Hanning
 8.6 Similar to a Blackman
 . . plot : :
 : width : 80%
 : include - source :
 from pylab import plot , legend , xlim
 from spectrum import window _ kaiser
 N = 64
 for beta in [ 1,2,4,8,16 ] :
 plot ( window _ kaiser ( N , beta ) , label = ' beta = ' + str ( beta ) )
 xlim ( 0 , N )
 legend ( )
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' kaiser ' , beta = 8 . )
 . . seealso : : numpy . kaiser , : func : ` spectrum . window . create _ window `"""

implementation:
if N == 1 :
    return ones ( 1 )
if method == 'numpy' :
    from numpy import kaiser
    return kaiser ( N , beta )
else :
    return _kaiser ( N , beta )


signature:
def window_blackman ( N , alpha = 0.16 ) :
    r"""Blackman window
 : param N : window length
 . . math : : a _ 0 - a _ 1 \ cos ( \ frac { 2 \ pi n } { N - 1 } ) + a _ 2 \ cos ( \ frac { 4 \ pi n } { N - 1 } )
 with
 . . math : :
 a _ 0 = ( 1 - \ alpha ) / 2 , a _ 1 = 0.5 , a _ 2 = \ alpha / 2 \ rm { \ ; and \ ; \ alpha } = 0.16
 When : math : ` \ alpha = 0.16 ` , this is the unqualified Blackman window with
 : math : ` a _ 0 = 0.48 ` and : math : ` a _ 2 = 0.08 ` .
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' blackman ' )
 . . note : : Although Numpy implements a blackman window for : math : ` \ alpha = 0.16 ` ,
 this implementation is valid for any : math : ` \ alpha ` .
 . . seealso : : numpy . blackman , : func : ` create _ window ` , : class : ` Window `"""

implementation:
a0 = ( 1. - alpha ) / 2.
a1 = 0.5
a2 = alpha / 2.
if ( N == 1 ) :
    win = array ( [ 1. ] )
else :
    k = arange ( 0 , N ) / float ( N - 1. )
    win = a0 - a1 * cos ( 2 * pi * k ) + a2 * cos ( 4 * pi * k )
return win


signature:
def window_gaussian ( N , alpha = 2.5 ) :
    r"""Gaussian window
 : param N : window length
 . . math : : \ exp ^ { - 0.5 \ left ( \ sigma \ frac { n } { N / 2 } \ right ) ^ 2}
 with : math : ` \ frac { N - 1 } { 2 } \ leq n \ leq \ frac { N - 1 } { 2 } ` .
 . . note : : N - 1 is used to be in agreement with octave convention . The ENBW of
 1.4 is also in agreement with [ Harris ] _
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' gaussian ' , alpha = 2.5)
 . . seealso : : scipy . signal . gaussian , : func : ` create _ window `"""


signature:
def window_chebwin ( N , attenuation = 50 ) :
    """Cheb window
 : param N : window length
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' chebwin ' , attenuation = 50)
 . . seealso : : scipy . signal . chebwin , : func : ` create _ window ` , : class : ` Window `"""

implementation:
import scipy . signal
return scipy . signal . chebwin ( N , attenuation )


signature:
def window_cosine ( N ) :
    r"""Cosine tapering window also known as sine window .
 : param N : window length
 . . math : : w ( n ) = \ cos \ left ( \ frac { \ pi n } { N - 1 } - \ frac { \ pi } { 2 } \ right ) = \ sin \ left ( \ frac { \ pi n } { N - 1 } \ right )
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' cosine ' )
 . . seealso : : : func : ` create _ window ` , : class : ` Window `"""

implementation:
if N == 1 :
    return ones ( 1 )
n = arange ( 0 , N )
win = sin ( pi * n / ( N - 1. ) )
return win


signature:
def window_lanczos ( N ) :
    r"""Lanczos window also known as sinc window .
 : param N : window length
 . . math : : w ( n ) = sinc \ left ( \ frac { 2n } { N - 1 } - 1 \ right )
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' lanczos ' )
 . . seealso : : : func : ` create _ window ` , : class : ` Window `"""

implementation:
if N == 1 :
    return ones ( 1 )
n = linspace ( - N / 2. , N / 2. , N )
win = sinc ( 2 * n / ( N - 1. ) )
return win


signature:
def window_bartlett_hann ( N ) :
    r"""Bartlett - Hann window
 : param N : window length
 . . math : : w ( n ) = a _ 0 + a _ 1 \ left | \ frac { n } { N - 1 } - \ frac { 1 } { 2 } \ right | - a _ 2 \ cos \ left ( \ frac { 2 \ pi n } { N - 1 } \ right )
 with : math : ` a _ 0 = 0.62 ` , : math : ` a _ 1 = 0.48 ` and : math : ` a _ 2 = 0.38 `
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' bartlett _ hann ' )
 . . seealso : : : func : ` create _ window ` , : class : ` Window `"""


signature:
def _coeff4 ( N , a0 , a1 , a2 , a3 ) :
    """a common internal function to some window functions with 4 coeffs
 For the blackmna harris for instance , the results are identical to octave if N is odd
 but not for even values . . . if n = 0 whatever N is , the w ( 0 ) must be equal to a0 - a1 + a2 - a3 , which
 is the case here , but not in octave . . ."""

implementation:
if N == 1 :
    return ones ( 1 )
n = arange ( 0 , N )
N1 = N - 1.
w = a0 - a1 * cos ( 2. * pi * n / N1 ) + a2 * cos ( 4. * pi * n / N1 ) - a3 * cos ( 6. * pi * n / N1 )
return w


signature:
def window_nuttall ( N ) :
    r"""Nuttall tapering window
 : param N : window length
 . . math : : w ( n ) = a _ 0 - a _ 1 \ cos \ left ( \ frac { 2 \ pi n } { N - 1 } \ right ) + a _ 2 \ cos \ left ( \ frac { 4 \ pi n } { N - 1 } \ right ) - a _ 3 \ cos \ left ( \ frac { 6 \ pi n } { N - 1 } \ right )
 with : math : ` a _ 0 = 0.355768 ` , : math : ` a _ 1 = 0.487396 ` , : math : ` a _ 2 = 0.144232 ` and : math : ` a _ 3 = 0.012604 `
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' nuttall ' , mindB = - 80)
 . . seealso : : : func : ` create _ window ` , : class : ` Window `"""


signature:
def window_blackman_nuttall ( N ) :
    r"""Blackman Nuttall window
 returns a minimum , 4 - term Blackman - Harris window . The window is minimum in the sense that its maximum sidelobes are minimized .
 The coefficients for this window differ from the Blackman - Harris window coefficients and produce slightly lower sidelobes .
 : param N : window length
 . . math : : w ( n ) = a _ 0 - a _ 1 \ cos \ left ( \ frac { 2 \ pi n } { N - 1 } \ right ) + a _ 2 \ cos \ left ( \ frac { 4 \ pi n } { N - 1 } \ right ) - a _ 3 \ cos \ left ( \ frac { 6 \ pi n } { N - 1 } \ right )
 with : math : ` a _ 0 = 0.3635819 ` , : math : ` a _ 1 = 0.4891775 ` , : math : ` a _ 2 = 0.1365995 ` and : math : ` 0_3 = . 0106411 `
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' blackman _ nuttall ' , mindB = - 80)
 . . seealso : : : func : ` spectrum . window . create _ window `
 . . seealso : : : func : ` create _ window ` , : class : ` Window `"""


signature:
def window_blackman_harris ( N ) :
    r"""Blackman Harris window
 : param N : window length
 . . math : : w ( n ) = a _ 0 - a _ 1 \ cos \ left ( \ frac { 2 \ pi n } { N - 1 } \ right ) + a _ 2 \ cos \ left ( \ frac { 4 \ pi n } { N - 1 } \ right ) - a _ 3 \ cos \ left ( \ frac { 6 \ pi n } { N - 1 } \ right )
 coeff value
 : math : ` a _ 0 ` 0.35875
 : math : ` a _ 1 ` 0.48829
 : math : ` a _ 2 ` 0.14128
 : math : ` a _ 3 ` 0.01168
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' blackman _ harris ' , mindB = - 80)
 . . seealso : : : func : ` spectrum . window . create _ window `
 . . seealso : : : func : ` create _ window ` , : class : ` Window `"""


signature:
def window_bohman ( N ) :
    r"""Bohman tapering window
 : param N : window length
 . . math : : w ( n ) = ( 1 - | x | ) \ cos ( \ pi | x | ) + \ frac { 1 } { \ pi } \ sin ( \ pi | x | )
 where x is a length N vector of linearly spaced values between
 -1 and 1.
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' bohman ' )
 . . seealso : : : func : ` create _ window ` , : class : ` Window `"""

implementation:
x = linspace ( - 1 , 1 , N )
w = ( 1. - abs ( x ) ) * cos ( pi * abs ( x ) ) + 1. / pi * sin ( pi * abs ( x ) )
return w


signature:
def window_tukey ( N , r = 0.5 ) :
    """Tukey tapering window ( or cosine - tapered window )
 : param N : window length
 : param r : defines the ratio between the constant section and the cosine
 section . It has to be between 0 and 1.
 The function returns a Hanning window for ` r = 0 ` and a full box for ` r = 1 ` .
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' tukey ' )
 window _ visu ( 64 , ' tukey ' , r = 1)
 . . math : : 0.5 ( 1 + cos ( 2pi / r ( x - r / 2 ) ) ) for 0 < = x < r / 2
 . . math : : 0.5 ( 1 + cos ( 2pi / r ( x - 1 + r / 2 ) ) ) for x > = 1 - r / 2
 . . seealso : : : func : ` create _ window ` , : class : ` Window `"""

implementation:
assert r >= 0 and r <= 1 , "r must be in [0,1]"
if N == 1 :
    return ones ( 1 )
if r == 0 :
    return ones ( N )
elif r == 1 :
    return window_hann ( N )
else :
    from numpy import flipud , concatenate , where
    # # cosine - tapered window
    x = linspace ( 0 , 1 , N )
    x1 = where ( x < r / 2. )
    w = 0.5 * ( 1 + cos ( 2 * pi / r * ( x [ x1 [ 0 ] ] - r / 2 ) ) )
    w = concatenate ( ( w , ones ( N - len ( w ) * 2 ) , flipud ( w ) ) )
    return w


signature:
def window_parzen ( N ) :
    r"""Parsen tapering window ( also known as de la Valle - Poussin )
 : param N : window length
 Parzen windows are piecewise cubic approximations
 of Gaussian windows . Parzen window sidelobes fall off as : math : ` 1 / \ omega ^ 4 ` .
 if : math : ` 0 \ leq | x | \ leq ( N - 1 ) / 4 ` :
 . . math : : w ( n ) = 1-6 \ left ( \ frac { | n | } { N / 2 } \ right ) ^ 2 + 6 \ left ( \ frac { | n | } { N / 2 } \ right ) ^ 3
 if : math : ` ( N - 1 ) / 4 \ leq | x | \ leq ( N - 1 ) / 2 `
 . . math : : w ( n ) = 2 \ left ( 1 - \ frac { | n | } { N / 2 } \ right ) ^ 3
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' parzen ' )
 . . seealso : : : func : ` create _ window ` , : class : ` Window `"""


signature:
def window_flattop ( N , mode = 'symmetric' , precision = None ) :
    r"""Flat - top tapering window
 Returns symmetric or periodic flat top window .
 : param N : window length
 : param mode : way the data are normalised . If mode is * symmetric * , then
 divide n by N - 1 . IF mode is * periodic * , divide by N ,
 to be consistent with octave code .
 When using windows for filter design , the * symmetric * mode
 should be used ( default ) . When using windows for spectral analysis , the * periodic *
 mode should be used . The mathematical form of the flat - top window in the symmetric
 case is :
 . . math : : w ( n ) = a _ 0
 - a _ 1 \ cos \ left ( \ frac { 2 \ pi n } { N - 1 } \ right )
 + a _ 2 \ cos \ left ( \ frac { 4 \ pi n } { N - 1 } \ right )
 - a _ 3 \ cos \ left ( \ frac { 6 \ pi n } { N - 1 } \ right )
 + a _ 4 \ cos \ left ( \ frac { 8 \ pi n } { N - 1 } \ right )
 coeff value
 a0 0.21557895
 a1 0.41663158
 a2 0.277263158
 a3 0.083578947
 a4 0.006947368
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' bohman ' )
 . . seealso : : : func : ` create _ window ` , : class : ` Window `"""

implementation:
assert mode in [ 'periodic' , 'symmetric' ]
t = arange ( 0 , N )
# FIXME : N = 1 for mode = periodic ?
if mode == 'periodic' :
    x = 2 * pi * t / float ( N )
else :
    if N == 1 :
        return ones ( 1 )
    x = 2 * pi * t / float ( N - 1 )
a0 = 0.21557895
a1 = 0.41663158
a2 = 0.277263158
a3 = 0.083578947
a4 = 0.006947368
if precision == 'octave' : # to compare with octave , same as above but less precise
    d = 4.6402
    a0 = 1. / d
    a1 = 1.93 / d
    a2 = 1.29 / d
    a3 = 0.388 / d
    a4 = 0.0322 / d
w = a0 - a1 * cos ( x ) + a2 * cos ( 2 * x ) - a3 * cos ( 3 * x ) + a4 * cos ( 4 * x )
return w


signature:
def window_taylor ( N , nbar = 4 , sll = - 30 ) :
    """Taylor tapering window
 Taylor windows allows you to make tradeoffs between the
 mainlobe width and sidelobe level ( sll ) .
 Implemented as described by Carrara , Goodman , and Majewski
 in ' Spotlight Synthetic Aperture Radar : Signal Processing Algorithms '
 Pages 512-513
 : param N : window length
 : param float nbar :
 : param float sll :
 The default values gives equal height
 sidelobes ( nbar ) and maximum sidelobe level ( sll ) .
 . . warning : : not implemented
 . . seealso : : : func : ` create _ window ` , : class : ` Window `"""

implementation:
B = 10 ** ( - sll / 20 )
A = log ( B + sqrt ( B ** 2 - 1 ) ) / pi
s2 = nbar ** 2 / ( A ** 2 + ( nbar - 0.5 ) ** 2 )
ma = arange ( 1 , nbar )
def calc_Fm ( m ) :
    numer = ( - 1 ) ** ( m + 1 ) * prod ( 1 - m ** 2 / s2 / ( A ** 2 + ( ma - 0.5 ) ** 2 ) )
    denom = 2 * prod ( [ 1 - m ** 2 / j ** 2 for j in ma if j != m ] )
    return numer / denom
Fm = array ( [ calc_Fm ( m ) for m in ma ] )
def W ( n ) :
    return 2 * np . sum ( Fm * cos ( 2 * pi * ma * ( n - N / 2 + 1 / 2 ) / N ) ) + 1
w = array ( [ W ( n ) for n in range ( N ) ] )
# normalize ( Note that this is not described in the original text )
scale = W ( ( N - 1 ) / 2 )
w /= scale
return w


signature:
def window_riesz ( N ) :
    r"""Riesz tapering window
 : param N : window length
 . . math : : w ( n ) = 1 - \ left | \ frac { n } { N / 2 } \ right | ^ 2
 with : math : ` - N / 2 \ leq n \ leq N / 2 ` .
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' riesz ' )
 . . seealso : : : func : ` create _ window ` , : class : ` Window `"""


signature:
def window_riemann ( N ) :
    r"""Riemann tapering window
 : param int N : window length
 . . math : : w ( n ) = 1 - \ left | \ frac { n } { N / 2 } \ right | ^ 2
 with : math : ` - N / 2 \ leq n \ leq N / 2 ` .
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' riesz ' )
 . . seealso : : : func : ` create _ window ` , : class : ` Window `"""

implementation:
n = linspace ( - N / 2. , ( N ) / 2. , N )
w = sin ( n / float ( N ) * 2. * pi ) / ( n / float ( N ) * 2. * pi )
return w


signature:
def window_poisson ( N , alpha = 2 ) :
    r"""Poisson tapering window
 : param int N : window length
 . . math : : w ( n ) = \ exp ^ { - \ alpha \ frac { | n | } { N / 2 } }
 with : math : ` - N / 2 \ leq n \ leq N / 2 ` .
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' poisson ' )
 window _ visu ( 64 , ' poisson ' , alpha = 3)
 window _ visu ( 64 , ' poisson ' , alpha = 4)
 . . seealso : : : func : ` create _ window ` , : class : ` Window `"""

implementation:
n = linspace ( - N / 2. , ( N ) / 2. , N )
w = exp ( - alpha * abs ( n ) / ( N / 2. ) )
return w


signature:
def window_poisson_hanning ( N , alpha = 2 ) :
    r"""Hann - Poisson tapering window
 This window is constructed as the product of the Hanning and Poisson
 windows . The parameter * * alpha * * is the Poisson parameter .
 : param int N : window length
 : param float alpha : parameter of the poisson window
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' poisson _ hanning ' , alpha = 0.5)
 window _ visu ( 64 , ' poisson _ hanning ' , alpha = 1)
 window _ visu ( 64 , ' poisson _ hanning ' )
 . . seealso : : : func : ` window _ poisson ` , : func : ` window _ hann `"""

implementation:
w1 = window_hann ( N )
w2 = window_poisson ( N , alpha = alpha )
return w1 * w2


signature:
def window_cauchy ( N , alpha = 3 ) :
    r"""Cauchy tapering window
 : param int N : window length
 : param float alpha : parameter of the poisson window
 . . math : : w ( n ) = \ frac { 1 } { 1 + \ left ( \ frac { \ alpha * n } { N / 2 } \ right ) * * 2}
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import window _ visu
 window _ visu ( 64 , ' cauchy ' , alpha = 3)
 window _ visu ( 64 , ' cauchy ' , alpha = 4)
 window _ visu ( 64 , ' cauchy ' , alpha = 5)
 . . seealso : : : func : ` window _ poisson ` , : func : ` window _ hann `"""


signature:
def compute_response ( self , ** kargs ) :
    """Compute the window data frequency response
 : param norm : True by default . normalised the frequency data .
 : param int NFFT : total length of the final data sets ( 2048 by default .
 if less than data length , then NFFT is set to the data length * 2 ) .
 The response is stored in : attr : ` response ` .
 . . note : : Units are dB ( 20 log10 ) since we plot the frequency response )"""

implementation:
from numpy . fft import fft , fftshift
norm = kargs . get ( 'norm' , self . norm )
# do some padding . Default is max ( 2048 , data . len * 2)
NFFT = kargs . get ( 'NFFT' , 2048 )
if NFFT < len ( self . data ) :
    NFFT = self . data . size * 2
# compute the fft modulus
A = fft ( self . data , NFFT )
mag = abs ( fftshift ( A ) )
# do we want to normalise the data
if norm is True :
    mag = mag / max ( mag )
response = 20. * stools . log10 ( mag )
# factor 20 we are looking at the response
# not the powe
# response = clip ( response , mindB , 100)
self . __response = response


signature:
def plot_frequencies ( self , mindB = None , maxdB = None , norm = True ) :
    """Plot the window in the frequency domain
 : param mindB : change the default lower y bound
 : param maxdB : change the default upper lower bound
 : param bool norm : if True , normalise the frequency response .
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum . window import Window
 w = Window ( 64 , name = ' hamming ' )
 w . plot _ frequencies ( )"""

implementation:
from pylab import plot , title , xlim , grid , ylim , xlabel , ylabel
# recompute the response
self . compute_response ( norm = norm )
plot ( self . frequencies , self . response )
title ( "ENBW=%2.1f" % ( self . enbw ) )
ylabel ( 'Frequency response (dB)' )
xlabel ( 'Fraction of sampling frequency' )
# define the plot limits
xlim ( - 0.5 , 0.5 )
y0 , y1 = ylim ( )
if mindB :
    y0 = mindB
if maxdB is not None :
    y1 = maxdB
else :
    y1 = max ( self . response )
ylim ( y0 , y1 )
grid ( True )


signature:
def plot_window ( self ) :
    """Plot the window in the time domain
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum . window import Window
 w = Window ( 64 , name = ' hamming ' )
 w . plot _ window ( )"""

implementation:
from pylab import plot , xlim , grid , title , ylabel , axis
x = linspace ( 0 , 1 , self . N )
xlim ( 0 , 1 )
plot ( x , self . data )
grid ( True )
title ( '%s Window (%s points)' % ( self . name . capitalize ( ) , self . N ) )
ylabel ( 'Amplitude' )
axis ( [ 0 , 1 , 0 , 1.1 ] )


signature:
def plot_time_freq ( self , mindB = - 100 , maxdB = None , norm = True , yaxis_label_position = "right" ) :
    """Plotting method to plot both time and frequency domain results .
 See : meth : ` plot _ frequencies ` for the optional arguments .
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum . window import Window
 w = Window ( 64 , name = ' hamming ' )
 w . plot _ time _ freq ( )"""

implementation:
from pylab import subplot , gca
subplot ( 1 , 2 , 1 )
self . plot_window ( )
subplot ( 1 , 2 , 2 )
self . plot_frequencies ( mindB = mindB , maxdB = maxdB , norm = norm )
if yaxis_label_position == "left" :
    try :
        tight_layout ( )
    except :
        pass
else :
    ax = gca ( )
    ax . yaxis . set_label_position ( "right" )


signature:
def TOEPLITZ ( T0 , TC , TR , Z ) :
    """solve the general toeplitz linear equations
 Solve TX = Z
 : param T0 : zero lag value
 : param TC : r1 to rN
 : param TR : r1 to rN
 returns X
 requires 3M ^ 2 + M operations instead of M ^ 3 with gaussian elimination
 . . warning : : not used right now"""

implementation:
assert len ( TC ) > 0
assert len ( TC ) == len ( TR )
M = len ( TC )
X = numpy . zeros ( M + 1 , dtype = complex )
A = numpy . zeros ( M , dtype = complex )
B = numpy . zeros ( M , dtype = complex )
P = T0
if P == 0 :
    raise ValueError ( "P must be different from zero" )
if P == 0 :
    raise ValueError ( "P must be different from zero" )
X [ 0 ] = Z [ 0 ] / T0
for k in range ( 0 , M ) :
    save1 = TC [ k ]
    save2 = TR [ k ]
    beta = X [ 0 ] * TC [ k ]
    if k == 0 :
        temp1 = - save1 / P
        temp2 = - save2 / P
    else :
        for j in range ( 0 , k ) :
            save1 = save1 + A [ j ] * TC [ k - j - 1 ]
            save2 = save2 + B [ j ] * TR [ k - j - 1 ]
            beta = beta + X [ j + 1 ] * TC [ k - j - 1 ]
        temp1 = - save1 / P
        temp2 = - save2 / P
    P = P * ( 1. - ( temp1 * temp2 ) )
    if P <= 0 :
        raise ValueError ( "singular matrix" )
    A [ k ] = temp1
    B [ k ] = temp2
    alpha = ( Z [ k + 1 ] - beta ) / P
    if k == 0 :
        X [ k + 1 ] = alpha
        for j in range ( 0 , k + 1 ) :
            X [ j ] = X [ j ] + alpha * B [ k - j ]
        continue
    for j in range ( 0 , k ) :
        kj = k - j - 1
        save1 = A [ j ]
        A [ j ] = save1 + temp1 * B [ kj ]
        B [ kj ] = B [ kj ] + temp2 * save1
    X [ k + 1 ] = alpha
    for j in range ( 0 , k + 1 ) :
        X [ j ] = X [ j ] + alpha * B [ k - j ]
return X


signature:
def HERMTOEP ( T0 , T , Z ) :
    """solve Tx = Z by a variation of Levinson algorithm where T
 is a complex hermitian toeplitz matrix
 : param T0 : zero lag value
 : param T : r1 to rN
 : return : X
 used by eigen PSD method"""

implementation:
assert len ( T ) > 0
M = len ( T )
X = numpy . zeros ( M + 1 , dtype = complex )
A = numpy . zeros ( M , dtype = complex )
P = T0
if P == 0 :
    raise ValueError ( "P must be different from zero" )
X [ 0 ] = Z [ 0 ] / T0
for k in range ( 0 , M ) :
    save = T [ k ]
    beta = X [ 0 ] * T [ k ]
    if k == 0 :
        temp = - save / P
    else :
        for j in range ( 0 , k ) :
            save = save + A [ j ] * T [ k - j - 1 ]
            beta = beta + X [ j + 1 ] * T [ k - j - 1 ]
        temp = - save / P
    P = P * ( 1. - ( temp . real ** 2 + temp . imag ** 2 ) )
    if P <= 0 :
        raise ValueError ( "singular matrix" )
    A [ k ] = temp
    alpha = ( Z [ k + 1 ] - beta ) / P
    if k == 0 : # print ' skipping code for k = 0'
        X [ k + 1 ] = alpha
        for j in range ( 0 , k + 1 ) :
            X [ j ] = X [ j ] + alpha * A [ k - j ] . conjugate ( )
        continue
    khalf = ( k + 1 ) // 2
    for j in range ( 0 , khalf ) :
        kj = k - j - 1
        save = A [ j ]
        A [ j ] = save + temp * A [ kj ] . conjugate ( )
        if j != kj :
            A [ kj ] = A [ kj ] + temp * save . conjugate ( )
    X [ k + 1 ] = alpha
    for j in range ( 0 , k + 1 ) :
        X [ j ] = X [ j ] + alpha * A [ k - j ] . conjugate ( )
return X


signature:
def get_short_module_name ( module_name , obj_name ) :
    """Get the shortest possible module name"""

implementation:
parts = module_name . split ( '.' )
short_name = module_name
for i in range ( len ( parts ) - 1 , 0 , - 1 ) :
    short_name = '.' . join ( parts [ : i ] )
    try :
        exec ( 'from %s import %s' % ( short_name , obj_name ) )
    except ImportError : # get the last working module name
        short_name = '.' . join ( parts [ : ( i + 1 ) ] )
        break
return short_name


signature:
def identify_names ( code ) :
    """Builds a codeobj summary by identifying and resolving used names
 > > > code = ' ' '
 . . . from a . b import c
 . . . import d as e
 . . . print ( c )
 . . . e . HelloWorld ( ) . f . g
 > > > for name , o in sorted ( identify _ names ( code ) . items ( ) ) :
 . . . print ( name , o [ ' name ' ] , o [ ' module ' ] , o [ ' module _ short ' ] )
 c c a . b a . b
 e . HelloWorld HelloWorld d d"""

implementation:
finder = NameFinder ( )
finder . visit ( ast . parse ( code ) )
example_code_obj = { }
for name , full_name in finder . get_mapping ( ) : # name is as written in file ( e . g . np . asarray )
# full _ name includes resolved import path ( e . g . numpy . asarray )
    module , attribute = full_name . rsplit ( '.' , 1 )
    # get shortened module name
    module_short = get_short_module_name ( module , attribute )
    cobj = { 'name' : attribute , 'module' : module , 'module_short' : module_short }
    example_code_obj [ name ] = cobj
return example_code_obj


signature:
def _thumbnail_div ( full_dir , fname , snippet , is_backref = False ) :
    """Generates RST to place a thumbnail in a gallery"""

implementation:
thumb = os . path . join ( full_dir , 'images' , 'thumb' , 'sphx_glr_%s_thumb.png' % fname [ : - 3 ] )
ref_name = os . path . join ( full_dir , fname ) . replace ( os . path . sep , '_' )
template = BACKREF_THUMBNAIL_TEMPLATE if is_backref else THUMBNAIL_TEMPLATE
return template . format ( snippet = snippet , thumbnail = thumb , ref_name = ref_name )


signature:
def modcovar_marple ( X , IP ) :
    """Fast algorithm for the solution of the modified covariance least squares normal equations .
 This implementation is based on [ Marple ] _ . This code is far more
 complicated and slower than : func : ` modcovar ` function , which is now the official version .
 See : func : ` modcovar ` for a detailed description of Modified Covariance method .
 : param X : - Array of complex data samples X ( 1 ) through X ( N )
 : param int IP : - Order of linear prediction model ( integer )
 : return :
 * P - Real linear prediction variance at order IP
 * A - Array of complex linear prediction coefficients
 * ISTAT - Integer status indicator at time of exit
 0 . for normal exit ( no numerical ill - conditioning )
 1 . if P is not a positive value
 2 . if DELTA ' and GAMMA ' do not lie in the range 0 to 1
 3 . if P ' is not a positive value
 4 . if DELTA and GAMMA do not lie in the range 0 to 1
 : validation : the AR parameters are the same as those returned by
 a completely different function : func : ` modcovar ` .
 . . note : : validation . results similar to test example in Marple but
 starts to differ for ip ~ 8 . with ratio of 0.975 for ip = 15 probably
 due to precision .
 : References : [ Marple ] _"""

implementation:
Pv = [ ]
N = len ( X )
A = np . zeros ( N , dtype = complex )
D = np . zeros ( N , dtype = complex )
C = np . zeros ( N , dtype = complex )
R = np . zeros ( N , dtype = complex )
# Initialization
R1 = 0.
for K in range ( 1 , N - 1 ) :
    R1 = R1 + 2. * ( X [ K ] . real ** 2 + X [ K ] . imag ** 2 )
R2 = X [ 0 ] . real ** 2 + X [ 0 ] . imag ** 2
R3 = X [ N - 1 ] . real ** 2 + X [ N - 1 ] . imag ** 2
R4 = 1. / ( R1 + 2. * ( R2 + R3 ) )
P = R1 + R2 + R3
DELTA = 1. - R2 * R4
GAMMA = 1. - R3 * R4
LAMBDA = ( X [ 0 ] * X [ N - 1 ] ) . conjugate ( ) * R4
C [ 0 ] = X [ N - 1 ] * R4
D [ 0 ] = X [ 0 ] . conjugate ( ) * R4
M = 0
if ( IP == 0 ) :
    P = ( .5 * R1 + R2 + R3 ) / float ( N )
    return [ ] , P , [ ]
# C Main loop
for M in range ( 0 , IP ) : # print ' - - - - - ' , M
    SAVE1 = 0 + 0j
    for K in range ( M + 1 , N ) :
        SAVE1 = SAVE1 + X [ K ] * X [ K - M - 1 ] . conjugate ( )
    SAVE1 *= 2.
    R [ M ] = SAVE1 . conjugate ( )
    THETA = X [ N - 1 ] * D [ 0 ]
    PSI = X [ N - 1 ] * C [ 0 ]
    XI = X [ 0 ] . conjugate ( ) * D [ 0 ]
    if M == 0 :
        pass
    else :
        for K in range ( 0 , M ) :
            THETA = THETA + X [ N - K - 2 ] * D [ K + 1 ]
            # Eq . ( 8 . D . 45)
            PSI = PSI + X [ N - K - 2 ] * C [ K + 1 ]
            # Eq . ( 8 . D . 45)
            XI = XI + X [ K + 1 ] . conjugate ( ) * D [ K + 1 ]
            # Eq . ( 8 . D . 45)
            R [ K ] = R [ K ] - X [ N - M - 1 ] * X [ N + 1 - M + K - 1 ] . conjugate ( ) - X [ M ] . conjugate ( ) * X [ M - K - 1 ]
            # Eq . ( 8 . D . 37)
            SAVE1 = SAVE1 + R [ K ] . conjugate ( ) * A [ M - K - 1 ]
        # Eq . ( 8 . D . 24)
        # print ' withini loop ' , K , THETA , PSI , XI , R [ K ] , SAVE1
        # print ' x - - - - - ' , C [ K ] , C [ K + 1]
        # print ' x - - - - - ' , D [ K ] , D [ K + 1]
        # print ' x - - - - - ' , N - K - 2 , K + 1 , N - M - 1 , M , M - K - 1
        # print ' x - - - - - ' , X [ N - K - 2 ] , X [ K + 1 ] , X [ N - M - 1 ] , X [ M ] , X [ M - K - 1]
    # Order update of A vector
    C1 = - SAVE1 / P
    A [ M ] = C1
    # Eq . ( 8 . D . 23)
    P = P * ( 1. - C1 . real ** 2 - C1 . imag ** 2 )
    # Eq . ( 8 . D . 25)
    # print ' C1 = ' , C1 , ' A [ M ] ' , A [ M ] , ' P = ' , P
    if M == 0 :
        pass
    else : # for k in range ( 0 , ( M - 1 ) / 2 + 1 ) :
    # print ' BEFORE ' , A [ 0:4]
        for K in range ( 0 , ( M + 1 ) // 2 ) :
            MK = M - K - 1
            SAVE1 = A [ K ]
            A [ K ] = SAVE1 + C1 * A [ MK ] . conjugate ( )
            # Eq . ( 8 . D . 22)
            # print ' A uopdate - - - - > ' , M , K , MK , A [ K ] , A [ MK ] , SAVE1 , C1
            if ( K != MK ) : # print ' ENTERI IF ' , K , MK , A [ MK ] , C1 , SAVE1
                A [ MK ] = A [ MK ] + C1 * ( SAVE1 . conjugate ( ) )
                # Eq . ( 8 . D . 22)
        # print ' AFTER ' , A [ 0:4]
    # print M + 1 , IP
    if M + 1 == IP : # Pv . append ( P )
        P = .5 * P / float ( N - M - 1 )
        Pv . append ( P )
        return A , P , Pv
    else :
        Pv . append ( .5 * P / float ( N - M - 1 ) )
        # Pv . append ( P )
    # Time update of C , D vectors and GAMMA , DELTA , LAMBDA scalars
    R1 = 1. / ( DELTA * GAMMA - ( LAMBDA . real ) ** 2 - ( LAMBDA . imag ) ** 2 )
    C1 = ( THETA * ( LAMBDA . conjugate ( ) ) + PSI * DELTA ) * R1
    C2 = ( PSI * LAMBDA + THETA * GAMMA ) * R1
    C3 = ( XI * ( LAMBDA . conjugate ( ) ) + THETA * DELTA ) * R1
    C4 = ( THETA * LAMBDA + XI * GAMMA ) * R1
    # print ' C1C2C3C4 ' , C1 , C2 , C3 , C4
    # print M , M / 2 + 1 , ( M - 1 ) / 2 + 1
    for K in range ( 0 , ( M ) // 2 + 1 ) :
        MK = M - K
        SAVE1 = C [ K ] . conjugate ( )
        SAVE2 = D [ K ] . conjugate ( )
        SAVE3 = C [ MK ] . conjugate ( )
        SAVE4 = D [ MK ] . conjugate ( )
        C [ K ] = C [ K ] + C1 * SAVE3 + C2 * SAVE4
        # Eq . ( 8 . D . 43)
        D [ K ] = D [ K ] + C3 * SAVE3 + C4 * SAVE4
        # Eq . ( 8 . D . 44)
        # print ' - - - - - > ' , K , MK , SAVE1 , SAVE2 , SAVE3 , SAVE4
        if K != MK :
            C [ MK ] = C [ MK ] + C1 * SAVE1 + C2 * SAVE2
            # Eq . ( 8 . D . 43)
            D [ MK ] = D [ MK ] + C3 * SAVE1 + C4 * SAVE2
            # Eq . ( 8 . D . 44)
            # print ' witninloop loop ' , K , MK , C [ MK ] , D [ MK ]
    # print ' after loop '
    # print ' C = ' , C [ 0:2]
    # print ' D = ' , D [ 0:2]
    R2 = PSI . real ** 2 + PSI . imag ** 2
    R3 = THETA . real ** 2 + THETA . imag ** 2
    R4 = XI . real ** 2 + XI . imag ** 2
    R5 = GAMMA - ( R2 * DELTA + R3 * GAMMA + 2. * np . real ( PSI * LAMBDA * THETA . conjugate ( ) ) ) * R1
    R2 = DELTA - ( R3 * DELTA + R4 * GAMMA + 2. * np . real ( THETA * LAMBDA * XI . conjugate ( ) ) ) * R1
    GAMMA = R5
    # Eq . ( 8 . D . 46)
    DELTA = R2
    # Eq . ( 8 . D . 47)
    LAMBDA = LAMBDA + C3 * PSI . conjugate ( ) + C4 * THETA . conjugate ( )
    # Eq . ( 8 . D . 48)
    # print ' before time updaqt ' , R2 , R3 , R4 , R5 , LAMBDA
    if P <= 0. :
        raise ValueError ( 'Found an invalid negative value for P.' )
    if ( DELTA > 0. and DELTA <= 1. and GAMMA > 0. and GAMMA <= 1. ) :
        pass
    else :
        raise ValueError ( 'Found an invalid DELTA or GAMMA value.' )
    # Time update of A vector ; order updates of C , D vectors aand GAMMA ,
    # DELTA , LAMBDA scalars
    # print ' time update - - - - - '
    R1 = 1. / P
    R2 = 1. / ( DELTA * GAMMA - LAMBDA . real ** 2 - LAMBDA . imag ** 2 )
    # Eq . ( 8 . D . 41)
    EF = X [ M + 1 ]
    EB = X [ N - M - 2 ]
    # print ' efeb ' , EF , EB
    # print ' A = ' , A [ 0:2]
    for K in range ( 0 , M + 1 ) : # print ' efeb lloop ' , K , A [ K ] , X [ M - K ] , X [ N - M + K - 1]
        EF = EF + A [ K ] * X [ M - K ]
        # Eq . ( 8 . D . 1)
        EB = EB + A [ K ] . conjugate ( ) * X [ N - M + K - 1 ]
        # Eq . ( 8 . D . 2)
    C1 = EB * R1
    # Eq . ( 8 . D . 28)
    C2 = EF . conjugate ( ) * R1
    # Eq . ( 8 . D . 29)
    C3 = ( EB . conjugate ( ) * DELTA + EF * LAMBDA ) * R2
    C4 = ( EF * GAMMA + ( EB * LAMBDA ) . conjugate ( ) ) * R2
    # print ' c1c2c3c4 ' , C1 , C2 , C3 , C4
    # print ' efeb ' , EF , EB
    for K in range ( M , - 1 , - 1 ) :
        SAVE1 = A [ K ]
        A [ K ] = SAVE1 + C3 * C [ K ] + C4 * D [ K ]
        # Eq . ( 8 . D . 38)
        C [ K + 1 ] = C [ K ] + C1 * SAVE1
        # Eq . ( 8 . D . 26)
        D [ K + 1 ] = D [ K ] + C2 * SAVE1
        # Eq . ( 8 . D . 27)
        # print ' ACD ' , K , A [ K ] , C [ K + 1 ] , D [ K + 1]
    C [ 0 ] = C1
    D [ 0 ] = C2
    R3 = EB . real ** 2 + EB . imag ** 2
    R4 = EF . real ** 2 + EF . imag ** 2
    P = P - ( R3 * DELTA + R4 * GAMMA + 2. * np . real ( EF * EB * LAMBDA ) ) * R2
    # Eq . ( 8 . D . 42)
    DELTA = DELTA - R4 * R1
    # Eq . ( 8 . D . 32)
    GAMMA = GAMMA - R3 * R1
    # Eq . ( 8 . D . 33)
    LAMBDA = LAMBDA + ( EF * EB ) . conjugate ( ) * R1
    # Eq . ( 8 . D . 35)
    # print ' final '
    # print R3 , R4 , P , DELTA , GAMMA , LAMBDA
    # print ' Afinal = ' , A [ 0:2]
    # print ' P = ' , P
    if ( P > 0. ) :
        pass
    else :
        raise ValueError ( "Found a invalid negative P value " )
    if ( DELTA > 0. and DELTA <= 1. and GAMMA > 0. and GAMMA <= 1. ) :
        pass
    else :
        raise ValueError ( "Found an invalid negative GAMMA or DELTA value " )


signature:
def modcovar ( x , order ) :
    """Simple and fast implementation of the covariance AR estimate
 This code is 10 times faster than : func : ` modcovar _ marple ` and more importantly
 only 10 lines of code , compared to a 200 loc for : func : ` modcovar _ marple `
 : param X : Array of complex data samples
 : param int order : Order of linear prediction model
 : return :
 * P - Real linear prediction variance at order IP
 * A - Array of complex linear prediction coefficients
 . . plot : :
 : include - source :
 : width : 80%
 from spectrum import modcovar , marple _ data , arma2psd , cshift
 from pylab import log10 , linspace , axis , plot
 a , p = modcovar ( marple _ data , 15)
 PSD = arma2psd ( a )
 PSD = cshift ( PSD , len ( PSD ) / 2 ) # switch positive and negative freq
 plot ( linspace ( - 0.5 , 0.5 , 4096 ) , 10 * log10 ( PSD / max ( PSD ) ) )
 axis ( [ - 0.5,0.5 , - 60,0 ] )
 . . seealso : : : class : ` ~ spectrum . modcovar . pmodcovar `
 : validation : the AR parameters are the same as those returned by
 a completely different function : func : ` modcovar _ marple ` .
 : References : Mathworks"""

implementation:
from spectrum import corrmtx
import scipy . linalg
X = corrmtx ( x , order , 'modified' )
Xc = np . matrix ( X [ : , 1 : ] )
X1 = np . array ( X [ : , 0 ] )
# Coefficients estimated via the covariance method
# Here we use lstsq rathre than solve function because Xc is not square matrix
a , residues , rank , singular_values = scipy . linalg . lstsq ( - Xc , X1 )
# Estimate the input white noise variance
Cz = np . dot ( X1 . conj ( ) . transpose ( ) , Xc )
e = np . dot ( X1 . conj ( ) . transpose ( ) , X1 ) + np . dot ( Cz , a )
assert e . imag < 1e-4 , 'wierd behaviour'
e = float ( e . real )
# ignore imag part that should be small
return a , e


signature:
def aryule ( X , order , norm = 'biased' , allow_singularity = True ) :
    r"""Compute AR coefficients using Yule - Walker method
 : param X : Array of complex data values , X ( 1 ) to X ( N )
 : param int order : Order of autoregressive process to be fitted ( integer )
 : param str norm : Use a biased or unbiased correlation .
 : param bool allow _ singularity :
 : return :
 * AR coefficients ( complex )
 * variance of white noise ( Real )
 * reflection coefficients for use in lattice filter
 . . rubric : : Description :
 The Yule - Walker method returns the polynomial A corresponding to the
 AR parametric signal model estimate of vector X using the Yule - Walker
 ( autocorrelation ) method . The autocorrelation may be computed using a
 * * biased * * or * * unbiased * * estimation . In practice , the biased estimate of
 the autocorrelation is used for the unknown true autocorrelation . Indeed ,
 an unbiased estimate may result in nonpositive - definite autocorrelation
 matrix .
 So , a biased estimate leads to a stable AR filter .
 The following matrix form represents the Yule - Walker equations . The are
 solved by means of the Levinson - Durbin recursion :
 . . math : :
 \ left ( \ begin { array } { cccc }
 r ( 1 ) & r ( 2 ) ^ * & \ dots & r ( n ) ^ * \ \
 r ( 2 ) & r ( 1 ) ^ * & \ dots & r ( n - 1 ) ^ * \ \
 \ dots & \ dots & \ dots & \ dots \ \
 r ( n ) & \ dots & r ( 2 ) & r ( 1 ) \ end { array } \ right )
 \ left ( \ begin { array } { cccc }
 a ( 2 ) \ \
 a ( 3 ) \ \
 \ dots \ \
 a ( n + 1 ) \ end { array } \ right )
 \ left ( \ begin { array } { cccc }
 - r ( 2 ) \ \
 - r ( 3 ) \ \
 \ dots \ \
 - r ( n + 1 ) \ end { array } \ right )
 The outputs consists of the AR coefficients , the estimated variance of the
 white noise process , and the reflection coefficients . These outputs can be
 used to estimate the optimal order by using : mod : ` ~ spectrum . criteria ` .
 . . rubric : : Examples :
 From a known AR process or order 4 , we estimate those AR parameters using
 the aryule function .
 . . doctest : :
 > > > from scipy . signal import lfilter
 > > > from spectrum import *
 > > > from numpy . random import randn
 > > > A = [ 1 , - 2.7607 , 3.8106 , - 2.6535 , 0.9238]
 > > > noise = randn ( 1 , 1024)
 > > > y = lfilter ( [ 1 ] , A , noise ) ;
 > > > # filter a white noise input to create AR ( 4 ) process
 > > > [ ar , var , reflec ] = aryule ( y [ 0 ] , 4)
 > > > # ar should contains values similar to A
 The PSD estimate of a data samples is computed and plotted as follows :
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import *
 from pylab import *
 ar , P , k = aryule ( marple _ data , 15 , norm = ' biased ' )
 psd = arma2psd ( ar )
 plot ( linspace ( - 0.5 , 0.5 , 4096 ) , 10 * log10 ( psd / max ( psd ) ) )
 axis ( [ - 0.5 , 0.5 , - 60 , 0 ] )
 . . note : : The outputs have been double checked against ( 1 ) octave outputs
 ( octave has norm = ' biased ' by default ) and ( 2 ) Marple test code .
 . . seealso : : This function uses : func : ` ~ spectrum . levinson . LEVINSON ` and
 : func : ` ~ spectrum . correlation . CORRELATION ` . See the : mod : ` ~ spectrum . criteria `
 module for criteria to automatically select the AR order .
 : References : [ Marple ] _"""

implementation:
assert norm in [ 'biased' , 'unbiased' ]
r = CORRELATION ( X , maxlags = order , norm = norm )
A , P , k = LEVINSON ( r , allow_singularity = allow_singularity )
return A , P , k


signature:
def LEVINSON ( r , order = None , allow_singularity = False ) :
    r"""Levinson - Durbin recursion .
 Find the coefficients of a length ( r ) - 1 order autoregressive linear process
 : param r : autocorrelation sequence of length N + 1 ( first element being the zero - lag autocorrelation )
 : param order : requested order of the autoregressive coefficients . default is N .
 : param allow _ singularity : false by default . Other implementations may be True ( e . g . , octave )
 : return :
 * the ` N + 1 ` autoregressive coefficients : math : ` A = ( 1 , a _ 1 . . . a _ N ) `
 * the prediction errors
 * the ` N ` reflections coefficients values
 This algorithm solves the set of complex linear simultaneous equations
 using Levinson algorithm .
 . . math : :
 \ bold { T } _ M \ left ( \ begin { array } { c } 1 \ \ \ bold { a } _ M \ end { array } \ right ) =
 \ left ( \ begin { array } { c } \ rho _ M \ \ \ bold { 0 } _ M \ end { array } \ right )
 where : math : ` \ bold { T } _ M ` is a Hermitian Toeplitz matrix with elements
 : math : ` T _ 0 , T _ 1 , \ dots , T _ M ` .
 . . note : : Solving this equations by Gaussian elimination would
 require : math : ` M ^ 3 ` operations whereas the levinson algorithm
 requires : math : ` M ^ 2 + M ` additions and : math : ` M ^ 2 + M ` multiplications .
 This is equivalent to solve the following symmetric Toeplitz system of
 linear equations
 . . math : :
 \ left ( \ begin { array } { cccc }
 r _ 1 & r _ 2 ^ * & \ dots & r _ { n } ^ * \ \
 r _ 2 & r _ 1 ^ * & \ dots & r _ { n - 1 } ^ * \ \
 \ dots & \ dots & \ dots & \ dots \ \
 r _ n & \ dots & r _ 2 & r _ 1 \ end { array } \ right )
 \ left ( \ begin { array } { cccc }
 a _ 2 \ \
 a _ 3 \ \
 \ dots \ \
 a _ { N + 1 } \ end { array } \ right )
 \ left ( \ begin { array } { cccc }
 - r _ 2 \ \
 - r _ 3 \ \
 \ dots \ \
 - r _ { N + 1 } \ end { array } \ right )
 where : math : ` r = ( r _ 1 . . . r _ { N + 1 } ) ` is the input autocorrelation vector , and
 : math : ` r _ i ^ * ` denotes the complex conjugate of : math : ` r _ i ` . The input r is typically
 a vector of autocorrelation coefficients where lag 0 is the first
 element : math : ` r _ 1 ` .
 . . doctest : :
 > > > import numpy ; from spectrum import LEVINSON
 > > > T = numpy . array ( [ 3 . , - 2 + 0.5j , . 7-1j ] )
 > > > a , e , k = LEVINSON ( T )"""

implementation:
# from numpy import isrealobj
T0 = numpy . real ( r [ 0 ] )
T = r [ 1 : ]
M = len ( T )
if order is None :
    M = len ( T )
else :
    assert order <= M , 'order must be less than size of the input data'
    M = order
realdata = numpy . isrealobj ( r )
if realdata is True :
    A = numpy . zeros ( M , dtype = float )
    ref = numpy . zeros ( M , dtype = float )
else :
    A = numpy . zeros ( M , dtype = complex )
    ref = numpy . zeros ( M , dtype = complex )
P = T0
for k in range ( 0 , M ) :
    save = T [ k ]
    if k == 0 :
        temp = - save / P
    else : # save + = sum ( [ A [ j ] * T [ k - j - 1 ] for j in range ( 0 , k ) ] )
        for j in range ( 0 , k ) :
            save = save + A [ j ] * T [ k - j - 1 ]
        temp = - save / P
    if realdata :
        P = P * ( 1. - temp ** 2. )
    else :
        P = P * ( 1. - ( temp . real ** 2 + temp . imag ** 2 ) )
    if P <= 0 and allow_singularity == False :
        raise ValueError ( "singular matrix" )
    A [ k ] = temp
    ref [ k ] = temp
    # save reflection coeff at each step
    if k == 0 :
        continue
    khalf = ( k + 1 ) // 2
    if realdata is True :
        for j in range ( 0 , khalf ) :
            kj = k - j - 1
            save = A [ j ]
            A [ j ] = save + temp * A [ kj ]
            if j != kj :
                A [ kj ] += temp * save
    else :
        for j in range ( 0 , khalf ) :
            kj = k - j - 1
            save = A [ j ]
            A [ j ] = save + temp * A [ kj ] . conjugate ( )
            if j != kj :
                A [ kj ] = A [ kj ] + temp * save . conjugate ( )
return A , P , ref


signature:
def rlevinson ( a , efinal ) :
    """computes the autocorrelation coefficients , R based
 on the prediction polynomial A and the final prediction error Efinal ,
 using the stepdown algorithm .
 Works for real or complex data
 : param a :
 : param efinal :
 : return :
 * R , the autocorrelation
 * U prediction coefficient
 * kr reflection coefficients
 * e errors
 A should be a minimum phase polynomial and A ( 1 ) is assumed to be unity .
 : returns : ( P + 1 ) by ( P + 1 ) upper triangular matrix , U ,
 that holds the i ' th order prediction polynomials
 Ai , i = 1 : P , where P is the order of the input
 polynomial , A .
 [ 1 a1(1 ) * a2(2 ) * . . . . . aP ( P ) * ]
 [ 0 1 a2(1 ) * . . . . . aP ( P - 1 ) * ]
 [ 0 0 0 . . . . . 1 ]
 from which the i ' th order prediction polynomial can be extracted
 using Ai = U ( i + 1 : - 1:1 , i + 1 ) ' . The first row of U contains the
 conjugates of the reflection coefficients , and the K ' s may be
 extracted using , K = conj ( U ( 1,2 : end ) ) .
 . . todo : : remove the conjugate when data is real data , clean up the code
 test and doc ."""

implementation:
a = numpy . array ( a )
realdata = numpy . isrealobj ( a )
assert a [ 0 ] == 1 , 'First coefficient of the prediction polynomial must be unity'
p = len ( a )
if p < 2 :
    raise ValueError ( 'Polynomial should have at least two coefficients' )
if realdata == True :
    U = numpy . zeros ( ( p , p ) )
    # This matrix will have the prediction
    # polynomials of orders 1 : p
else :
    U = numpy . zeros ( ( p , p ) , dtype = complex )
U [ : , p - 1 ] = numpy . conj ( a [ - 1 : : - 1 ] )
# Prediction coefficients of order p
p = p - 1
e = numpy . zeros ( p )
# First we find the prediction coefficients of smaller orders and form the
# Matrix U
# Initialize the step down
e [ - 1 ] = efinal
# Prediction error of order p
# Step down
for k in range ( p - 1 , 0 , - 1 ) :
    [ a , e [ k - 1 ] ] = levdown ( a , e [ k ] )
    U [ : , k ] = numpy . concatenate ( ( numpy . conj ( a [ - 1 : : - 1 ] . transpose ( ) ) , [ 0 ] * ( p - k ) ) )
e0 = e [ 0 ] / ( 1. - abs ( a [ 1 ] ** 2 ) )
# % Because a [ 1 ] = 1 ( true polynomial )
U [ 0 , 0 ] = 1
# % Prediction coefficient of zeroth order
kr = numpy . conj ( U [ 0 , 1 : ] )
# % The reflection coefficients
kr = kr . transpose ( )
# % To make it into a column vector
# % Once we have the matrix U and the prediction error at various orders , we can
# % use this information to find the autocorrelation coefficients .
R = numpy . zeros ( 1 , dtype = complex )
# % Initialize recursion
k = 1
R0 = e0
# To take care of the zero indexing problem
R [ 0 ] = - numpy . conj ( U [ 0 , 1 ] ) * R0
# R [ 1 ] = - a1[1 ] * R [ 0]
# Actual recursion
for k in range ( 1 , p ) :
    r = - sum ( numpy . conj ( U [ k - 1 : : - 1 , k ] ) * R [ - 1 : : - 1 ] ) - kr [ k ] * e [ k - 1 ]
    R = numpy . insert ( R , len ( R ) , r )
# Include R ( 0 ) and make it a column vector . Note the dot transpose
# R = [ R0 R ] . ' ;
R = numpy . insert ( R , 0 , e0 )
return R , U , kr , e


signature:
def levdown ( anxt , enxt = None ) :
    """One step backward Levinson recursion
 : param anxt :
 : param enxt :
 : return :
 * acur the P ' th order prediction polynomial based on the P + 1 ' th order prediction polynomial , anxt .
 * ecur the the P ' th order prediction error based on the P + 1 ' th order prediction error , enxt .
 . . * knxt the P + 1 ' th order reflection coefficient ."""

implementation:
# % Some preliminaries first
# if nargout > = 2 & nargin < 2
# raise ValueError ( ' Insufficient number of input arguments ' ) ;
if anxt [ 0 ] != 1 :
    raise ValueError ( 'At least one of the reflection coefficients is equal to one.' )
anxt = anxt [ 1 : ]
# Drop the leading 1 , it is not needed
# in the step down
# Extract the k + 1 ' th reflection coefficient
knxt = anxt [ - 1 ]
if knxt == 1.0 :
    raise ValueError ( 'At least one of the reflection coefficients is equal to one.' )
# A Matrix formulation from Stoica is used to avoid looping
acur = ( anxt [ 0 : - 1 ] - knxt * numpy . conj ( anxt [ - 2 : : - 1 ] ) ) / ( 1. - abs ( knxt ) ** 2 )
ecur = None
if enxt is not None :
    ecur = enxt / ( 1. - numpy . dot ( knxt . conj ( ) . transpose ( ) , knxt ) )
acur = numpy . insert ( acur , 0 , 1 )
return acur , ecur


signature:
def levup ( acur , knxt , ecur = None ) :
    """LEVUP One step forward Levinson recursion
 : param acur :
 : param knxt :
 : return :
 * anxt the P + 1 ' th order prediction polynomial based on the P ' th order prediction polynomial , acur , and the
 P + 1 ' th order reflection coefficient , Knxt .
 * enxt the P + 1 ' th order prediction prediction error , based on the P ' th order prediction error , ecur .
 : References : P . Stoica R . Moses , Introduction to Spectral Analysis Prentice Hall , N . J . , 1997 , Chapter 3."""

implementation:
if acur [ 0 ] != 1 :
    raise ValueError ( 'At least one of the reflection coefficients is equal to one.' )
acur = acur [ 1 : ]
# Drop the leading 1 , it is not needed
# Matrix formulation from Stoica is used to avoid looping
anxt = numpy . concatenate ( ( acur , [ 0 ] ) ) + knxt * numpy . concatenate ( ( numpy . conj ( acur [ - 1 : : - 1 ] ) , [ 1 ] ) )
enxt = None
if ecur is not None : # matlab version enxt = ( 1 - knxt ' . * knxt ) * ecur
    enxt = ( 1. - numpy . dot ( numpy . conj ( knxt ) , knxt ) ) * ecur
anxt = numpy . insert ( anxt , 0 , 1 )
return anxt , enxt


signature:
def arcovar ( x , order ) :
    r"""Simple and fast implementation of the covariance AR estimate
 This code is 10 times faster than : func : ` arcovar _ marple ` and more importantly
 only 10 lines of code , compared to a 200 loc for : func : ` arcovar _ marple `
 : param array X : Array of complex data samples
 : param int oder : Order of linear prediction model
 : return :
 * a - Array of complex forward linear prediction coefficients
 * e - error
 The covariance method fits a Pth order autoregressive ( AR ) model to the
 input signal , which is assumed to be the output of
 an AR system driven by white noise . This method minimizes the forward
 prediction error in the least - squares sense . The output vector
 contains the normalized estimate of the AR system parameters
 The white noise input variance estimate is also returned .
 If is the power spectral density of y ( n ) , then :
 . . math : : \ frac { e } { \ left | A ( e ^ { jw } ) \ right | ^ 2 } = \ frac { e } { \ left | 1 + \ sum _ { k - 1 } ^ P a ( k ) e ^ { - jwk } \ right | ^ 2}
 Because the method characterizes the input data using an all - pole model ,
 the correct choice of the model order p is important .
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import arcovar , marple _ data , arma2psd
 from pylab import plot , log10 , linspace , axis
 ar _ values , error = arcovar ( marple _ data , 15)
 psd = arma2psd ( ar _ values , sides = ' centerdc ' )
 plot ( linspace ( - 0.5 , 0.5 , len ( psd ) ) , 10 * log10 ( psd / max ( psd ) ) )
 axis ( [ - 0.5 , 0.5 , - 60 , 0 ] )
 . . seealso : : : class : ` pcovar `
 : validation : the AR parameters are the same as those returned by
 a completely different function : func : ` arcovar _ marple ` .
 : References : [ Mathworks ] _"""

implementation:
from spectrum import corrmtx
import scipy . linalg
X = corrmtx ( x , order , 'covariance' )
Xc = np . matrix ( X [ : , 1 : ] )
X1 = np . array ( X [ : , 0 ] )
# Coefficients estimated via the covariance method
# Here we use lstsq rathre than solve function because Xc is not square
# matrix
a , _residues , _rank , _singular_values = scipy . linalg . lstsq ( - Xc , X1 )
# Estimate the input white noise variance
Cz = np . dot ( X1 . conj ( ) . transpose ( ) , Xc )
e = np . dot ( X1 . conj ( ) . transpose ( ) , X1 ) + np . dot ( Cz , a )
assert e . imag < 1e-4 , 'wierd behaviour'
e = float ( e . real )
# ignore imag part that should be small
return a , e


signature:
def lpc ( x , N = None ) :
    """Linear Predictor Coefficients .
 : param x :
 : param int N : default is length ( X ) - 1
 : Details :
 Finds the coefficients : math : ` A = ( 1 , a ( 2 ) , \ dots a ( N + 1 ) ) ` , of an Nth order
 forward linear predictor that predicts the current value value of the
 real - valued time series x based on past samples :
 . . math : : \ hat { x } ( n ) = - a ( 2 ) * x ( n - 1 ) - a ( 3 ) * x ( n - 2 ) - . . . - a ( N + 1 ) * x ( n - N )
 such that the sum of the squares of the errors
 . . math : : err ( n ) = X ( n ) - Xp ( n )
 is minimized . This function uses the Levinson - Durbin recursion to
 solve the normal equations that arise from the least - squares formulation .
 . . seealso : : : func : ` levinson ` , : func : ` aryule ` , : func : ` prony ` , : func : ` stmcb `
 . . todo : : matrix case , references
 : Example :
 from scipy . signal import lfilter
 noise = randn ( 50000,1 ) ; % Normalized white Gaussian noise
 x = filter ( [ 1 ] , [ 1 1/2 1/3 1/4 ] , noise )
 x = x [ 45904:50000]
 x . reshape ( 4096 , 1)
 x = x [ 0]
 Compute the predictor coefficients , estimated signal , prediction error , and autocorrelation sequence of the prediction error :
 1.00000 + 0.00000i 0.51711 - 0.00000i 0.33908 - 0.00000i 0.24410 - 0.00000i
 a = lpc ( x , 3)
 est _ x = lfilter ( [ 0 - a ( 2 : end ) ] , 1 , x ) ; % Estimated signal
 e = x - est _ x ; % Prediction error
 [ acs , lags ] = xcorr ( e , ' coeff ' ) ; % ACS of prediction error"""

implementation:
m = len ( x )
if N is None :
    N = m - 1
    # default value if N is not provided
elif N > m - 1 : # disp ( ' Warning : zero - padding short input sequence ' )
    x . resize ( N + 1 )
    # todo : check this zero - padding .
X = fft ( x , 2 ** nextpow2 ( 2. * len ( x ) - 1 ) )
R = real ( ifft ( abs ( X ) ** 2 ) )
R = R / ( m - 1. )
# Biased autocorrelation estimate
a , e , ref = LEVINSON ( R , N )
return a , e


implementation:
errors . is_positive_integer ( n )
result = numpy . zeros ( ( n , n ) )
# fill the first row and column
for i in range ( 0 , n ) :
    result [ i , 0 ] = 1
    result [ 0 , i ] = 1
if n > 1 :
    for i in range ( 1 , n ) :
        for j in range ( 1 , n ) :
            result [ i , j ] = result [ i - 1 , j ] + result [ i , j - 1 ]
return result