1 //

2 // Copyright (c) 2003-2011, MIST Project, Nagoya University

3 // All rights reserved.

4 //

5 // Redistribution and use in source and binary forms, with or without modification,

6 // are permitted provided that the following conditions are met:

7 //

8 // 1. Redistributions of source code must retain the above copyright notice,

9 // this list of conditions and the following disclaimer.

10 //

11 // 2. Redistributions in binary form must reproduce the above copyright notice,

12 // this list of conditions and the following disclaimer in the documentation

13 // and/or other materials provided with the distribution.

14 //

15 // 3. Neither the name of the Nagoya University nor the names of its contributors

16 // may be used to endorse or promote products derived from this software

17 // without specific prior written permission.

18 //

19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR

20 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND

21 // FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR

22 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL

23 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,

24 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER

25 // IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF

26 // THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

27 

28 

29 #ifndef __INCLUDE_MIST_PARCOR_H__

30 #define __INCLUDE_MIST_PARCOR_H__

31 

32 #ifndef __INCLUDE_MIST_H__

33 #include "../mist.h"

34 #endif

35 

36 #include <cmath>

37 

38 _MIST_BEGIN

39 

40 

41 

50 template < class T1, class Allocator1, class T2, class Allocator2, class T3, class Allocator3 >

51 inline double parcor( const array1< T1, Allocator1 > &cor, array1< T2, Allocator2 > &alf, array1< T3, Allocator3 > &ref, int p )

52 {

53  // by Levinson-Durbin algorithm

54  double res = cor[ 0 ]; // residual

55 

56  alf.resize( p );

57  ref.resize( p );

58 

59  array1< double > alf_old;

60 

61  size_t length = p;

62 

63  for( size_t i = 0 ; i < length ; i++ )

64  {

65  double r = cor[ i + 1 ];

66 

67  for( size_t j = 0 ; j < i ; j++ )

68  {

69  r += alf_old[ j ] * cor[ i - j ];

70  }

71 

72  ref[ i ] = static_cast < T3 >( r / res );

73 

74  alf[ i ] = -ref[ i ];

75 

76  for( size_t j = 0 ; j < i ; j++ )

77  {

78  alf[ j ] = static_cast < T2 >( alf_old[ j ] - ref[ i ] * alf_old[ i - j - 1 ] );

79  }

80 

81  res *= ( 1.0 - ref[ i ] ) * ( 1.0 + ref[ i ] );

82 

83  alf_old = alf;

84 

85  }

86 

87  return( res );

88 }

89 

90 

91 _MIST_END

92 

93 

94 #endif

95