signature:
def p_moduleIdentityClause ( self , p ) :
    """moduleIdentityClause : LOWERCASE _ IDENTIFIER MODULE _ IDENTITY SubjectCategoriesPart LAST _ UPDATED ExtUTCTime ORGANIZATION Text CONTACT _ INFO Text DESCRIPTION Text RevisionPart COLON _ COLON _ EQUAL ' { ' objectIdentifier ' } '"""

implementation:
p [ 0 ] = ( 'moduleIdentityClause' , p [ 1 ] , # id
# p [ 2 ] , # MODULE _ IDENTITY
# XXX p [ 3 ] , # SubjectCategoriesPart
( p [ 4 ] , p [ 5 ] ) , # last updated
( p [ 6 ] , p [ 7 ] ) , # organization
( p [ 8 ] , p [ 9 ] ) , # contact info
( p [ 10 ] , p [ 11 ] ) , # description
p [ 12 ] , # RevisionPart
p [ 15 ] )


signature:
def p_ObjectSyntax ( self , p ) :
    """ObjectSyntax : SimpleSyntax
 | conceptualTable
 | row
 | entryType
 | ApplicationSyntax
 | typeTag SimpleSyntax"""


signature:
def p_SimpleSyntax ( self , p ) :
    """SimpleSyntax : INTEGER
 | INTEGER integerSubType
 | INTEGER enumSpec
 | INTEGER32
 | INTEGER32 integerSubType
 | UPPERCASE _ IDENTIFIER enumSpec
 | UPPERCASE _ IDENTIFIER integerSubType
 | OCTET STRING
 | OCTET STRING octetStringSubType
 | UPPERCASE _ IDENTIFIER octetStringSubType
 | OBJECT IDENTIFIER anySubType"""


signature:
def p_valueofSimpleSyntax ( self , p ) :
    """valueofSimpleSyntax : NUMBER
 | NEGATIVENUMBER
 | NUMBER64
 | NEGATIVENUMBER64
 | HEX _ STRING
 | BIN _ STRING
 | LOWERCASE _ IDENTIFIER
 | QUOTED _ STRING
 | ' { ' objectIdentifier _ defval ' } '"""

implementation:
# libsmi for objectIdentifier _ defval :
# This is only for some MIBs with invalid numerical
# OID notation for DEFVALs . We DO NOT parse them
# correctly . We just don ' t want to produce a
# parser error .
n = len ( p )
if n == 2 :
    p [ 0 ] = p [ 1 ]
elif n == 4 : # XXX
    pass


signature:
def p_sequenceSimpleSyntax ( self , p ) :
    """sequenceSimpleSyntax : INTEGER anySubType
 | INTEGER32 anySubType
 | OCTET STRING anySubType
 | OBJECT IDENTIFIER anySubType"""

implementation:
n = len ( p )
if n == 3 :
    p [ 0 ] = p [ 1 ]
    # XXX not supporting subtypes here
elif n == 4 :
    p [ 0 ] = p [ 1 ] + ' ' + p [ 2 ]


signature:
def p_ApplicationSyntax ( self , p ) :
    """ApplicationSyntax : IPADDRESS anySubType
 | COUNTER32
 | COUNTER32 integerSubType
 | GAUGE32
 | GAUGE32 integerSubType
 | UNSIGNED32
 | UNSIGNED32 integerSubType
 | TIMETICKS anySubType
 | OPAQUE
 | OPAQUE octetStringSubType
 | COUNTER64
 | COUNTER64 integerSubType"""

implementation:
# COUNTER32 and COUNTER64 was with anySubType in libsmi
n = len ( p )
if n == 2 :
    p [ 0 ] = ( 'ApplicationSyntax' , p [ 1 ] )
elif n == 3 :
    p [ 0 ] = ( 'ApplicationSyntax' , p [ 1 ] , p [ 2 ] )


signature:
def p_sequenceApplicationSyntax ( self , p ) :
    """sequenceApplicationSyntax : IPADDRESS anySubType
 | COUNTER32 anySubType
 | GAUGE32 anySubType
 | UNSIGNED32 anySubType
 | TIMETICKS anySubType
 | OPAQUE
 | COUNTER64 anySubType"""


signature:
def p_ranges ( self , p ) :
    """ranges : ranges ' | ' range
 | range"""


signature:
def p_range ( self , p ) :
    """range : value DOT _ DOT value
 | value"""


signature:
def p_enumItems ( self , p ) :
    """enumItems : enumItems ' , ' enumItem
 | enumItem"""


signature:
def p_IndexTypes ( self , p ) :
    """IndexTypes : IndexTypes ' , ' IndexType
 | IndexType"""


signature:
def p_IndexType ( self , p ) :
    """IndexType : IMPLIED Index
 | Index"""


signature:
def p_Value ( self , p ) :
    """Value : valueofObjectSyntax
 | ' { ' BitsValue ' } '"""


signature:
def p_BitNames ( self , p ) :
    """BitNames : BitNames ' , ' LOWERCASE _ IDENTIFIER
 | LOWERCASE _ IDENTIFIER"""


signature:
def p_Revisions ( self , p ) :
    """Revisions : Revisions Revision
 | Revision"""


signature:
def p_Objects ( self , p ) :
    """Objects : Objects ' , ' Object
 | Object"""


signature:
def p_Notifications ( self , p ) :
    """Notifications : Notifications ' , ' Notification
 | Notification"""


signature:
def p_subidentifiers ( self , p ) :
    """subidentifiers : subidentifiers subidentifier
 | subidentifier"""


signature:
def p_subidentifier ( self , p ) :
    """subidentifier : fuzzy _ lowercase _ identifier
 | NUMBER
 | LOWERCASE _ IDENTIFIER ' ( ' NUMBER ' ) '"""

implementation:
n = len ( p )
if n == 2 :
    p [ 0 ] = p [ 1 ]
elif n == 5 : # NOTE : we are not creating new symbol p [ 1 ] because formally
# it is not defined in * this * MIB
    p [ 0 ] = ( p [ 1 ] , p [ 3 ] )


signature:
def p_subidentifiers_defval ( self , p ) :
    """subidentifiers _ defval : subidentifiers _ defval subidentifier _ defval
 | subidentifier _ defval"""

implementation:
n = len ( p )
if n == 3 :
    p [ 0 ] = ( 'subidentifiers_defval' , p [ 1 ] [ 1 ] + [ p [ 2 ] ] )
elif n == 2 :
    p [ 0 ] = ( 'subidentifiers_defval' , [ p [ 1 ] ] )


signature:
def p_subidentifier_defval ( self , p ) :
    """subidentifier _ defval : LOWERCASE _ IDENTIFIER ' ( ' NUMBER ' ) '
 | NUMBER"""

implementation:
n = len ( p )
if n == 2 :
    p [ 0 ] = ( 'subidentifier_defval' , p [ 1 ] )
elif n == 5 :
    p [ 0 ] = ( 'subidentifier_defval' , p [ 1 ] , p [ 3 ] )


signature:
def p_objectGroupClause ( self , p ) :
    """objectGroupClause : LOWERCASE _ IDENTIFIER OBJECT _ GROUP ObjectGroupObjectsPart STATUS Status DESCRIPTION Text ReferPart COLON _ COLON _ EQUAL ' { ' objectIdentifier ' } '"""

implementation:
p [ 0 ] = ( 'objectGroupClause' , p [ 1 ] , # id
p [ 3 ] , # objects
p [ 5 ] , # status
( p [ 6 ] , p [ 7 ] ) , # description
p [ 8 ] , # reference
p [ 11 ] )


signature:
def p_notificationGroupClause ( self , p ) :
    """notificationGroupClause : LOWERCASE _ IDENTIFIER NOTIFICATION _ GROUP NotificationsPart STATUS Status DESCRIPTION Text ReferPart COLON _ COLON _ EQUAL ' { ' objectIdentifier ' } '"""

implementation:
p [ 0 ] = ( 'notificationGroupClause' , p [ 1 ] , # id
p [ 3 ] , # notifications
p [ 5 ] , # status
( p [ 6 ] , p [ 7 ] ) , # description
p [ 8 ] , # reference
p [ 11 ] )


signature:
def p_moduleComplianceClause ( self , p ) :
    """moduleComplianceClause : LOWERCASE _ IDENTIFIER MODULE _ COMPLIANCE STATUS Status DESCRIPTION Text ReferPart ComplianceModulePart COLON _ COLON _ EQUAL ' { ' objectIdentifier ' } '"""

implementation:
p [ 0 ] = ( 'moduleComplianceClause' , p [ 1 ] , # id
# p [ 2 ] , # MODULE _ COMPLIANCE
p [ 4 ] , # status
( p [ 5 ] , p [ 6 ] ) , # description
p [ 7 ] , # reference
p [ 8 ] , # ComplianceModules
p [ 11 ] )


signature:
def p_ComplianceModules ( self , p ) :
    """ComplianceModules : ComplianceModules ComplianceModule
 | ComplianceModule"""

implementation:
n = len ( p )
if n == 3 :
    p [ 0 ] = ( 'ComplianceModules' , p [ 1 ] [ 1 ] + [ p [ 2 ] ] )
elif n == 2 :
    p [ 0 ] = ( 'ComplianceModules' , [ p [ 1 ] ] )


signature:
def p_ComplianceModule ( self , p ) :
    """ComplianceModule : MODULE ComplianceModuleName MandatoryPart CompliancePart"""

implementation:
objects = p [ 3 ] and p [ 3 ] [ 1 ] or [ ]
objects += p [ 4 ] and p [ 4 ] [ 1 ] or [ ]
p [ 0 ] = ( p [ 2 ] , # ModuleName
objects )


signature:
def p_MandatoryGroups ( self , p ) :
    """MandatoryGroups : MandatoryGroups ' , ' MandatoryGroup
 | MandatoryGroup"""

implementation:
n = len ( p )
if n == 4 :
    p [ 0 ] = ( 'MandatoryGroups' , p [ 1 ] [ 1 ] + [ p [ 3 ] ] )
elif n == 2 :
    p [ 0 ] = ( 'MandatoryGroups' , [ p [ 1 ] ] )


signature:
def p_Compliances ( self , p ) :
    """Compliances : Compliances Compliance
 | Compliance"""

implementation:
n = len ( p )
if n == 3 :
    p [ 0 ] = p [ 1 ] and p [ 2 ] and ( 'Compliances' , p [ 1 ] [ 1 ] + [ p [ 2 ] ] ) or p [ 1 ]
elif n == 2 :
    p [ 0 ] = p [ 1 ] and ( 'Compliances' , [ p [ 1 ] ] ) or None


signature:
def p_agentCapabilitiesClause ( self , p ) :
    """agentCapabilitiesClause : LOWERCASE _ IDENTIFIER AGENT _ CAPABILITIES PRODUCT _ RELEASE Text STATUS Status DESCRIPTION Text ReferPart ModulePart _ Capabilities COLON _ COLON _ EQUAL ' { ' objectIdentifier ' } '"""

implementation:
p [ 0 ] = ( 'agentCapabilitiesClause' , p [ 1 ] , # id
# p [ 2 ] , # AGENT _ CAPABILITIES
( p [ 3 ] , p [ 4 ] ) , # product release
p [ 6 ] , # status
( p [ 7 ] , p [ 8 ] ) , # description
p [ 9 ] , # reference
# p [ 10 ] , # module capabilities
p [ 13 ] )


signature:
def p_Cells ( self , p ) :
    """Cells : Cells ' , ' Cell
 | Cell"""


signature:
def p_typeSMIv1 ( self , p ) :
    """typeSMIv1 : INTEGER
 | OCTET STRING
 | IPADDRESS
 | NETWORKADDRESS"""

implementation:
n = len ( p )
indextype = n == 3 and p [ 1 ] + ' ' + p [ 2 ] or p [ 1 ]
p [ 0 ] = indextype


signature:
def p_enumItems ( self , p ) :
    """enumItems : enumItems ' , ' enumItem
 | enumItem
 | enumItems enumItem
 | enumItems ' , '"""


signature:
def p_notificationTypeClause ( self , p ) :
    """notificationTypeClause : fuzzy _ lowercase _ identifier NOTIFICATION _ TYPE NotificationObjectsPart STATUS Status DESCRIPTION Text ReferPart COLON _ COLON _ EQUAL ' { ' NotificationName ' } '"""

implementation:
# some MIBs have uppercase and / or lowercase id
p [ 0 ] = ( 'notificationTypeClause' , p [ 1 ] , # id
# p [ 2 ] , # NOTIFICATION _ TYPE
p [ 3 ] , # NotificationObjectsPart
p [ 5 ] , # status
( p [ 6 ] , p [ 7 ] ) , # description
p [ 8 ] , # Reference
p [ 11 ] )


signature:
def p_trapTypeClause ( self , p ) :
    """trapTypeClause : fuzzy _ lowercase _ identifier TRAP _ TYPE EnterprisePart VarPart DescrPart ReferPart COLON _ COLON _ EQUAL NUMBER"""

implementation:
# libsmi : TODO : range of number ?
p [ 0 ] = ( 'trapTypeClause' , p [ 1 ] , # fuzzy _ lowercase _ identifier
# p [ 2 ] , # TRAP _ TYPE
p [ 3 ] , # EnterprisePart ( objectIdentifier )
p [ 4 ] , # VarPart
p [ 5 ] , # description
p [ 6 ] , # reference
p [ 8 ] )


signature:
def p_EnterprisePart ( self , p ) :
    """EnterprisePart : ENTERPRISE objectIdentifier
 | ENTERPRISE ' { ' objectIdentifier ' } '"""

implementation:
n = len ( p )
if n == 3 :
    p [ 0 ] = p [ 2 ]
elif n == 5 : # common mistake case
    p [ 0 ] = p [ 3 ]


signature:
def p_CreationPart ( self , p ) :
    """CreationPart : CREATION _ REQUIRES ' { ' Cells ' } '
 | CREATION _ REQUIRES ' { ' ' } '
 | empty"""


signature:
def AIC ( N , rho , k ) :
    r"""Akaike Information Criterion
 : param rho : rho at order k
 : param N : sample size
 : param k : AR order .
 If k is the AR order and N the size of the sample , then Akaike criterion is
 . . math : : AIC ( k ) = \ log ( \ rho _ k ) + 2 \ frac { k + 1 } { N }
 AIC ( 64 , [ 0.5,0.3,0.2 ] , [ 1,2,3 ] )
 : validation : double checked versus octave ."""

