# Source code for spectrum.linear_prediction

```"""Linear prediction tools

:References: [Kay]_
"""

from .levinson import LEVINSON, rlevinson
import numpy
from scipy.signal import deconvolve

__all__ = ['ac2poly', 'poly2ac', 'ac2rc', 'rc2poly', 'rc2ac', 'is2rc', 'rc2is',
'rc2lar', 'lar2rc', 'poly2rc', 'lsf2poly', 'poly2lsf']

"""
lpc Linear prediction filter coefficients
schurrc Compute reflection coefficients from autocorrelation sequence
"""

[docs]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])
>>> e
4.1895000000000007

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

[docs]def ac2rc(data):
"""Convert autocorrelation sequence to reflection coefficients

:param data: an autorrelation vector
:return: the reflection coefficient and data[0]

This is an alias to::

a, e, c = LEVINSON(data)
c, data[0]

"""
a, e, _c = LEVINSON(data)
return a, data[0]

[docs]def poly2ac(poly, efinal):
""" Convert prediction filter polynomial to autocorrelation sequence

:param array poly: the AR parameters
:param efinal: an estimate of the final error
:return: the autocorrelation  sequence in complex format.

.. doctest::

>>> from numpy import array
>>> from spectrum import poly2ac
>>> poly = [ 1. ,  0.38 , -0.05]
>>> efinal = 4.1895
>>> poly2ac(poly, efinal)
array([ 5.00+0.j, -2.00+0.j,  1.01-0.j])

"""
results = rlevinson(poly, efinal)
return results[0]

def ar2rc(ar):
""" Convert autoregressive parameters into reflection coefficients """
raise NotImplementedError

[docs]def poly2rc(a, efinal):
"""Convert prediction filter polynomial to reflection coefficients

:param a: AR parameters
:param efinal:
"""
results = rlevinson(a, efinal)
return results[2]

[docs]def rc2poly(kr, r0=None):
"""convert reflection coefficients to prediction filter polynomial

:param k: reflection coefficients

"""
# 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

[docs]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`.

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

[docs]def is2rc(inv_sin):
"""Convert inverse sine parameters to reflection coefficients.

:param inv_sin: inverse sine parameters
:return: reflection coefficients

.. seealso::  :func:`rc2is`, :func:`poly2rc`, :func:`ac2rc`, :func:`lar2rc`.

:Reference: J.R. Deller, J.G. Proakis, J.H.L. Hansen,
"Discrete-Time Processing of Speech Signals", Prentice Hall, Section 7.4.5.

"""
return numpy.sin(numpy.array(inv_sin)*numpy.pi/2);

[docs]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.

"""
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)

[docs]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.

"""
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))

[docs]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.

"""
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)

[docs]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.00000000e+00,   6.14837835e-01,   9.89884967e-01,
# 9.31594056e-05,   3.13713832e-03,  -8.12002261e-03 ])

.. seealso:: poly2lsf, rc2poly, ac2poly, rc2is
"""
#   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] # do not return last element

[docs]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
"""

#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

def schurrc():
raise NotImplementedError

```