Line data Source code
1 : /** \file linearPredictor.hpp
2 : * \brief Working with linear prediction.
3 : *
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : * \ingroup signal_processing_files
7 : *
8 : */
9 :
10 : //***********************************************************************//
11 : // Copyright 2015, 2016, 2017 Jared R. Males (jaredmales@gmail.com)
12 : //
13 : // This file is part of mxlib.
14 : //
15 : // mxlib is free software: you can redistribute it and/or modify
16 : // it under the terms of the GNU General Public License as published by
17 : // the Free Software Foundation, either version 3 of the License, or
18 : // (at your option) any later version.
19 : //
20 : // mxlib is distributed in the hope that it will be useful,
21 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
22 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 : // GNU General Public License for more details.
24 : //
25 : // You should have received a copy of the GNU General Public License
26 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
27 : //***********************************************************************//
28 :
29 : #ifndef linearPredictor_hpp
30 : #define linearPredictor_hpp
31 :
32 : #include <vector>
33 : #include <complex>
34 :
35 : #include "../mxlib.hpp"
36 : #include "../math/constants.hpp"
37 : #include "../math/eigenLapack.hpp"
38 :
39 : #include "levinsonRecursion.hpp"
40 :
41 : namespace mx
42 : {
43 : namespace sigproc
44 : {
45 :
46 : #define LP_BREADCRUMB
47 :
48 : //#define LP_BREADCRUMB std::cerr << __FILE__ << ' ' << __LINE__ << '\n';
49 :
50 :
51 : /// A class to support linear prediction.
52 : /** \ingroup signal_processing
53 : *
54 : * \todo document linearPredictor
55 : */
56 : template <typename _realT>
57 : struct linearPredictor
58 : {
59 : typedef _realT realT;
60 :
61 : std::vector<realT> m_c;
62 :
63 : realT _setCondition{ 0 };
64 : realT _actCondition{ 0 };
65 : int _nRejected{ 0 };
66 :
67 : /// Calculate the LP coefficients given an autocorrelation.
68 : /** If condition==0 then the levinson recursion is used.
69 : * Otherwise, SVD pseudo-inversion is used with the given condition number.
70 : */
71 121 : int calcCoefficients( const std::vector<realT> &ac,
72 : size_t Nc,
73 : size_t Npred = 1,
74 : realT condition = 0
75 : )
76 : {
77 121 : return calcCoefficients( ac.data(), ac.size(), Nc, Npred, condition);
78 : }
79 :
80 121 : int calcCoefficients( const realT *ac,
81 : size_t acSz,
82 : size_t Nc,
83 : size_t Npred = 1,
84 : realT condition = 0
85 : )
86 : {
87 :
88 121 : if( condition == 0 )
89 : {
90 : LP_BREADCRUMB;
91 121 : return calcCoefficientsLevinson( ac, acSz, Nc, Npred );
92 : }
93 :
94 0 : Eigen::Array<realT, -1, -1> Rmat, Rvec, PInv, LPcoeff;
95 :
96 0 : Rmat.resize( Nc, Nc );
97 0 : Rvec.resize( 1, Nc );
98 :
99 0 : for( int i = 0; i < Nc && i < acSz; ++i )
100 : {
101 0 : for( int j = 0; j < Nc && i < acSz; ++j )
102 : {
103 0 : Rmat( i, j ) = ac[abs( i - j )];
104 : }
105 :
106 0 : Rvec( 0, i ) = ac[i + Npred];
107 : }
108 :
109 0 : realT tmpCond = condition;
110 :
111 0 : _setCondition = condition;
112 0 : math::eigenPseudoInverse( PInv, tmpCond, _nRejected, Rmat, condition );
113 :
114 0 : _actCondition = tmpCond;
115 :
116 0 : m_c.resize(Nc);
117 0 : Eigen::Map<Eigen::Array<realT,-1,-1>> cmap(m_c.data(), 1, m_c.size());
118 0 : cmap = Rvec.matrix() * PInv.matrix();
119 :
120 0 : return 0;
121 0 : }
122 :
123 0 : int calcCoefficientsLevinson( const std::vector<realT> &ac, /**< [in] The autocorrelation, at least
124 : Nc+Npred in length */
125 : size_t Nc, /**< [in] The number of LP coefficients desired */
126 : size_t Npred = 1 /**< [in] [optional] The prediction length,
127 : default is 1 */
128 : )
129 : {
130 : LP_BREADCRUMB;
131 :
132 0 : return calcCoefficientsLevinson( ac.data(), ac.size(), Nc, Npred );
133 : }
134 :
135 121 : int calcCoefficientsLevinson( const realT *ac, /**< [in] The autocorrelation, at least Nc+Npred in length */
136 : size_t acSz, /**< [in] The length of the autocorrelation */
137 : size_t Nc, /**< [in] The number of LP coefficients desired */
138 : unsigned Npred = 1 /**< [in] [optional] The prediction length, default is 1 */
139 : )
140 : {
141 : LP_BREADCRUMB;
142 121 : if( acSz < Nc + Npred )
143 : {
144 0 : std::string msg = "too many coefficients for size and prediction length\n";
145 0 : msg += " acSz = " + std::to_string( acSz ) + "\n";
146 0 : msg += " Nc = " + std::to_string( Nc ) + "\n";
147 0 : msg += " Npred = " + std::to_string( Npred ) + "\n";
148 0 : throw(mx::exception( error_t::invalidarg ,msg));
149 0 : }
150 :
151 121 : if(Nc == 0)
152 : {
153 0 : std::string msg = "Nc can't be 0";
154 0 : throw(mx::exception( error_t::invalidarg , msg));
155 0 : }
156 :
157 121 : std::vector<realT> r, x, y;
158 :
159 : LP_BREADCRUMB;
160 121 : r.resize( 2. * Nc - 1 );
161 121 : m_c.resize( Nc );
162 121 : y.resize( Nc );
163 :
164 : LP_BREADCRUMB;
165 605 : for( size_t i = 0; i < Nc; ++i )
166 : {
167 484 : r[i] = ac[Nc - i - 1]; // this runs from Nc-1 to 0
168 : }
169 :
170 : LP_BREADCRUMB;
171 484 : for( size_t i = Nc; i < 2 * Nc - 1; ++i )
172 : {
173 363 : r[i] = ac[i - Nc + 1]; // this runs from 1 to Nc-1
174 : }
175 :
176 : LP_BREADCRUMB;
177 605 : for( size_t i = 0; i < Nc; ++i )
178 : {
179 484 : y[i] = ac[i + Npred]; // this runs from Npred to Nc-1 + Npred
180 : }
181 :
182 : LP_BREADCRUMB;
183 242 : return levinsonRecursion( r.data(), m_c.data(), y.data(), Nc );
184 121 : }
185 :
186 0 : realT c( size_t i )
187 : {
188 0 : return m_c[i];
189 : }
190 :
191 0 : size_t Nc()
192 : {
193 0 : return m_c.size();
194 : }
195 :
196 0 : realT predict( std::vector<realT> &hist, int idx )
197 : {
198 0 : realT x = 0;
199 :
200 0 : if( idx < m_c.size() )
201 : {
202 0 : return x;
203 : }
204 :
205 0 : for( int i = 0; i < m_c.size(); ++i )
206 : {
207 0 : x += m_c[i] * hist[idx - i];
208 : }
209 :
210 0 : return x;
211 : }
212 :
213 0 : realT spectralResponse( realT f, realT fs )
214 : {
215 0 : int n = m_c.size();
216 :
217 0 : std::complex<realT> He = 0;
218 0 : for( int j = 0; j < n; ++j )
219 : {
220 0 : realT s = ( j + 1.0 ) * math::two_pi<realT>();
221 0 : He += m_c[j] * exp( s * std::complex<realT>( 0, -1.0 ) * f / fs );
222 : }
223 :
224 0 : realT one = 1.0;
225 0 : return std::norm( one / ( one - He ) );
226 : }
227 : };
228 :
229 : } // namespace sigproc
230 : } // namespace mx
231 : #endif // linearPredictor_hpp
|