implementation:
from numpy import log , array
# k + 1 # todo check convention . agrees with octave
res = N * log ( array ( rho ) ) + 2. * ( array ( k ) + 1 )
return res


signature:
def AICc ( N , rho , k , norm = True ) :
    r"""corrected Akaike information criterion
 . . math : : AICc ( k ) = log ( \ rho _ k ) + 2 \ frac { k + 1 } { N - k - 2}
 : validation : double checked versus octave ."""

implementation:
from numpy import log , array
p = k
# todo check convention . agrees with octave
res = log ( rho ) + 2. * ( p + 1 ) / ( N - p - 2 )
return res


signature:
def KIC ( N , rho , k ) :
    r"""Kullback information criterion
 . . math : : KIC ( k ) = log ( \ rho _ k ) + 3 \ frac { k + 1 } { N }
 : validation : double checked versus octave ."""

implementation:
from numpy import log , array
res = log ( rho ) + 3. * ( k + 1. ) / float ( N )
return res


signature:
def AKICc ( N , rho , k ) :
    r"""approximate corrected Kullback information
 . . math : : AKICc ( k ) = log ( rho _ k ) + \ frac { p } { N * ( N - k ) } + ( 3 - \ frac { k + 2 } { N } ) * \ frac { k + 1 } { N - k - 2}"""

implementation:
from numpy import log , array
p = k
res = log ( rho ) + p / N / ( N - p ) + ( 3. - ( p + 2. ) / N ) * ( p + 1. ) / ( N - p - 2. )
return res


signature:
def FPE ( N , rho , k = None ) :
    r"""Final prediction error criterion
 . . math : : FPE ( k ) = \ frac { N + k + 1 } { N - k - 1 } \ rho _ k
 : validation : double checked versus octave ."""

implementation:
# k # todo check convention . agrees with octave
fpe = rho * ( N + k + 1. ) / ( N - k - 1 )
return fpe


signature:
def MDL ( N , rho , k ) :
    r"""Minimum Description Length
 . . math : : MDL ( k ) = N log \ rho _ k + p \ log N
 : validation : results"""

implementation:
from numpy import log
# p = arange ( 1 , len ( rho ) + 1)
mdl = N * log ( rho ) + k * log ( N )
return mdl


signature:
def CAT ( N , rho , k ) :
    r"""Criterion Autoregressive Transfer Function :
 . . math : : CAT ( k ) = \ frac { 1 } { N } \ sum _ { i = 1 } ^ k \ frac { 1 } { \ rho _ i } - \ frac { \ rho _ i } { \ rho _ k }
 . . todo : : validation"""

implementation:
from numpy import zeros , arange
cat = zeros ( len ( rho ) )
for p in arange ( 1 , len ( rho ) + 1 ) :
    rho_p = float ( N ) / ( N - p ) * rho [ p - 1 ]
    s = 0
    for j in range ( 1 , p + 1 ) :
        rho_j = float ( N ) / ( N - j ) * rho [ j - 1 ]
        s = s + 1. / rho_j
    # print ( s , s / float ( N ) , 1 . / rho _ p )
    cat [ p - 1 ] = s / float ( N ) - 1. / rho_p
return cat


signature:
def aic_eigen ( s , N ) :
    r"""AIC order - selection using eigen values
 : param s : a list of ` p ` sorted eigen values
 : param N : the size of the input data . To be defined precisely .
 : return :
 * an array containing the AIC values
 Given : math : ` n ` sorted eigen values : math : ` \ lambda _ i ` with
 : math : ` 0 < = i < n ` , the proposed criterion from Wax and Kailath ( 1985)
 is :
 . . math : : AIC ( k ) = - 2 ( n - k ) N \ ln \ frac { g ( k ) } { a ( k ) } + 2k ( 2n - k )
 where the arithmetic sum : math : ` a ( k ) ` is :
 . . math : : a ( k ) = \ sum _ { i = k + 1 } ^ { n } \ lambda _ i
 and the geometric sum : math : ` g ( k ) ` is :
 . . math : : g ( k ) = \ prod _ { i = k + 1 } ^ { n } \ lambda _ i ^ { - ( n - k ) }
 The number of relevant sinusoids in the signal subspace is determined by
 selecting the minimum of ` AIC ` .
 . . seealso : : : func : ` ~ spectrum . eigenfreq . eigen `
 . . todo : : define precisely the input parameter N . Should be the input
 data length but when using correlation matrix ( SVD ) , I suspect it
 should be the length of the correlation matrix rather than the
 original data .
 : References :
 * [ Marple ] _ Chap 13,
 * [ Wax ] _"""

implementation:
import numpy as np
kaic = [ ]
n = len ( s )
for k in range ( 0 , n - 1 ) :
    ak = 1. / ( n - k ) * np . sum ( s [ k + 1 : ] )
    gk = np . prod ( s [ k + 1 : ] ** ( 1. / ( n - k ) ) )
    kaic . append ( - 2. * ( n - k ) * N * np . log ( gk / ak ) + 2. * k * ( 2. * n - k ) )
return kaic


signature:
def mdl_eigen ( s , N ) :
    r"""MDL order - selection using eigen values
 : param s : a list of ` p ` sorted eigen values
 : param N : the size of the input data . To be defined precisely .
 : return :
 * an array containing the AIC values
 . . math : : MDL ( k ) = ( n - k ) N \ ln \ frac { g ( k ) } { a ( k ) } + 0.5k ( 2n - k ) log ( N )
 . . seealso : : : func : ` aic _ eigen ` for details
 : References :
 * [ Marple ] _ Chap 13,
 * [ Wax ] _"""

implementation:
import numpy as np
kmdl = [ ]
n = len ( s )
for k in range ( 0 , n - 1 ) :
    ak = 1. / ( n - k ) * np . sum ( s [ k + 1 : ] )
    gk = np . prod ( s [ k + 1 : ] ** ( 1. / ( n - k ) ) )
    kmdl . append ( - ( n - k ) * N * np . log ( gk / ak ) + 0.5 * k * ( 2. * n - k ) * np . log ( N ) )
return kmdl


signature:
def generate_gallery_rst ( app ) :
    """Generate the Main examples gallery reStructuredText
 Start the sphinx - gallery configuration and recursively scan the examples
 directories in order to populate the examples gallery"""

implementation:
try :
    plot_gallery = eval ( app . builder . config . plot_gallery )
except TypeError :
    plot_gallery = bool ( app . builder . config . plot_gallery )
gallery_conf . update ( app . config . sphinx_gallery_conf )
gallery_conf . update ( plot_gallery = plot_gallery )
gallery_conf . update ( abort_on_example_error = app . builder . config . abort_on_example_error )
# this assures I can call the config in other places
app . config . sphinx_gallery_conf = gallery_conf
app . config . html_static_path . append ( glr_path_static ( ) )
clean_gallery_out ( app . builder . outdir )
examples_dirs = gallery_conf [ 'examples_dirs' ]
gallery_dirs = gallery_conf [ 'gallery_dirs' ]
if not isinstance ( examples_dirs , list ) :
    examples_dirs = [ examples_dirs ]
if not isinstance ( gallery_dirs , list ) :
    gallery_dirs = [ gallery_dirs ]
mod_examples_dir = os . path . relpath ( gallery_conf [ 'mod_example_dir' ] , app . builder . srcdir )
seen_backrefs = set ( )
for examples_dir , gallery_dir in zip ( examples_dirs , gallery_dirs ) :
    examples_dir = os . path . relpath ( examples_dir , app . builder . srcdir )
    gallery_dir = os . path . relpath ( gallery_dir , app . builder . srcdir )
    for workdir in [ examples_dir , gallery_dir , mod_examples_dir ] :
        if not os . path . exists ( workdir ) :
            os . makedirs ( workdir )
    # we create an index . rst with all examples
    fhindex = open ( os . path . join ( gallery_dir , 'index.rst' ) , 'w' )
    # Here we don ' t use an os . walk , but we recurse only twice : flat is
    # better than nested .
    fhindex . write ( generate_dir_rst ( examples_dir , gallery_dir , gallery_conf , seen_backrefs ) )
    for directory in sorted ( os . listdir ( examples_dir ) ) :
        if os . path . isdir ( os . path . join ( examples_dir , directory ) ) :
            src_dir = os . path . join ( examples_dir , directory )
            target_dir = os . path . join ( gallery_dir , directory )
            fhindex . write ( generate_dir_rst ( src_dir , target_dir , gallery_conf , seen_backrefs ) )
    fhindex . flush ( )


signature:
def setup ( app ) :
    """Setup sphinx - gallery sphinx extension"""

implementation:
app . add_config_value ( 'plot_gallery' , True , 'html' )
app . add_config_value ( 'abort_on_example_error' , False , 'html' )
app . add_config_value ( 'sphinx_gallery_conf' , gallery_conf , 'html' )
app . add_stylesheet ( 'gallery.css' )
app . connect ( 'builder-inited' , generate_gallery_rst )
app . connect ( 'build-finished' , embed_code_links )


signature:
def CORRELATION ( x , y = None , maxlags = None , norm = 'unbiased' ) :
    r"""Correlation function
 This function should give the same results as : func : ` xcorr ` but it
 returns the positive lags only . Moreover the algorithm does not use
 FFT as compared to other algorithms .
 : param array x : first data array of length N
 : param array y : second data array of length N . If not specified , computes the
 autocorrelation .
 : param int maxlags : compute cross correlation between [ 0 : maxlags ]
 when maxlags is not specified , the range of lags is [ 0 : maxlags ] .
 : param str norm : normalisation in [ ' biased ' , ' unbiased ' , None , ' coeff ' ]
 * * biased * correlation = raw / N ,
 * * unbiased * correlation = raw / ( N - ` | lag | ` )
 * * coeff * correlation = raw / ( rms ( x ) . rms ( y ) ) / N
 * None correlation = raw
 : return :
 * a numpy . array correlation sequence , r [ 1 , N ]
 * a float for the zero - lag correlation , r [ 0]
 The * unbiased * correlation has the form :
 . . math : :
 \ hat { r } _ { xx } = \ frac { 1 } { N - m } T \ sum _ { n = 0 } ^ { N - m - 1 } x [ n + m ] x ^ * [ n ] T
 The * biased * correlation differs by the front factor only :
 . . math : :
 \ check { r } _ { xx } = \ frac { 1 } { N } T \ sum _ { n = 0 } ^ { N - m - 1 } x [ n + m ] x ^ * [ n ] T
 with : math : ` 0 \ leq m \ leq N - 1 ` .
 . . doctest : :
 > > > from spectrum import CORRELATION
 > > > x = [ 1,2,3,4,5]
 > > > res = CORRELATION ( x , x , maxlags = 0 , norm = ' biased ' )
 > > > res [ 0]
 11.0
 . . note : : this function should be replaced by : func : ` xcorr ` .
 . . seealso : : : func : ` xcorr `"""

implementation:
assert norm in [ 'unbiased' , 'biased' , 'coeff' , None ]
# transform lag into list if it is an integer
x = np . array ( x )
if y is None :
    y = x
else :
    y = np . array ( y )
# N is the max of x and y
N = max ( len ( x ) , len ( y ) )
if len ( x ) < N :
    x = y . copy ( )
    x . resize ( N )
if len ( y ) < N :
    y = y . copy ( )
    y . resize ( N )
# default lag is N - 1
if maxlags is None :
    maxlags = N - 1
assert maxlags < N , 'lag must be less than len(x)'
realdata = np . isrealobj ( x ) and np . isrealobj ( y )
# create an autocorrelation array with same length as lag
if realdata == True :
    r = np . zeros ( maxlags , dtype = float )
else :
    r = np . zeros ( maxlags , dtype = complex )
if norm == 'coeff' :
    rmsx = pylab_rms_flat ( x )
    rmsy = pylab_rms_flat ( y )
for k in range ( 0 , maxlags + 1 ) :
    nk = N - k - 1
    if realdata == True :
        sum = 0
        for j in range ( 0 , nk + 1 ) :
            sum = sum + x [ j + k ] * y [ j ]
    else :
        sum = 0. + 0j
        for j in range ( 0 , nk + 1 ) :
            sum = sum + x [ j + k ] * y [ j ] . conjugate ( )
    if k == 0 :
        if norm in [ 'biased' , 'unbiased' ] :
            r0 = sum / float ( N )
        elif norm is None :
            r0 = sum
        else :
            r0 = 1.
    else :
        if norm == 'unbiased' :
            r [ k - 1 ] = sum / float ( N - k )
        elif norm == 'biased' :
            r [ k - 1 ] = sum / float ( N )
        elif norm is None :
            r [ k - 1 ] = sum
        elif norm == 'coeff' :
            r [ k - 1 ] = sum / ( rmsx * rmsy ) / float ( N )
r = np . insert ( r , 0 , r0 )
return r


signature:
def xcorr ( x , y = None , maxlags = None , norm = 'biased' ) :
    """Cross - correlation using numpy . correlate
 Estimates the cross - correlation ( and autocorrelation ) sequence of a random
 process of length N . By default , there is no normalisation and the output
 sequence of the cross - correlation has a length 2 * N + 1.
 : param array x : first data array of length N
 : param array y : second data array of length N . If not specified , computes the
 autocorrelation .
 : param int maxlags : compute cross correlation between [ - maxlags : maxlags ]
 when maxlags is not specified , the range of lags is [ - N + 1 : N - 1 ] .
 : param str option : normalisation in [ ' biased ' , ' unbiased ' , None , ' coeff ' ]
 The true cross - correlation sequence is
 . . math : : r _ { xy } [ m ] = E ( x [ n + m ] . y ^ * [ n ] ) = E ( x [ n ] . y ^ * [ n - m ] )
 However , in practice , only a finite segment of one realization of the
 infinite - length random process is available .
 The correlation is estimated using numpy . correlate ( x , y , ' full ' ) .
 Normalisation is handled by this function using the following cases :
 * ' biased ' : Biased estimate of the cross - correlation function
 * ' unbiased ' : Unbiased estimate of the cross - correlation function
 * ' coeff ' : Normalizes the sequence so the autocorrelations at zero
 lag is 1.0.
 : return :
 * a numpy . array containing the cross - correlation sequence ( length 2 * N - 1)
 * lags vector
 . . note : : If x and y are not the same length , the shorter vector is
 zero - padded to the length of the longer vector .
 . . rubric : : Examples
 . . doctest : :
 > > > from spectrum import xcorr
 > > > x = [ 1,2,3,4,5]
 > > > c , l = xcorr ( x , x , maxlags = 0 , norm = ' biased ' )
 array ( [ 11 . ] )
 . . seealso : : : func : ` CORRELATION ` ."""

