Line data Source code
1 :
2 : //***********************************************************************//
3 : // Copyright 2015, 2016, 2017 Jared R. Males (jaredmales@gmail.com)
4 : //
5 : // This file is part of mxlib.
6 : //
7 : // mxlib is free software: you can redistribute it and/or modify
8 : // it under the terms of the GNU General Public License as published by
9 : // the Free Software Foundation, either version 3 of the License, or
10 : // (at your option) any later version.
11 : //
12 : // mxlib is distributed in the hope that it will be useful,
13 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 : // GNU General Public License for more details.
16 : //
17 : // You should have received a copy of the GNU General Public License
18 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
19 : //***********************************************************************//
20 :
21 : #ifndef levinsonRecursion_hpp
22 : #define levinsonRecursion_hpp
23 :
24 : /*Solves the Toeplitz system PN
25 : j=1 R(N+i−j)xj = yi (i = 0,...,N-1). The Toeplitz matrix need
26 : not be symmetric. y[0..n-1] and r[0..2*n-2] are input arrays; x[0..n-1] is the output array.*/
27 :
28 : template <typename floatT>
29 121 : int levinsonRecursion( floatT *_r, floatT *_x, floatT *_y, int n )
30 : {
31 :
32 : int j, k, m, m1, m2;
33 : floatT pp, pt1, pt2, qq, qt1, qt2, sd, sgd, sgn, shn, sxn;
34 :
35 : /*Have to adjust for the stupid NR convention of 1-counting*/
36 121 : floatT *r = _r - 1;
37 121 : floatT *x = _x - 1;
38 121 : floatT *y = _y - 1;
39 :
40 : floatT *g, *h, *_g, *_h;
41 :
42 121 : if( r[n] == 0.0 )
43 : {
44 : //std::cerr << "toeplz-1 singular principal minor" << "\n";
45 0 : return -1;
46 : }
47 :
48 : /*Have to adjust for the stupid NR convention of 1-counting*/
49 121 : _g = new floatT[n];
50 121 : g = _g - 1;
51 121 : _h = new floatT[n];
52 121 : h = _h - 1;
53 :
54 121 : x[1] = y[1] / r[n]; // Initialize for the recursion.
55 :
56 121 : if( n == 1 )
57 : {
58 0 : delete[] _g;
59 0 : delete[] _h;
60 0 : return 0;
61 : }
62 :
63 121 : g[1] = r[n - 1] / r[n];
64 :
65 121 : h[1] = r[n + 1] / r[n];
66 :
67 363 : for( m = 1; m <= n; ++m )
68 : { // Main loop over the recursion.
69 :
70 363 : m1 = m + 1;
71 363 : sxn = -y[m1]; // Compute numerator and denominator for x,
72 363 : sd = -r[n];
73 :
74 1089 : for( j = 1; j <= m; ++j )
75 : {
76 726 : sxn += r[n + m1 - j] * x[j];
77 726 : sd += r[n + m1 - j] * g[m - j + 1];
78 : }
79 :
80 363 : if( sd == 0.0 )
81 : {
82 : //std::cerr << "toeplz-2 singular principal minor" << "\n";
83 :
84 0 : delete[] _g;
85 0 : delete[] _h;
86 0 : return -3;
87 : }
88 :
89 363 : x[m1] = sxn / sd; // whence x.
90 :
91 1089 : for( j = 1; j <= m; ++j )
92 726 : x[j] -= x[m1] * g[m - j + 1];
93 :
94 363 : if( m1 == n )
95 : {
96 121 : delete[] _g;
97 121 : delete[] _h;
98 121 : return 0;
99 : }
100 :
101 242 : sgn = -r[n - m1]; // Compute numerator and denominator for G and H,
102 :
103 242 : shn = -r[n + m1];
104 :
105 242 : sgd = -r[n];
106 :
107 605 : for( j = 1; j <= m; ++j )
108 : {
109 363 : sgn += r[n + j - m1] * g[j];
110 363 : shn += r[n + m1 - j] * h[j];
111 363 : sgd += r[n + j - m1] * h[m - j + 1];
112 : }
113 :
114 242 : if( sd == 0.0 || sgd == 0.0 )
115 : {
116 : //std::cerr << "toeplz-3 singular principal minor" << "\n";
117 :
118 0 : delete[] _g;
119 0 : delete[] _h;
120 0 : return -5;
121 : }
122 :
123 242 : g[m1] = sgn / sgd; // whence G and H.
124 242 : h[m1] = shn / sd;
125 :
126 242 : k = m;
127 242 : m2 = ( m + 1 ) >> 1;
128 242 : pp = g[m1];
129 242 : qq = h[m1];
130 484 : for( j = 1; j <= m2; ++j )
131 : {
132 242 : pt1 = g[j];
133 242 : pt2 = g[k];
134 242 : qt1 = h[j];
135 242 : qt2 = h[k];
136 242 : g[j] = pt1 - pp * qt2;
137 242 : g[k] = pt2 - pp * qt1;
138 242 : h[j] = qt1 - qq * pt2;
139 242 : h[k--] = qt2 - qq * pt1;
140 : }
141 : } // Back for another recurrence.
142 :
143 : //std::cerr << "toeplz - should not arrive here!" << "\n";
144 :
145 0 : delete[] _g;
146 0 : delete[] _h;
147 0 : return -6;
148 : }
149 :
150 : #endif // levinsonRecursion_hpp
|
