Line data Source code
1 : /** \file psdVarMean.hpp
2 : * \brief Tools for calculating the variance of the mean of a PSD
3 : *
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : * \ingroup signal_processing_files
7 : *
8 : */
9 :
10 : //***********************************************************************//
11 : // Copyright 2018,2019,2020 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 psdVarMean_hpp
30 : #define psdVarMean_hpp
31 :
32 : #include <vector>
33 :
34 : #include <gsl/gsl_integration.h>
35 : #include <gsl/gsl_errno.h>
36 :
37 : #include "../math/func/bessel.hpp"
38 : #include "../math/constants.hpp"
39 : #include "../math/gslInterpolator.hpp"
40 :
41 : #include "psdFilter.hpp"
42 :
43 : namespace mx
44 : {
45 : namespace sigproc
46 : {
47 :
48 : /// Parameters for calculating the variance of the mean from a numerical PSD.
49 : /** Interpolates on the PSD for the integration.
50 : *
51 : * \tparam realT the real floating point type of the frequency and PSD vectors. [Must be double due to GSL]
52 : */
53 : template <typename _realT>
54 : struct psdVarMeanParams
55 : {
56 : typedef _realT realT;
57 :
58 : std::vector<realT> *freq{ nullptr }; ///< Pointer to the frequency vector
59 : std::vector<realT> *psd{ nullptr }; ///< Pointer to the psd vector
60 : realT T{ 0 }; ///< The sample length (units appropriate for freq)
61 :
62 : realT minf{ 0 }; ///< The minimum frequency
63 : realT maxf{ 0 }; ///< The maximum frequency
64 :
65 : realT extrapAlpha{ 0 };
66 :
67 : bool initialized{ false }; ///< Flag controlling whether the interpolator is initialized
68 :
69 : #ifdef MX_OLD_GSL
70 : math::gslInterpolator<math::gsl_interp_linear<realT>> terp; ///< Interpolator for the PSD
71 : #else
72 : math::gslInterpolator<math::gsl_interp_steffen<realT>> terp; ///< Interpolator for the PSD
73 : #endif
74 :
75 : /// Get the value of the PSD at k using interpolation
76 : /** Extrapolates outside the range of frequencies in freq.
77 : * Right now this is just constant.
78 : *
79 : * \todo implement better extrapolation past max-freq.
80 : */
81 0 : realT psdVal( realT k /**< [in] The transformed frequency coordinate */ )
82 : {
83 0 : if( !initialized )
84 : {
85 0 : if( freq == nullptr || psd == nullptr )
86 : {
87 0 : std::cerr << "Vectors not set.\n";
88 0 : exit( -1 );
89 : }
90 0 : if( T == 0 )
91 : {
92 0 : std::cerr << "T not set.\n";
93 0 : exit( -1 );
94 : }
95 :
96 0 : terp.setup( *freq, *psd );
97 : // #ifdef MX_OLD_GSL
98 : // terp.setup(gsl_interp_linear, *freq, *psd);
99 : // #else
100 : // terp.setup(gsl_interp_steffen, *freq, *psd);
101 : // #endif
102 :
103 0 : minf = ( *freq )[0];
104 0 : maxf = ( *freq )[freq->size() - 1];
105 :
106 0 : initialized = true;
107 : }
108 :
109 0 : realT f = ( 2.0 / T ) * k;
110 :
111 0 : if( f < minf )
112 0 : return ( *psd )[0];
113 :
114 0 : if( f > maxf )
115 : {
116 : // Enter extrapolation here
117 : // return 0;//(*psd)[psd->size()-1];
118 0 : if( extrapAlpha == 0 )
119 0 : return 0;
120 :
121 0 : return ( *psd )[psd->size() - 1] * pow( f / maxf, -1 * extrapAlpha );
122 : }
123 :
124 0 : return terp( f );
125 : }
126 : };
127 :
128 : /// Integration worker function for psdVarMean::operator()
129 : /**
130 : * \tparam paramsT the type of the parameter structure. See psdVarMean.
131 : *
132 : * \returns the value of the integrand at k
133 : *
134 : *
135 : */
136 : template <typename paramsT>
137 0 : typename paramsT::realT psdVarMeanFunc( typename paramsT::realT k, ///< [in] the scaled frequency coordinate
138 : void *params ///< [in] pointer to the parameter structure
139 : )
140 : {
141 : typedef typename paramsT::realT realT;
142 :
143 0 : paramsT *pvar = (paramsT *)params;
144 :
145 0 : realT J = math::func::bessel_j<realT, realT>( 0.5, math::two_pi<realT>() * k );
146 :
147 0 : return J * J / k * pvar->psdVal( k );
148 : }
149 :
150 : /// Calculate the variance of the mean for a process given its PSD
151 : /** Functor which manages the GSL integration workspace memory.
152 : *
153 : * \todo document the paramsT interface
154 : *
155 : * \tparam paramsT is a parameter structure with the correct interface.
156 : *
157 : * \ingroup psds
158 : */
159 : template <typename paramsT>
160 : struct psdVarMean
161 : {
162 : typedef typename paramsT::realT realT;
163 :
164 : realT extrapAlpha{ 0 }; ///< The extrapolation exponent. No extrapolation if 0.
165 :
166 : protected:
167 : gsl_integration_workspace *w{ nullptr }; ///< The GSL integration workspace memory
168 :
169 : size_t wSize{ 10000 }; ///< The size of the GSL workspace
170 :
171 : /// Initializes the gsl library.
172 0 : void init()
173 : {
174 : ///\todo We'll get the occasional failure to reach tolerance error, just ignore them all for now.
175 0 : gsl_set_error_handler_off();
176 0 : }
177 :
178 : public:
179 : /// Default c'tor
180 0 : psdVarMean()
181 0 : {
182 0 : init();
183 0 : }
184 :
185 : /// Constructor which initializes workspace size
186 0 : psdVarMean( size_t wm /**< [in] Size of the GSL integration workspace*/ )
187 0 : {
188 0 : init();
189 0 : wSize = wm;
190 0 : }
191 :
192 : /// d'tor cleans up the workspace memory
193 0 : ~psdVarMean()
194 : {
195 0 : if( w != nullptr )
196 : {
197 0 : gsl_integration_workspace_free( w );
198 : }
199 0 : }
200 :
201 : /// Calculate the variance of the mean for a psd over a sample length
202 : /**
203 : * \returns the variance of the mean for the process governed by the PSD.
204 : */
205 0 : realT operator()( realT &error, ///< [out] The error estimate returned by gsl_integration_qagiu
206 : std::vector<realT> &freq, ///< [in] The frequency vector
207 : std::vector<realT> &psd, ///< [in] The psd vector
208 : realT T ///< [in] The sample length (units approprate for freq)
209 : )
210 : {
211 0 : if( w == nullptr )
212 : {
213 0 : w = gsl_integration_workspace_alloc( wSize );
214 : }
215 :
216 : // Setup the integration parameters structure
217 0 : paramsT pvar;
218 :
219 0 : pvar.freq = &freq;
220 0 : pvar.psd = &psd;
221 0 : pvar.T = T;
222 :
223 0 : pvar.extrapAlpha = extrapAlpha;
224 :
225 : // Setup the function
226 : gsl_function func;
227 :
228 0 : func.function = &psdVarMeanFunc<paramsT>;
229 0 : func.params = &pvar;
230 :
231 : realT result;
232 :
233 0 : int ec = gsl_integration_qagiu( &func, 0.0, 1e-14, 1e-8, wSize, w, &result, &error );
234 : static_cast<void>( ec );
235 :
236 0 : return result / ( 2 * T );
237 0 : }
238 : };
239 :
240 : } // namespace sigproc
241 : } // namespace mx
242 :
243 : #endif // psdVarMean_hpp
|