implementation:
N = len ( x )
if y is None :
    y = x
assert len ( x ) == len ( y ) , 'x and y must have the same length. Add zeros if needed'
if maxlags is None :
    maxlags = N - 1
    lags = np . arange ( 0 , 2 * N - 1 )
else :
    assert maxlags <= N , 'maxlags must be less than data length'
    lags = np . arange ( N - maxlags - 1 , N + maxlags )
res = np . correlate ( x , y , mode = 'full' )
if norm == 'biased' :
    Nf = float ( N )
    res = res [ lags ] / float ( N )
    # do not use / = ! !
elif norm == 'unbiased' :
    res = res [ lags ] / ( float ( N ) - abs ( np . arange ( - N + 1 , N ) ) ) [ lags ]
elif norm == 'coeff' :
    Nf = float ( N )
    rms = pylab_rms_flat ( x ) * pylab_rms_flat ( y )
    res = res [ lags ] / rms / Nf
else :
    res = res [ lags ]
lags = np . arange ( - maxlags , maxlags + 1 )
return res , lags


signature:
def MINEIGVAL ( T0 , T , TOL ) :
    """Finds the minimum eigenvalue of a Hermitian Toeplitz matrix
 The classical power method is used together with a fast Toeplitz
 equation solution routine . The eigenvector is normalized to unit length .
 : param T0 : Scalar corresponding to real matrix element t ( 0)
 : param T : Array of M complex matrix elements t ( 1 ) , . . . , t ( M ) C from the left column of the Toeplitz matrix
 : param TOL : Real scalar tolerance ; routine exits when [ EVAL ( k ) - EVAL ( k - 1 ) ] / EVAL ( k - 1 ) < TOL , where the index k denotes the iteration number .
 : return :
 * EVAL - Real scalar denoting the minimum eigenvalue of matrix
 * EVEC - Array of M complex eigenvector elements associated
 . . note : :
 * External array T must be dimensioned > = M
 * array EVEC must be > = M + 1
 * Internal array E must be dimensioned > = M + 1 .
 * * * dependencies * *
 * : meth : ` spectrum . toeplitz . HERMTOEP `"""

implementation:
M = len ( T )
eigval = 10
eigvalold = 1
eigvec = numpy . zeros ( M + 1 , dtype = complex )
for k in range ( 0 , M + 1 ) :
    eigvec [ k ] = 1 + 0j
it = 0
# print ' initialisation ' , T0 , T , eigval , eigvec
maxit = 15
while abs ( eigvalold - eigval ) > TOL * eigvalold and it < maxit :
    it = it + 1
    eigvalold = eigval
    # print ' iteration ' , it , ' eigvalold = ' , eigvalold , ' eigval = ' , eigval
    eig = toeplitz . HERMTOEP ( T0 , T , eigvec )
    SUM = 0
    save = 0. + 0j
    for k in range ( 0 , M + 1 ) :
        SUM = SUM + eig [ k ] . real ** 2 + eig [ k ] . imag ** 2
        save = save + eig [ k ] * eigvec [ k ] . conjugate ( )
    SUM = 1. / SUM
    eigval = save . real * SUM
    for k in range ( 0 , M + 1 ) :
        eigvec [ k ] = SUM * eig [ k ]
if it == maxit :
    print ( 'warning reached max number of iteration (%s)' % maxit )
return eigval , eigvec


signature:
def morlet ( lb , ub , n ) :
    r"""Generate the Morlet waveform
 The Morlet waveform is defined as follows :
 . . math : : w [ x ] = \ cos { 5x } \ exp ^ { - x ^ 2/2}
 : param lb : lower bound
 : param ub : upper bound
 : param int n : waveform data samples
 . . plot : :
 : include - source :
 : width : 80%
 from spectrum import morlet
 from pylab import plot
 plot ( morlet ( 0,10,100 ) )"""

implementation:
if n <= 0 :
    raise ValueError ( "n must be strictly positive" )
x = numpy . linspace ( lb , ub , n )
psi = numpy . cos ( 5 * x ) * numpy . exp ( - x ** 2 / 2. )
return psi


signature:
def chirp ( t , f0 = 0. , t1 = 1. , f1 = 100. , form = 'linear' , phase = 0 ) :
    r"""Evaluate a chirp signal at time t .
 A chirp signal is a frequency swept cosine wave .
 . . math : : a = \ pi ( f _ 1 - f _ 0 ) / t _ 1
 . . math : : b = 2 \ pi f _ 0
 . . math : : y = \ cos \ left ( \ pi \ frac { f _ 1 - f _ 0 } { t _ 1 } t ^ 2 + 2 \ pi f _ 0 t + \ rm { phase } \ right )
 : param array t : times at which to evaluate the chirp signal
 : param float f0 : frequency at time t = 0 ( Hz )
 : param float t1 : time t1
 : param float f1 : frequency at time t = t1 ( Hz )
 : param str form : shape of frequency sweep in [ ' linear ' , ' quadratic ' , ' logarithmic ' ]
 : param float phase : phase shift at t = 0
 The parameter * * form * * can be :
 * ' linear ' : math : ` f ( t ) = ( f _ 1 - f _ 0 ) ( t / t _ 1 ) + f _ 0 `
 * ' quadratic ' : math : ` f ( t ) = ( f _ 1 - f _ 0 ) ( t / t _ 1 ) ^ 2 + f _ 0 `
 * ' logarithmic ' : math : ` f ( t ) = ( f _ 1 - f _ 0 ) ^ { ( t / t _ 1 ) } + f _ 0 `
 Example :
 . . plot : :
 : include - source :
 : width : 80%
 from spectrum import chirp
 from pylab import linspace , plot
 t = linspace ( 0 , 1 , 1000)
 y = chirp ( t , form = ' linear ' )
 plot ( y )
 y = chirp ( t , form = ' quadratic ' )
 plot ( y , ' r ' )"""

implementation:
valid_forms = [ 'linear' , 'quadratic' , 'logarithmic' ]
if form not in valid_forms :
    raise ValueError ( "Invalid form. Valid form are %s" % valid_forms )
t = numpy . array ( t )
phase = 2. * pi * phase / 360.
if form == "linear" :
    a = pi * ( f1 - f0 ) / t1
    b = 2. * pi * f0
    y = numpy . cos ( a * t ** 2 + b * t + phase )
elif form == "quadratic" :
    a = ( 2 / 3. * pi * ( f1 - f0 ) / t1 / t1 )
    b = 2. * pi * f0
    y = numpy . cos ( a * t ** 3 + b * t + phase )
elif form == "logarithmic" :
    a = 2. * pi * t1 / numpy . log ( f1 - f0 )
    b = 2. * pi * f0
    x = ( f1 - f0 ) ** ( 1. / t1 )
    y = numpy . cos ( a * x ** t + b * t + phase )
return y


signature:
def mexican ( lb , ub , n ) :
    r"""Generate the mexican hat wavelet
 The Mexican wavelet is :
 . . math : : w [ x ] = \ cos { 5x } \ exp ^ { - x ^ 2/2}
 : param lb : lower bound
 : param ub : upper bound
 : param int n : waveform data samples
 : return : the waveform
 . . plot : :
 : include - source :
 : width : 80%
 from spectrum import mexican
 from pylab import plot
 plot ( mexican ( 0 , 10 , 100 ) )"""

implementation:
if n <= 0 :
    raise ValueError ( "n must be strictly positive" )
x = numpy . linspace ( lb , ub , n )
psi = ( 1. - x ** 2. ) * ( 2. / ( numpy . sqrt ( 3. ) * pi ** 0.25 ) ) * numpy . exp ( - x ** 2 / 2. )
return psi


signature:
def ac2poly ( data ) :
    """Convert autocorrelation sequence to prediction polynomial
 : param array data : input data ( list or numpy . array )
 : return :
 * AR parameters
 * noise variance
 This is an alias to : :
 a , e , c = LEVINSON ( data )
 : Example :
 . . doctest : :
 > > > from spectrum import ac2poly
 > > > from numpy import array
 > > > r = [ 5 , - 2 , 1.01]
 > > > ar , e = ac2poly ( r )
 > > > ar
 array ( [ 1 . , 0.38 , - 0.05 ] )
 4.1895000007"""

implementation:
a , e , _c = LEVINSON ( data )
a = numpy . insert ( a , 0 , 1 )
return a , e


signature:
def rc2poly ( kr , r0 = None ) :
    """convert reflection coefficients to prediction filter polynomial
 : param k : reflection coefficients"""

implementation:
# Initialize the recursion
from . levinson import levup
p = len ( kr )
# % p is the order of the prediction polynomial .
a = numpy . array ( [ 1 , kr [ 0 ] ] )
# % a is a true polynomial .
e = numpy . zeros ( len ( kr ) )
if r0 is None :
    e0 = 0
else :
    e0 = r0
e [ 0 ] = e0 * ( 1. - numpy . conj ( numpy . conjugate ( kr [ 0 ] ) * kr [ 0 ] ) )
# Continue the recursion for k = 2,3 , . . . , p , where p is the order of the
# prediction polynomial .
for k in range ( 1 , p ) :
    [ a , e [ k ] ] = levup ( a , kr [ k ] , e [ k - 1 ] )
efinal = e [ - 1 ]
return a , efinal


signature:
def rc2ac ( k , R0 ) :
    """Convert reflection coefficients to autocorrelation sequence .
 : param k : reflection coefficients
 : param R0 : zero - lag autocorrelation
 : returns : the autocorrelation sequence
 . . seealso : : : func : ` ac2rc ` , : func : ` poly2rc ` , : func : ` ac2poly ` , : func : ` poly2rc ` , : func : ` rc2poly ` ."""

implementation:
[ a , efinal ] = rc2poly ( k , R0 )
R , u , kr , e = rlevinson ( a , efinal )
return R


signature:
def rc2is ( k ) :
    """Convert reflection coefficients to inverse sine parameters .
 : param k : reflection coefficients
 : return : inverse sine parameters
 . . seealso : : : func : ` is2rc ` , : func : ` rc2poly ` , : func : ` rc2acC ` , : func : ` rc2lar ` .
 Reference : J . R . Deller , J . G . Proakis , J . H . L . Hansen , " Discrete - Time
 Processing of Speech Signals " , Prentice Hall , Section 7.4.5."""

implementation:
assert numpy . isrealobj ( k ) , 'Inverse sine parameters not defined for complex reflection coefficients.'
if max ( numpy . abs ( k ) ) >= 1 :
    raise ValueError ( 'All reflection coefficients should have magnitude less than unity.' )
return ( 2 / numpy . pi ) * numpy . arcsin ( k )


signature:
def rc2lar ( k ) :
    """Convert reflection coefficients to log area ratios .
 : param k : reflection coefficients
 : return : inverse sine parameters
 The log area ratio is defined by G = log ( ( 1 + k ) / ( 1 - k ) ) , where the K
 parameter is the reflection coefficient .
 . . seealso : : : func : ` lar2rc ` , : func : ` rc2poly ` , : func : ` rc2ac ` , : func : ` rc2ic ` .
 : References :
 [1 ] J . Makhoul , " Linear Prediction : A Tutorial Review , " Proc . IEEE , Vol . 63 , No . 4 , pp . 561-580 , Apr 1975."""

implementation:
assert numpy . isrealobj ( k ) , 'Log area ratios not defined for complex reflection coefficients.'
if max ( numpy . abs ( k ) ) >= 1 :
    raise ValueError ( 'All reflection coefficients should have magnitude less than unity.' )
# Use the relation , atanh ( x ) = ( 1/2 ) * log ( ( 1 + k ) / ( 1 - k ) )
return - 2 * numpy . arctanh ( - numpy . array ( k ) )


signature:
def lar2rc ( g ) :
    """Convert log area ratios to reflection coefficients .
 : param g : log area ratios
 : returns : the reflection coefficients
 . . seealso : : func : ` rc2lar ` , : func : ` poly2rc ` , : func : ` ac2rc ` , : func : ` is2rc ` .
 : References :
 [1 ] J . Makhoul , " Linear Prediction : A Tutorial Review , " Proc . IEEE , Vol . 63 , No . 4 , pp . 561-580 , Apr 1975."""

implementation:
assert numpy . isrealobj ( g ) , 'Log area ratios not defined for complex reflection coefficients.'
# Use the relation , tanh ( x ) = ( 1 - exp ( 2x ) ) / ( 1 + exp ( 2x ) )
return - numpy . tanh ( - numpy . array ( g ) / 2 )


signature:
def lsf2poly ( lsf ) :
    """Convert line spectral frequencies to prediction filter coefficients
 returns a vector a containing the prediction filter coefficients from a vector lsf of line spectral frequencies .
 . . doctest : :
 > > > from spectrum import lsf2poly
 > > > lsf = [ 0.7842 , 1.5605 , 1.8776 , 1.8984 , 2.3593]
 > > > a = lsf2poly ( lsf )
 # array ( [ 1.00000e + 00 , 6.14837835e - 01 , 9.89884967e - 01,
 # 9.31594056e - 05 , 3.13713832e - 03 , - 8.12002261e - 03 ] )
 . . seealso : : poly2lsf , rc2poly , ac2poly , rc2is"""

implementation:
# Reference : A . M . Kondoz , " Digital Speech : Coding for Low Bit Rate Communications
# Systems " John Wiley & Sons 1994 , Chapter 4
# Line spectral frequencies must be real .
lsf = numpy . array ( lsf )
if max ( lsf ) > numpy . pi or min ( lsf ) < 0 :
    raise ValueError ( 'Line spectral frequencies must be between 0 and pi.' )
p = len ( lsf )
# model order
# Form zeros using the LSFs and unit amplitudes
z = numpy . exp ( 1.j * lsf )
# Separate the zeros to those belonging to P and Q
rQ = z [ 0 : : 2 ]
rP = z [ 1 : : 2 ]
# Include the conjugates as well
rQ = numpy . concatenate ( ( rQ , rQ . conjugate ( ) ) )
rP = numpy . concatenate ( ( rP , rP . conjugate ( ) ) )
# Form the polynomials P and Q , note that these should be real
Q = numpy . poly ( rQ ) ;
P = numpy . poly ( rP ) ;
# Form the sum and difference filters by including known roots at z = 1 and
# z = - 1
if p % 2 : # Odd order : z = + 1 and z = - 1 are roots of the difference filter , P1 ( z )
    P1 = numpy . convolve ( P , [ 1 , 0 , - 1 ] )
    Q1 = Q
else : # Even order : z = - 1 is a root of the sum filter , Q1 ( z ) and z = 1 is a
# root of the difference filter , P1 ( z )
    P1 = numpy . convolve ( P , [ 1 , - 1 ] )
    Q1 = numpy . convolve ( Q , [ 1 , 1 ] )
# Prediction polynomial is formed by averaging P1 and Q1
a = .5 * ( P1 + Q1 )
return a [ 0 : - 1 : 1 ]


signature:
def poly2lsf ( a ) :
    """Prediction polynomial to line spectral frequencies .
 converts the prediction polynomial specified by A ,
 into the corresponding line spectral frequencies , LSF .
 normalizes the prediction polynomial by A ( 1 ) .
 . . doctest : :
 > > > from spectrum import poly2lsf
 > > > a = [ 1.0000 , 0.6149 , 0.9899 , 0.0000 , 0.0031 , - 0.0082]
 > > > lsf = poly2lsf ( a )
 > > > lsf = array ( [ 0.7842 , 1.5605 , 1.8776 , 1.8984 , 2.3593 ] )
 . . seealso : : lsf2poly , poly2rc , poly2qc , rc2is"""

implementation:
# Line spectral frequencies are not defined for complex polynomials .
# Normalize the polynomial
a = numpy . array ( a )
if a [ 0 ] != 1 :
    a /= a [ 0 ]
if max ( numpy . abs ( numpy . roots ( a ) ) ) >= 1.0 :
    error ( 'The polynomial must have all roots inside of the unit circle.' ) ;
# Form the sum and differnce filters
p = len ( a ) - 1
# The leading one in the polynomial is not used
a1 = numpy . concatenate ( ( a , numpy . array ( [ 0 ] ) ) )
a2 = a1 [ - 1 : : - 1 ]
P1 = a1 - a2
# Difference filter
Q1 = a1 + a2
# Sum Filter
# If order is even , remove the known root at z = 1 for P1 and z = - 1 for Q1
# If odd , remove both the roots from P1
if p % 2 : # Odd order
    P , r = deconvolve ( P1 , [ 1 , 0 , - 1 ] )
    Q = Q1
else : # Even order
    P , r = deconvolve ( P1 , [ 1 , - 1 ] )
    Q , r = deconvolve ( Q1 , [ 1 , 1 ] )
rP = numpy . roots ( P )
rQ = numpy . roots ( Q )
aP = numpy . angle ( rP [ 1 : : 2 ] )
aQ = numpy . angle ( rQ [ 1 : : 2 ] )
lsf = sorted ( numpy . concatenate ( ( - aP , - aQ ) ) )
return lsf


signature:
def _swapsides ( data ) :
    """todo is it really useful ?
 Swap sides
 . . doctest : :
 > > > from spectrum import swapsides
 > > > x = [ - 2 , - 1 , 1 , 2]
 > > > swapsides ( x )
 array ( [ 2 , - 2 , - 1 ] )"""

implementation:
N = len ( data )
return np . concatenate ( ( data [ N // 2 + 1 : ] , data [ 0 : N // 2 ] ) )


signature:
def twosided_2_onesided ( data ) :
    """Convert a one - sided PSD to a twosided PSD
 In order to keep the power in the onesided PSD the same
 as in the twosided version , the onesided values are twice
 as much as in the input data ( except for the zero - lag value ) .
 > > > twosided _ 2 _ onesided ( [ 10 , 2,3,3,2,8 ] )
 array ( [ 10 . , 4 . , 6 . , 8 . ] )"""

implementation:
assert len ( data ) % 2 == 0
N = len ( data )
psd = np . array ( data [ 0 : N // 2 + 1 ] ) * 2.
psd [ 0 ] /= 2.
psd [ - 1 ] = data [ - 1 ]
return psd


signature:
def onesided_2_twosided ( data ) :
    """Convert a two - sided PSD to a one - sided PSD
 In order to keep the power in the twosided PSD the same
 as in the onesided version , the twosided values are 2 times
 lower than the input data ( except for the zero - lag and N - lag
 values ) .
 > > > twosided _ 2 _ onesided ( [ 10 , 4 , 6 , 8 ] )
 array ( [ 10 . , 2 . , 3 . , 3 . , 2 . , 8 . ] )"""

implementation:
psd = np . concatenate ( ( data [ 0 : - 1 ] , cshift ( data [ - 1 : 0 : - 1 ] , - 1 ) ) ) / 2.
psd [ 0 ] *= 2.
psd [ - 1 ] *= 2.
return psd


signature:
def twosided_2_centerdc ( data ) :
    """Convert a two - sided PSD to a center - dc PSD"""

implementation:
N = len ( data )
# could us int ( ) or / / in python 3
newpsd = np . concatenate ( ( cshift ( data [ N // 2 : ] , 1 ) , data [ 0 : N // 2 ] ) )
newpsd [ 0 ] = data [ - 1 ]
return newpsd


signature:
def centerdc_2_twosided ( data ) :
    """Convert a center - dc PSD to a twosided PSD"""

implementation:
N = len ( data )
newpsd = np . concatenate ( ( data [ N // 2 : ] , ( cshift ( data [ 0 : N // 2 ] , - 1 ) ) ) )
return newpsd


signature:
def _twosided_zerolag ( data , zerolag ) :
    """Build a symmetric vector out of stricly positive lag vector and zero - lag
 . . doctest : :
 > > > data = [ 3,2,1]
 > > > zerolag = 4
 > > > twosided _ zerolag ( data , zerolag )
 array ( [ 1 , 2 , 3 , 4 , 3 , 2 , 1 ] )
 . . seealso : : Same behaviour as : func : ` twosided _ zerolag `"""

implementation:
res = twosided ( np . insert ( data , 0 , zerolag ) )
return res


signature:
def cshift ( data , offset ) :
    """Circular shift to the right ( within an array ) by a given offset
 : param array data : input data ( list or numpy . array )
 : param int offset : shift the array with the offset
 . . doctest : :
 > > > from spectrum import cshift
 > > > cshift ( [ 0 , 1 , 2 , 3 , - 2 , - 1 ] , 2)
 array ( [ - 2 , - 1 , 0 , 1 , 2 , 3 ] )"""

implementation:
# the deque method is suppose to be optimal when using rotate to shift the
# data that playing with the data to build a new list .
if isinstance ( offset , float ) :
    offset = int ( offset )
a = deque ( data )
a . rotate ( offset )
return np . array ( a )


signature:
def data_cosine ( N = 1024 , A = 0.1 , sampling = 1024. , freq = 200 ) :
    r"""Return a noisy cosine at a given frequency .
 : param N : the final data size
 : param A : the strength of the noise
 : param float sampling : sampling frequency of the input : attr : ` data ` .
 : param float freq : the frequency : math : ` f _ 0 ` of the cosine .
 . . math : : x [ t ] = cos ( 2 \ pi t * f _ 0 ) + A w [ t ]
 where w [ t ] is a white noise of variance 1.
 . . doctest : :
 > > > from spectrum import data _ cosine
 > > > a = data _ cosine ( N = 1024 , sampling = 1024 , A = 0.5 , freq = 100)"""

implementation:
t = arange ( 0 , float ( N ) / sampling , 1. / sampling )
x = cos ( 2. * pi * t * freq ) + A * randn ( t . size )
return x


signature:
def data_two_freqs ( N = 200 ) :
    """A simple test example with two close frequencies"""

implementation:
nn = arange ( N )
xx = cos ( 0.257 * pi * nn ) + sin ( 0.2 * pi * nn ) + 0.01 * randn ( nn . size )
return xx


signature:
def spectrum_data ( filename ) :
    """Simple utilities to retrieve data sets from"""

implementation:
import os
import pkg_resources
info = pkg_resources . get_distribution ( 'spectrum' )
location = info . location
# first try develop mode
share = os . sep . join ( [ location , "spectrum" , 'data' ] )
filename2 = os . sep . join ( [ share , filename ] )
if os . path . exists ( filename2 ) :
    return filename2
else :
    raise Exception ( 'unknown file %s' % filename2 )


signature:
def plot ( self , ** kargs ) :
    """Plot the data set , using the sampling information to set the x - axis
 correctly ."""

implementation:
from pylab import plot , linspace , xlabel , ylabel , grid
time = linspace ( 1 * self . dt , self . N * self . dt , self . N )
plot ( time , self . data , ** kargs )
xlabel ( 'Time' )
ylabel ( 'Amplitude' )
grid ( True )


signature:
def readwav ( filename ) :
    """Read a WAV file and returns the data and sample rate
 from spectrum . io import readwav
 readwav ( )"""

implementation:
from scipy . io . wavfile import read as readwav
samplerate , signal = readwav ( filename )
return signal , samplerate


signature:
def pmtm ( x , NW = None , k = None , NFFT = None , e = None , v = None , method = 'adapt' , show = False ) :
    """Multitapering spectral estimation
 : param array x : the data
 : param float NW : The time half bandwidth parameter ( typical values are
 2.5,3,3.5,4 ) . Must be provided otherwise the tapering windows and
 eigen values ( outputs of dpss ) must be provided
 : param int k : uses the first k Slepian sequences . If * k * is not provided ,
 * k * is set to * NW * 2 * .
 : param NW :
 : param e : the window concentrations ( eigenvalues )
 : param v : the matrix containing the tapering windows
 : param str method : set how the eigenvalues are used . Must be
 in [ ' unity ' , ' adapt ' , ' eigen ' ]
 : param bool show : plot results
 : return : Sk ( complex ) , weights , eigenvalues
 Usually in spectral estimation the mean to reduce bias is to use tapering
 window . In order to reduce variance we need to average different spectrum .
 The problem is that we have only one set of data . Thus we need to
 decompose a set into several segments . Such method are well - known : simple
 daniell ' s periodogram , Welch ' s method and so on . The drawback of such
 methods is a loss of resolution since the segments used to compute the
 spectrum are smaller than the data set .
 The interest of multitapering method is to keep a good resolution while
 reducing bias and variance .
 How does it work ? First we compute different simple periodogram with the
 whole data set ( to keep good resolution ) but each periodgram is computed
 with a differenttapering windows . Then , we average all these spectrum .
 To avoid redundancy and bias due to the tapers mtm use special tapers .
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import data _ cosine , dpss , pmtm
 data = data _ cosine ( N = 2048 , A = 0.1 , sampling = 1024 , freq = 200)
 # If you already have the DPSS windows
 [ tapers , eigen ] = dpss ( 2048 , 2.5 , 4)
 res = pmtm ( data , e = eigen , v = tapers , show = False )
 # You do not need to compute the DPSS before end
 res = pmtm ( data , NW = 2.5 , show = False )
 res = pmtm ( data , NW = 2.5 , k = 4 , show = True )
 . . versionchanged : : 0.6.2
 APN modified method to return each Sk as complex values , the eigenvalues
 and the weights"""

implementation:
assert method in [ 'adapt' , 'eigen' , 'unity' ]
N = len ( x )
# if dpss not provided , compute them
if e is None and v is None :
    if NW is not None :
        [ tapers , eigenvalues ] = dpss ( N , NW , k = k )
    else :
        raise ValueError ( "NW must be provided (e.g. 2.5, 3, 3.5, 4" )
elif e is not None and v is not None :
    eigenvalues = e [ : ]
    tapers = v [ : ]
else :
    raise ValueError ( "if e provided, v must be provided as well and viceversa." )
nwin = len ( eigenvalues )
# length of the eigen values vector to be used later
# set the NFFT
if NFFT == None :
    NFFT = max ( 256 , 2 ** nextpow2 ( N ) )
Sk_complex = np . fft . fft ( np . multiply ( tapers . transpose ( ) , x ) , NFFT )
Sk = abs ( Sk_complex ) ** 2
# si nfft smaller thqn N , cut otherwise add wero .
# compute
if method in [ 'eigen' , 'unity' ] :
    if method == 'unity' :
        weights = np . ones ( ( nwin , 1 ) )
    elif method == 'eigen' : # The S _ k spectrum can be weighted by the eigenvalues , as in Park et al .
        weights = np . array ( [ _x / float ( i + 1 ) for i , _x in enumerate ( eigenvalues ) ] )
        weights = weights . reshape ( nwin , 1 )
elif method == 'adapt' : # This version uses the equations from [ 2 ] ( P & W pp 368-370 ) .
# Wrap the data modulo nfft if N > nfft
    sig2 = np . dot ( x , x ) / float ( N )
    Sk = abs ( np . fft . fft ( np . multiply ( tapers . transpose ( ) , x ) , NFFT ) ) ** 2
    Sk = Sk . transpose ( )
    S = ( Sk [ : , 0 ] + Sk [ : , 1 ] ) / 2
    # Initial spectrum estimate
    S = S . reshape ( NFFT , 1 )
    Stemp = np . zeros ( ( NFFT , 1 ) )
    S1 = np . zeros ( ( NFFT , 1 ) )
    # Set tolerance for acceptance of spectral estimate :
    tol = 0.0005 * sig2 / float ( NFFT )
    i = 0
    a = sig2 * ( 1 - eigenvalues )
    # converges very quickly but for safety ; set i < 100
    while sum ( np . abs ( S - S1 ) ) / NFFT > tol and i < 100 :
        i = i + 1
        # calculate weights
        b1 = np . multiply ( S , np . ones ( ( 1 , nwin ) ) )
        b2 = np . multiply ( S , eigenvalues . transpose ( ) ) + np . ones ( ( NFFT , 1 ) ) * a . transpose ( )
        b = b1 / b2
        # calculate new spectral estimate
        wk = ( b ** 2 ) * ( np . ones ( ( NFFT , 1 ) ) * eigenvalues . transpose ( ) )
        S1 = sum ( wk . transpose ( ) * Sk . transpose ( ) ) / sum ( wk . transpose ( ) )
        S1 = S1 . reshape ( NFFT , 1 )
        Stemp = S1
        S1 = S
        S = Stemp
        # swap S and S1
    weights = wk
if show is True :
    from pylab import semilogy
    if method == "adapt" :
        Sk = np . mean ( Sk * weights , axis = 1 )
    else :
        Sk = np . mean ( Sk * weights , axis = 0 )
    semilogy ( Sk )
return Sk_complex , weights , eigenvalues


signature:
def dpss ( N , NW = None , k = None ) :
    r"""Discrete prolate spheroidal ( Slepian ) sequences
 Calculation of the Discrete Prolate Spheroidal Sequences also known as the
 slepian sequences , and the corresponding eigenvalues .
 : param int N : desired window length
 : param float NW : The time half bandwidth parameter ( typical values are
 2.5,3,3.5,4 ) .
 : param int k : returns the first k Slepian sequences . If * k * is not
 provided , * k * is set to * NW * 2 * .
 : return :
 * tapers , a matrix of tapering windows . Matrix is a N by * k * ( k
 is the number of windows )
 * eigen , a vector of eigenvalues of length * k *
 The discrete prolate spheroidal or Slepian sequences derive from the following
 time - frequency concentration problem . For all finite - energy sequences index
 limited to some set , which sequence maximizes the following ratio :
 . . math : :
 \ lambda = \ frac { \ int _ { - W } ^ { W } \ left | X ( f ) \ right | ^ 2 df }
 { \ int _ { - F _ s / 2 } ^ { F _ s / 2 } \ left | X ( f ) \ right | ^ 2 df }
 where : math : ` F _ s ` is the sampling frequency and : math : ` | W | < F _ s / 2 ` .
 This ratio determines which index - limited sequence has the largest proportion of its
 energy in the band : math : ` [ - W , W ] ` with : math : ` 0 < \ lambda < 1 ` .
 The sequence maximizing the ratio is the first
 discrete prolate spheroidal or Slepian sequence . The second Slepian sequence
 maximizes the ratio and is orthogonal to the first Slepian sequence . The third
 Slepian sequence maximizes the ratio of integrals and is orthogonal to both
 the first and second Slepian sequences and so on .
 . . note : : Note about the implementation . Since the slepian generation is
 computationally expensive , we use a C implementation based on the C
 code written by Lees as published in :
 Lees , J . M . and J . Park ( 1995 ) : Multiple - taper spectral analysis : A stand - alone
 C - subroutine : Computers & Geology : 21 , 199-236.
 However , the original C code has been trimmed . Indeed , we only require the
 multitap function ( that depends on jtridib , jtinvit functions only ) .
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import *
 from pylab import *
 N = 512
 [ w , eigens ] = dpss ( N , 2.5 , 4)
 plot ( w )
 title ( ' Slepian Sequences N = % s , NW = 2.5 ' % N )
 axis ( [ 0 , N , - 0.15 , 0.15 ] )
 legend ( [ ' 1st window ' , ' 2nd window ' , ' 3rd window ' , ' 4th window ' ] )
 Windows are normalised :
 . . math : : \ sum _ k h _ k h _ k = 1
 : references : [ Percival ] _
 Slepian , D . Prolate spheroidal wave functions , Fourier analysis , and
 uncertainty V : The discrete case . Bell System Technical Journal ,
 Volume 57 ( 1978 ) , 1371430
 . . note : : the C code to create the slepian windows is extracted from original C code
 from Lees and Park ( 1995 ) and uses the conventions of Percival and Walden ( 1993 ) .
 Functions that are not used here were removed ."""

implementation:
assert NW < N / 2 , "NW ({}) must be stricly less than N/2 ({}/2)" . format ( NW , N )
if k is None :
    k = min ( round ( 2 * NW ) , N )
    k = int ( max ( k , 1 ) )
from numpy import dot , zeros , arange , sqrt
mtspeclib . multitap . restype = None
lam = zeros ( k , dtype = float )
tapers = zeros ( k * N , dtype = float )
tapsum = zeros ( k , dtype = float )
res = mtspeclib . multitap ( c_int ( N ) , c_int ( k ) , lam . ctypes . data_as ( c_void_p ) , c_float ( NW ) , tapers . ctypes . data_as ( c_void_p ) , tapsum . ctypes . data_as ( c_void_p ) , )
# normalisation by sqtr ( N ) . It is required to have normalised windows
tapers = tapers . reshape ( k , N ) . transpose ( ) / sqrt ( N )
for i in range ( k ) : # By convention ( Percival and Walden , 1993 pg 379)
# * symmetric tapers ( k = 0,2,4 , . . . ) should have a positive average .
# * antisymmetric tapers should begin with a positive lobe
    if i % 2 == 0 :
        if tapsum [ i ] < 0 :
            tapsum [ i ] *= - 1
            tapers [ : , i ] *= - 1
    else :
        if tapers [ 0 , i ] < 0 :
            tapsum [ i ] *= - 1
            tapers [ : , i ] *= - 1
# Now find the eigenvalues of the original
# Use the autocovariance sequence technique from Percival and Walden , 1993
# pg 390 to get the eigenvalues more precisely ( same as matlab output )
# The values returned in lam are not exacly the same as in the following methods .
acvs = _autocov ( tapers . transpose ( ) , debias = False ) * N
nidx = arange ( N )
W = float ( NW ) / N
r = 4 * W * np . sinc ( 2 * W * nidx )
r [ 0 ] = 2 * W
eigvals = dot ( acvs , r )
# return ( tapers , lam )
return [ tapers , eigvals ]


signature:
def _other_dpss_method ( N , NW , Kmax ) :
    """Returns the Discrete Prolate Spheroidal Sequences of orders [ 0 , Kmax - 1]
 for a given frequency - spacing multiple NW and sequence length N .
 See dpss function that is the official version . This version is indepedant
 of the C code and relies on Scipy function . However , it is slower by a factor 3
 Tridiagonal form of DPSS calculation from :"""

implementation:
# here we want to set up an optimization problem to find a sequence
# whose energy is maximally concentrated within band [ - W , W ] .
# Thus , the measure lambda ( T , W ) is the ratio between the energy within
# that band , and the total energy . This leads to the eigen - system
# ( A - ( l1 ) I ) v = 0 , where the eigenvector corresponding to the largest
# eigenvalue is the sequence with maximally concentrated energy . The
# collection of eigenvectors of this system are called Slepian sequences ,
# or discrete prolate spheroidal sequences ( DPSS ) . Only the first K ,
# K = 2NW / dt orders of DPSS will exhibit good spectral concentration
# [ see http : / / en . wikipedia . org / wiki / Spectral _ concentration _ problem ]
# Here I set up an alternative symmetric tri - diagonal eigenvalue problem
# such that
# ( B - ( l2 ) I ) v = 0 , and v are our DPSS ( but eigenvalues l2 ! = l1)
# the main diagonal = ( [ N - 1-2 * t ] / 2 ) * * 2 cos ( 2PIW ) , t = [ 0,1,2 , . . . , N - 1]
# and the first off - diangonal = t ( N - t ) / 2 , t = [ 1,2 , . . . , N - 1]
# [ see Percival and Walden , 1993]
from scipy import linalg as la
Kmax = int ( Kmax )
W = float ( NW ) / N
ab = np . zeros ( ( 2 , N ) , 'd' )
nidx = np . arange ( N )
ab [ 0 , 1 : ] = nidx [ 1 : ] * ( N - nidx [ 1 : ] ) / 2.
ab [ 1 ] = ( ( N - 1 - 2 * nidx ) / 2. ) ** 2 * np . cos ( 2 * np . pi * W )
# only calculate the highest Kmax - 1 eigenvectors
l , v = la . eig_banded ( ab , select = 'i' , select_range = ( N - Kmax , N - 1 ) )
dpss = v . transpose ( ) [ : : - 1 ]
# By convention ( Percival and Walden , 1993 pg 379)
# * symmetric tapers ( k = 0,2,4 , . . . ) should have a positive average .
# * antisymmetric tapers should begin with a positive lobe
fix_symmetric = ( dpss [ 0 : : 2 ] . sum ( axis = 1 ) < 0 )
for i , f in enumerate ( fix_symmetric ) :
    if f :
        dpss [ 2 * i ] *= - 1
fix_skew = ( dpss [ 1 : : 2 , 1 ] < 0 )
for i , f in enumerate ( fix_skew ) :
    if f :
        dpss [ 2 * i + 1 ] *= - 1
# Now find the eigenvalues of the original
# Use the autocovariance sequence technique from Percival and Walden , 1993
# pg 390
# XXX : why debias false ? it ' s all messed up o . w . , even with means
# on the order of 1e - 2
acvs = _autocov ( dpss , debias = False ) * N
r = 4 * W * np . sinc ( 2 * W * nidx )
r [ 0 ] = 2 * W
eigvals = np . dot ( acvs , r )
return dpss , eigvals


signature:
def _autocov ( s , ** kwargs ) :
    """Returns the autocovariance of signal s at all lags .
 Adheres to the definition
 sxx [ k ] = E { S [ n ] S [ n + k ] } = cov { S [ n ] , S [ n + k ] }
 where E { } is the expectation operator , and S is a zero mean process"""

implementation:
# only remove the mean once , if needed
debias = kwargs . pop ( 'debias' , True )
axis = kwargs . get ( 'axis' , - 1 )
if debias :
    s = _remove_bias ( s , axis )
kwargs [ 'debias' ] = False
return _crosscov ( s , s , ** kwargs )


signature:
def _crosscov ( x , y , axis = - 1 , all_lags = False , debias = True ) :
    """Returns the crosscovariance sequence between two ndarrays .
 This is performed by calling fftconvolve on x , y [ : : - 1]
 Parameters
 x : ndarray
 y : ndarray
 axis : time axis
 all _ lags : { True / False }
 whether to return all nonzero lags , or to clip the length of s _ xy
 to be the length of x and y . If False , then the zero lag covariance
 is at index 0 . Otherwise , it is found at ( len ( x ) + len ( y ) - 1 ) / 2
 debias : { True / False }
 Always removes an estimate of the mean along the axis , unless
 told not to .
 cross covariance is defined as
 sxy [ k ] : = E { X [ t ] * Y [ t + k ] } , where X , Y are zero mean random processes"""

implementation:
if x . shape [ axis ] != y . shape [ axis ] :
    raise ValueError ( 'crosscov() only works on same-length sequences for now' )
if debias :
    x = _remove_bias ( x , axis )
    y = _remove_bias ( y , axis )
slicing = [ slice ( d ) for d in x . shape ]
slicing [ axis ] = slice ( None , None , - 1 )
sxy = _fftconvolve ( x , y [ tuple ( slicing ) ] , axis = axis , mode = 'full' )
N = x . shape [ axis ]
sxy /= N
if all_lags :
    return sxy
slicing [ axis ] = slice ( N - 1 , 2 * N - 1 )
return sxy [ tuple ( slicing ) ]


signature:
def _crosscorr ( x , y , ** kwargs ) :
    """Returns the crosscorrelation sequence between two ndarrays .
 This is performed by calling fftconvolve on x , y [ : : - 1]
 Parameters
 x : ndarray
 y : ndarray
 axis : time axis
 all _ lags : { True / False }
 whether to return all nonzero lags , or to clip the length of r _ xy
 to be the length of x and y . If False , then the zero lag correlation
 is at index 0 . Otherwise , it is found at ( len ( x ) + len ( y ) - 1 ) / 2
 Notes
 cross correlation is defined as
 rxy [ k ] : = E { X [ t ] * Y [ t + k ] } / ( E { X * X } E { Y * Y } ) * * . 5,
 where X , Y are zero mean random processes . It is the noramlized cross
 covariance ."""

implementation:
sxy = _crosscov ( x , y , ** kwargs )
# estimate sigma _ x , sigma _ y to normalize
sx = np . std ( x )
sy = np . std ( y )
return sxy / ( sx * sy )


signature:
def _remove_bias ( x , axis ) :
    "Subtracts an estimate of the mean from signal x at axis"

implementation:
padded_slice = [ slice ( d ) for d in x . shape ]
padded_slice [ axis ] = np . newaxis
mn = np . mean ( x , axis = axis )
return x - mn [ tuple ( padded_slice ) ]


signature:
def get_docstring_and_rest ( filename ) :
    """Separate ` filename ` content between docstring and the rest
 Strongly inspired from ast . get _ docstring .
 Returns
 docstring : str
 docstring of ` filename `
 rest : str
 ` filename ` content without the docstring"""

implementation:
with open ( filename ) as f :
    content = f . read ( )
node = ast . parse ( content )
if not isinstance ( node , ast . Module ) :
    raise TypeError ( "This function only supports modules. " "You provided {0}" . format ( node . __class__ . __name__ ) )
if node . body and isinstance ( node . body [ 0 ] , ast . Expr ) and isinstance ( node . body [ 0 ] . value , ast . Str ) :
    docstring_node = node . body [ 0 ]
    docstring = docstring_node . value . s
    # This get the content of the file after the docstring last line
    # Note : ' maxsplit ' argument is not a keyword argument in python2
    rest = content . split ( '\n' , docstring_node . lineno ) [ - 1 ]
    return docstring , rest
else :
    raise ValueError ( ( 'Could not find docstring in file "{0}". ' 'A docstring is required by sphinx-gallery' ) . format ( filename ) )


signature:
def split_code_and_text_blocks ( source_file ) :
    """Return list with source file separated into code and text blocks .
 Returns
 blocks : list of ( label , content )
 List where each element is a tuple with the label ( ' text ' or ' code ' ) ,
 and content string of block ."""

implementation:
docstring , rest_of_content = get_docstring_and_rest ( source_file )
blocks = [ ( 'text' , docstring ) ]
pattern = re . compile ( r'(?P<header_line>^#{20,}.*)\s(?P<text_content>(?:^#.*\s)*)' , flags = re . M )
pos_so_far = 0
for match in re . finditer ( pattern , rest_of_content ) :
    match_start_pos , match_end_pos = match . span ( )
    code_block_content = rest_of_content [ pos_so_far : match_start_pos ]
    text_content = match . group ( 'text_content' )
    sub_pat = re . compile ( '^#' , flags = re . M )
    text_block_content = dedent ( re . sub ( sub_pat , '' , text_content ) )
    if code_block_content . strip ( ) :
        blocks . append ( ( 'code' , code_block_content ) )
    if text_block_content . strip ( ) :
        blocks . append ( ( 'text' , text_block_content ) )
    pos_so_far = match_end_pos
remaining_content = rest_of_content [ pos_so_far : ]
if remaining_content . strip ( ) :
    blocks . append ( ( 'code' , remaining_content ) )
return blocks


signature:
def codestr2rst ( codestr , lang = 'python' ) :
    """Return reStructuredText code block from code string"""

implementation:
code_directive = "\n.. code-block:: {0}\n\n" . format ( lang )
indented_block = indent ( codestr , ' ' * 4 )
return code_directive + indented_block


signature:
def extract_intro ( filename ) :
    """Extract the first paragraph of module - level docstring . max : 95 char"""

implementation:
docstring , _ = get_docstring_and_rest ( filename )
# lstrip is just in case docstring has a ' \ n \ n ' at the beginning
paragraphs = docstring . lstrip ( ) . split ( '\n\n' )
if len ( paragraphs ) > 1 :
    first_paragraph = re . sub ( '\n' , ' ' , paragraphs [ 1 ] )
    first_paragraph = ( first_paragraph [ : 95 ] + '...' if len ( first_paragraph ) > 95 else first_paragraph )
else :
    raise ValueError ( "Example docstring should have a header for the example title " "and at least a paragraph explaining what the example is about. " "Please check the example file:\n {}\n" . format ( filename ) )
return first_paragraph


signature:
def get_md5sum ( src_file ) :
    """Returns md5sum of file"""

implementation:
with open ( src_file , 'r' ) as src_data :
    src_content = src_data . read ( )
    # data needs to be encoded in python3 before hashing
    if sys . version_info [ 0 ] == 3 :
        src_content = src_content . encode ( 'utf-8' )
    src_md5 = hashlib . md5 ( src_content ) . hexdigest ( )
return src_md5


signature:
def check_md5sum_change ( src_file ) :
    """Returns True if src _ file has a different md5sum"""

implementation:
src_md5 = get_md5sum ( src_file )
src_md5_file = src_file + '.md5'
src_file_changed = True
if os . path . exists ( src_md5_file ) :
    with open ( src_md5_file , 'r' ) as file_checksum :
        ref_md5 = file_checksum . read ( )
    if src_md5 == ref_md5 :
        src_file_changed = False
if src_file_changed :
    with open ( src_md5_file , 'w' ) as file_checksum :
        file_checksum . write ( src_md5 )
return src_file_changed


signature:
def _plots_are_current ( src_file , image_file ) :
    """Test existence of image file and no change in md5sum of
 example"""

implementation:
first_image_file = image_file . format ( 1 )
has_image = os . path . exists ( first_image_file )
src_file_changed = check_md5sum_change ( src_file )
return has_image and not src_file_changed


signature:
def save_figures ( image_path , fig_count , gallery_conf ) :
    """Save all open matplotlib figures of the example code - block
 Parameters
 image _ path : str
 Path where plots are saved ( format string which accepts figure number )
 fig _ count : int
 Previous figure number count . Figure number add from this number
 Returns
 list of strings containing the full path to each figure"""

implementation:
figure_list = [ ]
fig_managers = matplotlib . _pylab_helpers . Gcf . get_all_fig_managers ( )
for fig_mngr in fig_managers : # Set the fig _ num figure as the current figure as we can ' t
# save a figure that ' s not the current figure .
    fig = plt . figure ( fig_mngr . num )
    kwargs = { }
    to_rgba = matplotlib . colors . colorConverter . to_rgba
    for attr in [ 'facecolor' , 'edgecolor' ] :
        fig_attr = getattr ( fig , 'get_' + attr ) ( )
        default_attr = matplotlib . rcParams [ 'figure.' + attr ]
        if to_rgba ( fig_attr ) != to_rgba ( default_attr ) :
            kwargs [ attr ] = fig_attr
    current_fig = image_path . format ( fig_count + fig_mngr . num )
    fig . savefig ( current_fig , ** kwargs )
    figure_list . append ( current_fig )
if gallery_conf . get ( 'find_mayavi_figures' , False ) :
    from mayavi import mlab
    e = mlab . get_engine ( )
    last_matplotlib_fig_num = len ( figure_list )
    total_fig_num = last_matplotlib_fig_num + len ( e . scenes )
    mayavi_fig_nums = range ( last_matplotlib_fig_num , total_fig_num )
    for scene , mayavi_fig_num in zip ( e . scenes , mayavi_fig_nums ) :
        current_fig = image_path . format ( mayavi_fig_num )
        mlab . savefig ( current_fig , figure = scene )
        # make sure the image is not too large
        scale_image ( current_fig , current_fig , 850 , 999 )
        figure_list . append ( current_fig )
    mlab . close ( all = True )
return figure_list


signature:
def scale_image ( in_fname , out_fname , max_width , max_height ) :
    """Scales an image with the same aspect ratio centered in an
 image with a given max _ width and max _ height
 if in _ fname = = out _ fname the image can only be scaled down"""

implementation:
# local import to avoid testing dependency on PIL :
try :
    from PIL import Image
except ImportError :
    import Image
img = Image . open ( in_fname )
width_in , height_in = img . size
scale_w = max_width / float ( width_in )
scale_h = max_height / float ( height_in )
if height_in * scale_w <= max_height :
    scale = scale_w
else :
    scale = scale_h
if scale >= 1.0 and in_fname == out_fname :
    return
width_sc = int ( round ( scale * width_in ) )
height_sc = int ( round ( scale * height_in ) )
# resize the image
img . thumbnail ( ( width_sc , height_sc ) , Image . ANTIALIAS )
# insert centered
thumb = Image . new ( 'RGB' , ( max_width , max_height ) , ( 255 , 255 , 255 ) )
pos_insert = ( ( max_width - width_sc ) // 2 , ( max_height - height_sc ) // 2 )
thumb . paste ( img , pos_insert )
thumb . save ( out_fname )
# Use optipng to perform lossless compression on the resized image if
# software is installed
if os . environ . get ( 'SKLEARN_DOC_OPTIPNG' , False ) :
    try :
        subprocess . call ( [ "optipng" , "-quiet" , "-o" , "9" , out_fname ] )
    except Exception :
        warnings . warn ( 'Install optipng to reduce the size of the \
                          generated images' )


signature:
def save_thumbnail ( image_path , base_image_name , gallery_conf ) :
    """Save the thumbnail image"""

implementation:
first_image_file = image_path . format ( 1 )
thumb_dir = os . path . join ( os . path . dirname ( first_image_file ) , 'thumb' )
if not os . path . exists ( thumb_dir ) :
    os . makedirs ( thumb_dir )
thumb_file = os . path . join ( thumb_dir , 'sphx_glr_%s_thumb.png' % base_image_name )
if os . path . exists ( first_image_file ) :
    scale_image ( first_image_file , thumb_file , 400 , 280 )
elif not os . path . exists ( thumb_file ) : # create something to replace the thumbnail
    default_thumb_file = os . path . join ( glr_path_static ( ) , 'no_image.png' )
    default_thumb_file = gallery_conf . get ( "default_thumb_file" , default_thumb_file )
    scale_image ( default_thumb_file , thumb_file , 200 , 140 )


signature:
def generate_dir_rst ( src_dir , target_dir , gallery_conf , seen_backrefs ) :
    """Generate the gallery reStructuredText for an example directory"""

implementation:
if not os . path . exists ( os . path . join ( src_dir , 'README.txt' ) ) :
    print ( 80 * '_' )
    print ( 'Example directory %s does not have a README.txt file' % src_dir )
    print ( 'Skipping this directory' )
    print ( 80 * '_' )
    return ""
# because string is an expected return type
fhindex = open ( os . path . join ( src_dir , 'README.txt' ) ) . read ( )
if not os . path . exists ( target_dir ) :
    os . makedirs ( target_dir )
sorted_listdir = [ fname for fname in sorted ( os . listdir ( src_dir ) ) if fname . endswith ( '.py' ) ]
entries_text = [ ]
for fname in sorted_listdir :
    amount_of_code = generate_file_rst ( fname , target_dir , src_dir , gallery_conf )
    new_fname = os . path . join ( src_dir , fname )
    intro = extract_intro ( new_fname )
    write_backreferences ( seen_backrefs , gallery_conf , target_dir , fname , intro )
    this_entry = _thumbnail_div ( target_dir , fname , intro ) + """

.. toctree::
   :hidden:

   /%s/%s\n""" % ( target_dir , fname [ : - 3 ] )
    entries_text . append ( ( amount_of_code , this_entry ) )
# sort to have the smallest entries in the beginning
entries_text . sort ( )
for _ , entry_text in entries_text :
    fhindex += entry_text
# clear at the end of the section
fhindex += """.. raw:: html\n

\n\n"""
return fhindex


signature:
def execute_script ( code_block , example_globals , image_path , fig_count , src_file , gallery_conf ) :
    """Executes the code block of the example file"""

implementation:
time_elapsed = 0
stdout = ''
# We need to execute the code
print ( 'plotting code blocks in %s' % src_file )
plt . close ( 'all' )
cwd = os . getcwd ( )
# Redirect output to stdout and
orig_stdout = sys . stdout
try : # First cd in the original example dir , so that any file
# created by the example get created in this directory
    os . chdir ( os . path . dirname ( src_file ) )
    my_buffer = StringIO ( )
    my_stdout = Tee ( sys . stdout , my_buffer )
    sys . stdout = my_stdout
    t_start = time ( )
    exec ( code_block , example_globals )
    time_elapsed = time ( ) - t_start
    sys . stdout = orig_stdout
    my_stdout = my_buffer . getvalue ( ) . strip ( ) . expandtabs ( )
    if my_stdout :
        stdout = CODE_OUTPUT . format ( indent ( my_stdout , ' ' * 4 ) )
    os . chdir ( cwd )
    figure_list = save_figures ( image_path , fig_count , gallery_conf )
    # Depending on whether we have one or more figures , we ' re using a
    # horizontal list or a single rst call to ' image ' .
    image_list = ""
    if len ( figure_list ) == 1 :
        figure_name = figure_list [ 0 ]
        image_list = SINGLE_IMAGE % figure_name . lstrip ( '/' )
    elif len ( figure_list ) > 1 :
        image_list = HLIST_HEADER
        for figure_name in figure_list :
            image_list += HLIST_IMAGE_TEMPLATE % figure_name . lstrip ( '/' )
except Exception :
    formatted_exception = traceback . format_exc ( )
    print ( 80 * '_' )
    print ( '%s is not compiling:' % src_file )
    print ( formatted_exception )
    print ( 80 * '_' )
    figure_list = [ ]
    image_list = codestr2rst ( formatted_exception , lang = 'pytb' )
    # Overrides the output thumbnail in the gallery for easy identification
    broken_img = os . path . join ( glr_path_static ( ) , 'broken_example.png' )
    shutil . copyfile ( broken_img , os . path . join ( cwd , image_path . format ( 1 ) ) )
    fig_count += 1
    # raise count to avoid overwriting image
    # Breaks build on first example error
    if gallery_conf [ 'abort_on_example_error' ] :
        raise
finally :
    os . chdir ( cwd )
    sys . stdout = orig_stdout
print ( " - time elapsed : %.2g sec" % time_elapsed )
code_output = "\n{0}\n\n{1}\n\n" . format ( image_list , stdout )
return code_output , time_elapsed , fig_count + len ( figure_list )


signature:
def generate_file_rst ( fname , target_dir , src_dir , gallery_conf ) :
    """Generate the rst file for a given example .
 Returns the amout of code ( in characters ) of the corresponding
 files ."""

implementation:
src_file = os . path . join ( src_dir , fname )
example_file = os . path . join ( target_dir , fname )
shutil . copyfile ( src_file , example_file )
image_dir = os . path . join ( target_dir , 'images' )
if not os . path . exists ( image_dir ) :
    os . makedirs ( image_dir )
base_image_name = os . path . splitext ( fname ) [ 0 ]
image_fname = 'sphx_glr_' + base_image_name + '_{0:03}.png'
image_path = os . path . join ( image_dir , image_fname )
script_blocks = split_code_and_text_blocks ( example_file )
amount_of_code = sum ( [ len ( bcontent ) for blabel , bcontent in script_blocks if blabel == 'code' ] )
if _plots_are_current ( example_file , image_path ) :
    return amount_of_code
time_elapsed = 0
ref_fname = example_file . replace ( os . path . sep , '_' )
example_rst = """\n\n.. _sphx_glr_{0}:\n\n""" . format ( ref_fname )
example_nb = Notebook ( fname , target_dir )
filename_pattern = gallery_conf . get ( 'filename_pattern' )
if re . search ( filename_pattern , src_file ) and gallery_conf [ 'plot_gallery' ] : # A lot of examples contains ' print ( _ _ doc _ _ ) ' for example in
# scikit - learn so that running the example prints some useful
# information . Because the docstring has been separated from
# the code blocks in sphinx - gallery , _ _ doc _ _ is actually
# _ _ builtin _ _ . _ _ doc _ _ in the execution context and we do not
# want to print it
    example_globals = { '__doc__' : '' }
    fig_count = 0
    # A simple example has two blocks : one for the
    # example introduction / explanation and one for the code
    is_example_notebook_like = len ( script_blocks ) > 2
    for blabel , bcontent in script_blocks :
        if blabel == 'code' :
            code_output , rtime , fig_count = execute_script ( bcontent , example_globals , image_path , fig_count , src_file , gallery_conf )
            time_elapsed += rtime
            example_nb . add_code_cell ( bcontent )
            if is_example_notebook_like :
                example_rst += codestr2rst ( bcontent ) + '\n'
                example_rst += code_output
            else :
                example_rst += code_output
                example_rst += codestr2rst ( bcontent ) + '\n'
        else :
            example_rst += text2string ( bcontent ) + '\n'
            example_nb . add_markdown_cell ( text2string ( bcontent ) )
else :
    for blabel , bcontent in script_blocks :
        if blabel == 'code' :
            example_rst += codestr2rst ( bcontent ) + '\n'
            example_nb . add_code_cell ( bcontent )
        else :
            example_rst += bcontent + '\n'
            example_nb . add_markdown_cell ( text2string ( bcontent ) )
save_thumbnail ( image_path , base_image_name , gallery_conf )
time_m , time_s = divmod ( time_elapsed , 60 )
example_nb . save_file ( )
with open ( os . path . join ( target_dir , base_image_name + '.rst' ) , 'w' ) as f :
    example_rst += CODE_DOWNLOAD . format ( time_m , time_s , fname , example_nb . file_name )
    f . write ( example_rst )
return amount_of_code


signature:
def _arburg2 ( X , order ) :
    """This version is 10 times faster than arburg , but the output rho is not correct .
 returns [ 1 a0 , a1 , an - 1]"""

implementation:
x = np . array ( X )
N = len ( x )
if order <= 0. :
    raise ValueError ( "order must be > 0" )
# Initialisation
# - - - - - rho , den
rho = sum ( abs ( x ) ** 2. ) / N
# Eq 8.21 [ Marple ] _
den = rho * 2. * N
# - - - - - backward and forward errors
ef = np . zeros ( N , dtype = complex )
eb = np . zeros ( N , dtype = complex )
for j in range ( 0 , N ) : # eq 8.11
    ef [ j ] = x [ j ]
    eb [ j ] = x [ j ]
# AR order to be stored
a = np . zeros ( 1 , dtype = complex )
a [ 0 ] = 1
# - - - - rflection coeff to be stored
ref = np . zeros ( order , dtype = complex )
temp = 1.
E = np . zeros ( order + 1 )
E [ 0 ] = rho
for m in range ( 0 , order ) : # print m
# Calculate the next order reflection ( parcor ) coefficient
    efp = ef [ 1 : ]
    ebp = eb [ 0 : - 1 ]
    # print efp , ebp
    num = - 2. * np . dot ( ebp . conj ( ) . transpose ( ) , efp )
    den = np . dot ( efp . conj ( ) . transpose ( ) , efp )
    den += np . dot ( ebp , ebp . conj ( ) . transpose ( ) )
    ref [ m ] = num / den
    # Update the forward and backward prediction errors
    ef = efp + ref [ m ] * ebp
    eb = ebp + ref [ m ] . conj ( ) . transpose ( ) * efp
    # Update the AR coeff .
    a . resize ( len ( a ) + 1 )
    a = a + ref [ m ] * np . flipud ( a ) . conjugate ( )
    # Update the prediction error
    E [ m + 1 ] = ( 1 - ref [ m ] . conj ( ) . transpose ( ) * ref [ m ] ) * E [ m ]
    # print ' REF ' , ref , num , den
return a , E [ - 1 ] , ref


signature:
def arburg ( X , order , criteria = None ) :
    r"""Estimate the complex autoregressive parameters by the Burg algorithm .
 . . math : : x ( n ) = \ sqrt { ( v } ) e ( n ) + \ sum _ { k = 1 } ^ { P + 1 } a ( k ) x ( n - k )
 : param x : Array of complex data samples ( length N )
 : param order : Order of autoregressive process ( 0 < order < N )
 : param criteria : select a criteria to automatically select the order
 : return :
 * A Array of complex autoregressive parameters A ( 1 ) to A ( order ) . First
 value ( unity ) is not included ! !
 * P Real variable representing driving noise variance ( mean square
 of residual noise ) from the whitening operation of the Burg
 filter .
 * reflection coefficients defining the filter of the model .
 . . plot : :
 : width : 80%
 : include - source :
 from pylab import plot , log10 , linspace , axis
 from spectrum import *
 AR , P , k = arburg ( marple _ data , 15)
 PSD = arma2psd ( AR , sides = ' centerdc ' )
 plot ( linspace ( - 0.5 , 0.5 , len ( PSD ) ) , 10 * log10 ( PSD / max ( PSD ) ) )
 axis ( [ - 0.5,0.5 , - 60,0 ] )
 . . note : :
 1 . no detrend . Should remove the mean trend to get PSD . Be careful if
 presence of large mean .
 2 . If you don ' t know what the order value should be , choose the
 criterion = ' AKICc ' , which has the least bias and best
 resolution of model - selection criteria .
 . . note : : real and complex results double - checked versus octave using
 complex 64 samples stored in marple _ data . It does not agree with Marple
 fortran routine but this is due to the simplex precision of complex
 data in fortran .
 : reference : [ Marple ] _ [ octave ] _"""

implementation:
if order <= 0. :
    raise ValueError ( "order must be > 0" )
if order > len ( X ) :
    raise ValueError ( "order must be less than length input - 2" )
x = np . array ( X )
N = len ( x )
# Initialisation
# - - - - - rho , den
rho = sum ( abs ( x ) ** 2. ) / float ( N )
# Eq 8.21 [ Marple ] _
den = rho * 2. * N
# - - - - criteria
if criteria :
    from spectrum import Criteria
    crit = Criteria ( name = criteria , N = N )
    crit . data = rho
    logging . debug ( 'Step {}. old criteria={} new one={}.  rho={}' . format ( 0 , crit . old_data , crit . data , rho ) )
# p = 0
a = np . zeros ( 0 , dtype = complex )
ref = np . zeros ( 0 , dtype = complex )
ef = x . astype ( complex )
eb = x . astype ( complex )
temp = 1.
# Main recursion
for k in range ( 0 , order ) : # calculate the next order reflection coefficient Eq 8.14 Marple
    num = sum ( [ ef [ j ] * eb [ j - 1 ] . conjugate ( ) for j in range ( k + 1 , N ) ] )
    den = temp * den - abs ( ef [ k ] ) ** 2 - abs ( eb [ N - 1 ] ) ** 2
    kp = - 2. * num / den
    # eq 8.14
    temp = 1. - abs ( kp ) ** 2.
    new_rho = temp * rho
    if criteria :
        logging . debug ( 'Step {}. old criteria={} new one={}. rho={}' . format ( k + 1 , crit . old_data , crit . data , new_rho ) )
        # k + 1 because order goes from 1 to P whereas k starts at 0.
        status = crit ( rho = temp * rho , k = k + 1 )
        if status is False :
            logging . debug ( 'Stop criteria reached %s %s ' % ( crit . data , crit . old_data ) )
            break
    # this should be after the criteria
    rho = new_rho
    if rho <= 0 :
        raise ValueError ( "Found a negative value (expected positive stricly) %s. Decrease the order" % rho )
    a . resize ( a . size + 1 )
    a [ k ] = kp
    if k == 0 :
        for j in range ( N - 1 , k , - 1 ) :
            save2 = ef [ j ]
            ef [ j ] = save2 + kp * eb [ j - 1 ]
            # Eq . ( 8.7)
            eb [ j ] = eb [ j - 1 ] + kp . conjugate ( ) * save2
    else : # update the AR coeff
        khalf = ( k + 1 ) // 2
        # FIXME here khalf must be an integer
        for j in range ( 0 , khalf ) :
            ap = a [ j ]
            # previous value
            a [ j ] = ap + kp * a [ k - j - 1 ] . conjugate ( )
            # Eq . ( 8.2)
            if j != k - j - 1 :
                a [ k - j - 1 ] = a [ k - j - 1 ] + kp * ap . conjugate ( )
        # Eq . ( 8.2)
        # update the prediction error
        for j in range ( N - 1 , k , - 1 ) :
            save2 = ef [ j ]
            ef [ j ] = save2 + kp * eb [ j - 1 ]
            # Eq . ( 8.7)
            eb [ j ] = eb [ j - 1 ] + kp . conjugate ( ) * save2
    # save the reflection coefficient
    ref . resize ( ref . size + 1 )
    ref [ k ] = kp
return a , rho , ref


signature:
def _numpy_cholesky ( A , B ) :
    """Solve Ax = B using numpy cholesky solver
 A = LU
 in the case where A is square and Hermitian , A = L . L * where L * is
 transpoed and conjugate matrix
 Ly = b
 where
 Ux = y
 so x = U ^ { - 1 } y
 where U = L *
 and y = L ^ { - 1 } B"""

implementation:
L = numpy . linalg . cholesky ( A )
# A = L * numpy . transpose ( L ) . conjugate ( )
# Ly = b
y = numpy . linalg . solve ( L , B )
# Ux = y
x = numpy . linalg . solve ( L . transpose ( ) . conjugate ( ) , y )
return x , L


signature:
def _numpy_solver ( A , B ) :
    """This function solve Ax = B directly without taking care of the input
 matrix properties ."""

implementation:
x = numpy . linalg . solve ( A , B )
return x


signature:
def CHOLESKY ( A , B , method = 'scipy' ) :
    """Solve linear system ` AX = B ` using CHOLESKY method .
 : param A : an input Hermitian matrix
 : param B : an array
 : param str method : a choice of method in [ numpy , scipy , numpy _ solver ]
 * ` numpy _ solver ` relies entirely on numpy . solver ( no cholesky decomposition )
 * ` numpy ` relies on the numpy . linalg . cholesky for the decomposition and
 numpy . linalg . solve for the inversion .
 * ` scipy ` uses scipy . linalg . cholesky for the decomposition and
 scipy . linalg . cho _ solve for the inversion .
 . . rubric : : Description
 When a matrix is square and Hermitian ( symmetric with lower part being
 the complex conjugate of the upper one ) , then the usual triangular
 factorization takes on the special form :
 . . math : : A = R R ^ H
 where : math : ` R ` is a lower triangular matrix with nonzero real principal
 diagonal element . The input matrix can be made of complex data . Then , the
 inversion to find : math : ` x ` is made as follows :
 . . math : : Ry = B
 and
 . . math : : Rx = y
 . . doctest : :
 > > > import numpy
 > > > from spectrum import CHOLESKY
 > > > A = numpy . array ( [ [ 2.0 + 0 . j , 0.5-0.5j , - 0.2 + 0.1j ] ,
 . . . [ 0.5 + 0.5j , 1.0 + 0 . j , 0.3-0.2j ] ,
 . . . [ - 0.2-0.1j , 0.3 + 0.2j , 0.5 + 0 . j ] ] )
 > > > B = numpy . array ( [ 1.0 + 3 . j , 2.0-1 . j , 0.5 + 0.8j ] )
 > > > CHOLESKY ( A , B )
 array ( [ 0.95945946 + 5.25675676j , 4.41891892-7.04054054j ,
 -5.13513514 + 6.35135135j ] )"""

implementation:
if method == 'numpy_solver' :
    X = _numpy_solver ( A , B )
    return X
elif method == 'numpy' :
    X , _L = _numpy_cholesky ( A , B )
    return X
elif method == 'scipy' :
    import scipy . linalg
    L = scipy . linalg . cholesky ( A )
    X = scipy . linalg . cho_solve ( ( L , False ) , B )
else :
    raise ValueError ( 'method must be numpy_solver, numpy_cholesky or cholesky_inplace' )
return X


signature:
def music ( X , IP , NSIG = None , NFFT = default_NFFT , threshold = None , criteria = 'aic' , verbose = False ) :
    """Eigen value pseudo spectrum estimate . See : func : ` eigenfre `"""

implementation: return eigen ( X , IP , NSIG = NSIG , method = 'music' , NFFT = NFFT , threshold = threshold , criteria = criteria , verbose = verbose )


signature:
def eigen ( X , P , NSIG = None , method = 'music' , threshold = None , NFFT = default_NFFT , criteria = 'aic' , verbose = False ) :
    r"""Pseudo spectrum using eigenvector method ( EV or Music )
 This function computes either the Music or EigenValue ( EV ) noise
 subspace frequency estimator .
 First , an autocorrelation matrix of order ` P ` is computed from
 the data . Second , this matrix is separated into vector subspaces ,
 one a signal subspace and the other a noise
 subspace using a SVD method to obtain the eigen values and vectors .
 From the eigen values : math : ` \ lambda _ i ` , and eigen vectors : math : ` v _ k ` ,
 the * * pseudo spectrum * * ( see note below ) is computed as follows :
 . . math : : P _ { ev } ( f ) = \ frac { 1 } { e ^ H ( f ) \ left ( \ sum \ limits _ { k = M + 1 } ^ { p } \ frac { 1 } { \ lambda _ k } v _ kv _ k ^ H \ right ) e ( f ) }
 The separation of the noise and signal subspaces requires expertise
 of the signal . However , AIC and MDL criteria may be used to automatically
 perform this task .
 You still need to provide the parameter ` P ` to indicate the maximum number
 of eigen values to be computed . The criteria will just select a subset
 to estimate the pseudo spectrum ( see : func : ` ~ spectrum . criteria . aic _ eigen `
 and : func : ` ~ spectrum . criteria . mdl _ eigen ` for details .
 . . note : : * * pseudo spectrum * * . func : ` eigen ` does not compute a PSD estimate .
 Indeed , the method does not preserve the measured process power .
 : param X : Array data samples
 : param int P : maximum number of eigen values to compute . NSIG ( if
 specified ) must therefore be less than P .
 : param str method : ' music ' or ' ev ' .
 : param int NSIG : If specified , the signal sub space uses NSIG eigen values .
 : param float threshold : If specified , the signal sub space is made of the
 eigen values larger than : math : ` \ rm { threshold } \ times \ lambda _ { min } ` ,
 where : math : ` \ lambda _ { min } ` is the minimum eigen values .
 : param int NFFT : total length of the final data sets ( padded with zero
 if needed ; default is 4096)
 : return :
 * PSD : Array of real frequency estimator values ( two sided for
 complex data and one sided for real data )
 * S , the eigen values
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import eigen , marple _ data
 from pylab import plot , log10 , linspace , legend , axis
 psd , ev = eigen ( marple _ data , 15 , NSIG = 11)
 f = linspace ( - 0.5 , 0.5 , len ( psd ) )
 plot ( f , 10 * log10 ( psd / max ( psd ) ) , label = ' User defined ' )
 psd , ev = eigen ( marple _ data , 15 , threshold = 2)
 plot ( f , 10 * log10 ( psd / max ( psd ) ) , label = ' threshold method ( 100 ) ' )
 psd , ev = eigen ( marple _ data , 15)
 plot ( f , 10 * log10 ( psd / max ( psd ) ) , label = ' AIC method ( 8 ) ' )
 legend ( )
 axis ( [ - 0.5 , 0.5 , - 120 , 0 ] )
 . . seealso : :
 : func : ` pev ` ,
 : func : ` pmusic ` ,
 : func : ` ~ spectrum . criteria . aic _ eigen `
 : References : [ Marple ] _ , Chap 13
 . . todo : : for developers :
 * what should be the second argument of the criteria N , N - P , P . . . ?
 * what should be the max value of NP"""

implementation:
if method not in [ 'music' , 'ev' ] :
    raise ValueError ( "method must be 'music' or 'ev'" )
if NSIG != None and threshold != None :
    raise ValueError ( "NSIG and threshold cannot be provided together" )
if NSIG is not None :
    if NSIG < 0 :
        raise ValueError ( 'NSIG must be positive' )
    if NSIG >= P :
        raise ValueError ( "NSIG must be stricly less than IP" )
N = len ( X )
NP = N - P
assert 2 * NP > P - 1 , 'decrease the second argument'
if NP > 100 :
    NP = 100
FB = np . zeros ( ( 2 * NP , P ) , dtype = complex )
# FB = numpy . zeros ( ( MAXU , IP ) , dtype = complex )
Z = np . zeros ( NFFT , dtype = complex )
PSD = np . zeros ( NFFT )
# These loops can surely be replaced by a function that create such matrix
for I in range ( 0 , NP ) :
    for K in range ( 0 , P ) :
        FB [ I , K ] = X [ I - K + P - 1 ]
        FB [ I + NP , K ] = X [ I + K + 1 ] . conjugate ( )
# This commented line produces the correct FB , as the 2 for loops above
# It is more elegant but slower . . . corrmtx needs to be optimised ( 20/4/11)
# FB2 = spectrum . linalg . corrmtx ( X , P - 1 , method = ' modified ' )
# Compute the eigen values / vectors
_U , S , V = svd ( FB )
# U and V are not the same as in Marple . Real or Imaginary absolute values
# are correct but signs are not . This is wierd because the svd function
# gives the same result as cvsd in Marple . Is FB correct ? it seems so .
# The following operation has to be done . Otherwise , the resulting PSD is
# not corect
V = - V . transpose ( )
NSIG = _get_signal_space ( S , 2 * NP , verbose = verbose , threshold = threshold , NSIG = NSIG , criteria = criteria )
# C AI or Expert Knowledge to choose " signal " singular values , or input
# C NSIG at this point
for I in range ( NSIG , P ) :
    Z [ 0 : P ] = V [ 0 : P , I ]
    Z [ P : NFFT ] = 0
    Z = fft ( Z , NFFT )
    if method == 'music' :
        PSD = PSD + abs ( Z ) ** 2.
    elif method == 'ev' :
        PSD = PSD + abs ( Z ) ** 2. / S [ I ]
PSD = 1. / PSD
# for some reasons , we need to rearrange the output . this is related to
# the way U and V are order in the routine svd
nby2 = int ( NFFT / 2 )
# return PSD , S
newpsd = np . append ( PSD [ nby2 : 0 : - 1 ] , PSD [ nby2 * 2 - 1 : nby2 - 1 : - 1 ] )
return newpsd , S


signature:
def _get_signal_space ( S , NP , verbose = False , threshold = None , NSIG = None , criteria = 'aic' ) :
    """todo"""

implementation:
from . criteria import aic_eigen , mdl_eigen
# This section selects automatically the noise and signal subspaces .
# NSIG being the number of eigenvalues corresponding to signals .
if NSIG is None :
    if threshold is None :
        logging . debug ( 'computing NSIG using AIC method' )
        # get the minimum index of the AIC vector
        if criteria == 'aic' :
            aic = aic_eigen ( S , NP * 2 )
        elif criteria == 'mdl' :
            aic = mdl_eigen ( S , NP * 2 )
        # get the minimum index of the AIC vector , add 1 to get the NSIG
        NSIG = np . argmin ( aic ) + 1
        logging . debug ( 'NSIG=' , NSIG , ' found as the number of pertinent sinusoids' )
    else :
        logging . debug ( 'computing NSIG using user threshold ' )
        # following an idea from Matlab , pmusic , we look at the minimum
        # eigen value , and split the eigen values above and below
        # K times min eigen value , where K is > 1
        m = threshold * min ( S )
        new_s = S [ np . where ( S > m ) ]
        NSIG = len ( new_s )
        logging . debug ( 'found' , NSIG )
        if NSIG == 0 :
            NSIG = 1
return NSIG


signature:
def speriodogram ( x , NFFT = None , detrend = True , sampling = 1. , scale_by_freq = True , window = 'hamming' , axis = 0 ) :
    """Simple periodogram , but matrices accepted .
 : param x : an array or matrix of data samples .
 : param NFFT : length of the data before FFT is computed ( zero padding )
 : param bool detrend : detrend the data before co , puteing the FFT
 : param float sampling : sampling frequency of the input : attr : ` data ` .
 : param scale _ by _ freq :
 : param str window :
 : return : 2 - sided PSD if complex data , 1 - sided if real .
 if a matrix is provided ( using numpy . matrix ) , then a periodogram
 is computed for each row . The returned matrix has the same shape as the input
 matrix .
 The mean of the input data is also removed from the data before computing
 the psd .
 . . plot : :
 : width : 80%
 : include - source :
 from pylab import grid , semilogy
 from spectrum import data _ cosine , speriodogram
 data = data _ cosine ( N = 1024 , A = 0.1 , sampling = 1024 , freq = 200)
 semilogy ( speriodogram ( data , detrend = False , sampling = 1024 ) , marker = ' o ' )
 grid ( True )
 . . plot : :
 : width : 80%
 : include - source :
 import numpy
 from spectrum import speriodogram , data _ cosine
 from pylab import figure , semilogy , figure , imshow
 # create N data sets and make the frequency dependent on the time
 N = 100
 m = numpy . concatenate ( [ data _ cosine ( N = 1024 , A = 0.1 , sampling = 1024 , freq = x )
 for x in range ( 1 , N ) ] ) ;
 m . resize ( N , 1024)
 res = speriodogram ( m )
 figure ( 1)
 semilogy ( res )
 figure ( 2)
 imshow ( res . transpose ( ) , aspect = ' auto ' )
 . . todo : : a proper spectrogram class / function that takes care of normalisation"""

implementation:
x = np . array ( x )
# array with 1 dimension case
if x . ndim == 1 :
    axis = 0
    r = x . shape [ 0 ]
    w = Window ( r , window )
    # same size as input data
    w = w . data
# matrix case
elif x . ndim == 2 :
    logging . debug ( '2D array. each row is a 1D array' )
    [ r , c ] = x . shape
    w = np . array ( [ Window ( r , window ) . data for this in range ( c ) ] ) . reshape ( r , c )
if NFFT is None :
    NFFT = len ( x )
isreal = np . isrealobj ( x )
if detrend == True :
    m = np . mean ( x , axis = axis )
else :
    m = 0
if isreal == True :
    if x . ndim == 2 :
        res = ( abs ( rfft ( x * w - m , NFFT , axis = 0 ) ) ) ** 2. / r
    else :
        res = ( abs ( rfft ( x * w - m , NFFT , axis = - 1 ) ) ) ** 2. / r
else :
    if x . ndim == 2 :
        res = ( abs ( fft ( x * w - m , NFFT , axis = 0 ) ) ) ** 2. / r
    else :
        res = ( abs ( fft ( x * w - m , NFFT , axis = - 1 ) ) ) ** 2. / r
if scale_by_freq is True :
    df = sampling / float ( NFFT )
    res *= 2 * np . pi / df
if x . ndim == 1 :
    return res . transpose ( )
else :
    return res


signature:
def WelchPeriodogram ( data , NFFT = None , sampling = 1. , ** kargs ) :
    r"""Simple periodogram wrapper of numpy . psd function .
 : param A : the input data
 : param int NFFT : total length of the final data sets ( padded
 with zero if needed ; default is 4096)
 : param str window :
 : Technical documentation :
 When we calculate the periodogram of a set of data we get an estimation
 of the spectral density . In fact as we use a Fourier transform and a
 truncated segments the spectrum is the convolution of the data with a
 rectangular window which Fourier transform is
 . . math : :
 W ( s ) = \ frac { 1 } { N ^ 2 } \ left [ \ frac { \ sin ( \ pi s ) } { \ sin ( \ pi s / N ) } \ right ] ^ 2
 Thus oscillations and sidelobes appears around the main frequency . One aim of t he tapering is to reduced this effects . We multiply data by a window whose sidelobes are much smaller than the main lobe . Classical window is hanning window . But other windows are available . However we must take into account this energy and divide the spectrum by energy of taper used . Thus periodogram becomes :
 . . math : :
 D _ k \ equiv \ sum _ { j = 0 } ^ { N - 1 } c _ jw _ j \ ; e ^ { 2 \ pi ijk / N } \ qquad k = 0 , . . . , N - 1
 . . math : :
 P ( 0 ) = P ( f _ 0 ) = \ frac { 1 } { 2 \ pi W _ { ss } } \ arrowvert { D _ 0 } \ arrowvert ^ 2
 . . math : :
 P ( f _ k ) = \ frac { 1 } { 2 \ pi W _ { ss } } \ left [ \ arrowvert { D _ k } \ arrowvert ^ 2 + \ arrowvert { D _ { N - k } } \ arrowvert ^ 2 \ right ] \ qquad k = 0,1 , . . . , \ left ( \ frac { 1 } { 2 } - 1 \ right )
 . . math : :
 P ( f _ c ) = P ( f _ { N / 2 } ) = \ frac { 1 } { 2 \ pi W _ { ss } } \ arrowvert { D _ { N / 2 } } \ arrowvert ^ 2
 with
 . . math : :
 { W _ { ss } } \ equiv N \ sum _ { j = 0 } ^ { N - 1 } w _ j ^ 2
 . . plot : :
 : width : 80%
 : include - source :
 from spectrum import WelchPeriodogram , marple _ data
 psd = WelchPeriodogram ( marple _ data , 256)"""

implementation:
from pylab import psd
spectrum = Spectrum ( data , sampling = 1. )
P = psd ( data , NFFT , Fs = sampling , ** kargs )
spectrum . psd = P [ 0 ]
# spectrum . _ _ Spectrum _ sides = ' twosided '
return P , spectrum


signature:
def DaniellPeriodogram ( data , P , NFFT = None , detrend = 'mean' , sampling = 1. , scale_by_freq = True , window = 'hamming' ) :
    r"""Return Daniell ' s periodogram .
 To reduce fast fluctuations of the spectrum one idea proposed by daniell
 is to average each value with points in its neighboorhood . It ' s like
 a low filter .
 . . math : : \ hat { P } _ D [ f _ i ] = \ frac { 1 } { 2P + 1 } \ sum _ { n = i - P } ^ { i + P } \ tilde { P } _ { xx } [ f _ n ]
 where P is the number of points to average .
 Daniell ' s periodogram is the convolution of the spectrum with a low filter :
 . . math : : \ hat { P } _ D ( f ) = \ hat { P } _ { xx } ( f ) * H ( f )
 Example : :
 > > > DaniellPeriodogram ( data , 8)
 if N / P is not integer , the final values of the original PSD are not used .
 using DaniellPeriodogram ( data , 0 ) should give the original PSD ."""

implementation:
psd = speriodogram ( data , NFFT = NFFT , detrend = detrend , sampling = sampling , scale_by_freq = scale_by_freq , window = window )
if len ( psd ) % 2 == 1 :
    datatype = 'real'
else :
    datatype = 'complex'
N = len ( psd )
_slice = 2 * P + 1
if datatype == 'real' : # must get odd value
    newN = np . ceil ( psd . size / float ( _slice ) )
    if newN % 2 == 0 :
        newN = psd . size / _slice
else :
    newN = np . ceil ( psd . size / float ( _slice ) )
    if newN % 2 == 1 :
        newN = psd . size / _slice
newpsd = np . zeros ( int ( newN ) )
# keep integer division
for i in range ( 0 , newpsd . size ) :
    count = 0
    # needed to know the number of valid averaged values
    for n in range ( i * _slice - P , i * _slice + P + 1 ) : # + 1 to have P values on each sides
        if n > 0 and n < N : # needed to start the average
            count += 1
            newpsd [ i ] += psd [ n ]
    newpsd [ i ] /= float ( count )
# todo : check this
if datatype == 'complex' :
    freq = np . linspace ( 0 , sampling , len ( newpsd ) )
else :
    df = 1. / sampling
    freq = np . linspace ( 0 , sampling / 2. , len ( newpsd ) )
# psd . refreq ( 2 * psd . size ( ) / A . freq ( ) ) ;
# psd . retime ( - 1 . / psd . freq ( ) + 1 . / A . size ( ) ) ;
return newpsd , freq


signature:
def centerdc_gen ( self ) :
    """Return the centered frequency range as a generator .
 > > > print ( list ( Range ( 8 ) . centerdc _ gen ( ) ) )
 [ - 0.5 , - 0.375 , - 0.25 , - 0.125 , 0.0 , 0.125 , 0.25 , 0.375]"""

implementation:
for a in range ( 0 , self . N ) :
    yield ( a - self . N / 2 ) * self . df


signature:
def onesided_gen ( self ) :
    """Return the one - sided frequency range as a generator .
 If : attr : ` N ` is even , the length is N / 2 + 1.
 If : attr : ` N ` is odd , the length is ( N + 1 ) / 2.
 > > > print ( list ( Range ( 8 ) . onesided ( ) ) )
 [0.0 , 0.125 , 0.25 , 0.375 , 0.5]
 > > > print ( list ( Range ( 9 ) . onesided ( ) ) )
 [0.0 , 0.1111 , 0.2222 , 0.3333 , 0.4444]"""

implementation:
if self . N % 2 == 0 :
    for n in range ( 0 , self . N // 2 + 1 ) :
        yield n * self . df
else :
    for n in range ( 0 , ( self . N + 1 ) // 2 ) :
        yield n * self . df


signature:
def frequencies ( self , sides = None ) :
    """Return the frequency vector according to : attr : ` sides `"""

implementation:
# use the attribute sides except if a valid sides argument is provided
if sides is None :
    sides = self . sides
if sides not in self . _sides_choices :
    raise errors . SpectrumChoiceError ( sides , self . _sides_choices )
if sides == 'onesided' :
    return self . _range . onesided ( )
if sides == 'twosided' :
    return self . _range . twosided ( )
if sides == 'centerdc' :
    return self . _range . centerdc ( )