Line data Source code
1 : /** \file psdUtils.hpp
2 : * \brief Tools for working with PSDs
3 : *
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : * \ingroup signal_processing_files
7 : *
8 : */
9 :
10 : //***********************************************************************//
11 : // Copyright 2015-2021 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 psdUtils_hpp
30 : #define psdUtils_hpp
31 :
32 : #ifndef EIGEN_NO_CUDA
33 : #define EIGEN_NO_CUDA
34 : #endif
35 :
36 : #include <type_traits>
37 :
38 : #include <Eigen/Dense>
39 :
40 : #include <iostream>
41 :
42 : #include "../mxlib.hpp"
43 :
44 : #include "../math/ft/fftT.hpp"
45 : #include "../math/vectorUtils.hpp"
46 :
47 : namespace mx
48 : {
49 : namespace sigproc
50 : {
51 :
52 : /** \ingroup psds
53 : * @{
54 : */
55 :
56 : /// Calculate the variance of a 1-D, 1-sided PSD
57 : /** By default uses trapezoid rule integration. This can be changed to mid-point integration.
58 : *
59 : * \returns the variance of a PSD (the integral).
60 : *
61 : * \tparam realT the real floating point type
62 : *
63 : * \test Scenario: calculating variance from a 1D PSD \ref tests_sigproc_psdUtils_psdVar_1D "[test doc]"
64 : */
65 : template <typename realT>
66 3 : realT psdVar1sided( realT df, ///< [in] the frequency scale of the PSD
67 : const realT *PSD, ///< [in] the PSD to integrate.
68 : size_t sz, ///< [in] the size of the PSD vector
69 : realT half = 0.5 ///< [in] [optional] controls if trapezoid (0.5) or mid-point (1.0) integration is
70 : ///< used. Do not use other values.
71 : )
72 : {
73 3 : realT var = 0;
74 :
75 3 : var = half * PSD[0];
76 :
77 12 : for( size_t i = 1; i < sz - 1; ++i )
78 : {
79 9 : var += PSD[i];
80 : }
81 :
82 3 : var += half * PSD[sz - 1];
83 :
84 3 : var *= df;
85 :
86 3 : return var;
87 : }
88 :
89 : /// Calculate the variance of a 1-D, 2-sided PSD
90 : /** By default uses trapezoid rule integration. This can be changed to mid-point integration.
91 : *
92 : * Assumes the 2-sided PSD is in standard FFT storage order, and that sz is even.
93 : *
94 : * \returns the variance of a PSD (the integral).
95 : *
96 : * \tparam realT the real floating point type
97 : *
98 : * \test Scenario: calculating variance from a 1D PSD. \ref tests_sigproc_psdUtils_psdVar_1D "[test doc]"
99 : */
100 : template <typename realT>
101 3 : realT psdVar2sided( realT df, ///< [in] the frequency scale of the PSD
102 : const realT *PSD, ///< [in] the PSD to integrate.
103 : size_t sz, ///< [in] the size of the PSD vector
104 : realT half = 0.5 ///< [in] [optional] controls if trapezoid (0.5) or mid-point (1.0) integration is
105 : ///< used. Do not use other values.
106 : )
107 : {
108 3 : realT var = 0;
109 :
110 3 : var = PSD[0];
111 :
112 3 : size_t i = 1;
113 10 : for( ; i < sz / 2; ++i )
114 : {
115 7 : var += PSD[i];
116 7 : var += PSD[sz - i];
117 : }
118 3 : var += half * PSD[i]; // The mid-point is double. It is also the endpoint of integration from each side, so it
119 : // would enter twice, hence once here.
120 :
121 3 : var *= df;
122 :
123 3 : return var;
124 : }
125 :
126 : /// Calculate the variance of a 1-D PSD
127 : /** By default uses trapezoid rule integration. This can be changed to mid-point integration.
128 : *
129 : * If f.back() < 0, then a 2-sided PSD in FFT storage order is assumed. Otherwise, PSD is treated as 1-sided.
130 : *
131 : * \returns the variance of a PSD (the integral).
132 : *
133 : * \tparam realT the real floating point type
134 : *
135 : * \test Scenario: calculating variance from a 1D PSD. \ref tests_sigproc_psdUtils_psdVar_1D "[test doc]"
136 : */
137 : template <typename realT>
138 6 : realT psdVar( const std::vector<realT> &f, ///< [in] the frequency scale of the PSD.
139 : const std::vector<realT> &PSD, ///< [in] the PSD to integrate.
140 : realT half = 0.5 ///< [in] [optional] controls if trapezoid (0.5) or mid-point (1.0) integration is used.
141 : ///< Do not use other values.
142 : )
143 : {
144 6 : if( f.back() < 0 )
145 : {
146 3 : return psdVar2sided( f[1] - f[0], PSD.data(), PSD.size(), half );
147 : }
148 : else
149 : {
150 3 : return psdVar1sided( f[1] - f[0], PSD.data(), PSD.size(), half );
151 : }
152 : }
153 :
154 : /// Calculate the variance of a PSD
155 : /** By default uses trapezoid rule integration. This can be changed to mid-point integration.
156 : *
157 : * \overload
158 : *
159 : * \returns the variance of a PSD (the integral).
160 : *
161 : * \tparam realT the real floating point type
162 : */
163 : template <typename eigenArrT>
164 : typename eigenArrT::Scalar psdVarDisabled(
165 : eigenArrT &freq, ///< [in] the frequency scale of the PSD
166 : eigenArrT &PSD, ///< [in] the PSD to integrate.
167 : bool trap = true ///< [in] [optional] controls if trapezoid (true) or mid-point (false) integration is used.
168 : )
169 : {
170 : typename eigenArrT::Scalar half = 0.5;
171 : if( !trap )
172 : half = 1.0;
173 :
174 : typename eigenArrT::Scalar var = 0;
175 :
176 : var = half * PSD( 0, 0 );
177 :
178 : for( int i = 1; i < freq.rows() - 1; ++i )
179 : var += PSD( i, 0 );
180 :
181 : var += half * PSD( freq.rows() - 1, 0 );
182 :
183 : var *= ( freq( 1, 0 ) - freq( 0, 0 ) );
184 :
185 : return var;
186 : }
187 :
188 : /// Calculates the frequency sampling for a grid given maximum dimension and maximum frequency.
189 : /** The freq_sampling is
190 : * @f$ \Delta f = f_{max}/ (0.5*dim) @f$
191 : * where @f$ f_{max} = 1/(2\Delta t) @f$ is the maximum frequency and @f$ dim @f$ is the size of the grid.
192 : *
193 : * \param [in] dim is the size of the grid
194 : * \param [in] f_max is the maximum frequency of the grid
195 : *
196 : * \returns the sampling interval @f$ \delta f @f$
197 : *
198 : * \tparam realT is the real floating point type used for calculations.
199 : *
200 : */
201 : template <class realT>
202 4 : realT freq_sampling( size_t dim, realT f_max )
203 : {
204 4 : return ( f_max / ( 0.5 * dim ) );
205 : }
206 :
207 : #if 0
208 : ///Create a 1-D frequency grid
209 : /**
210 : * \param [out] vec the pre-allocated Eigen-type 1xN or Nx1 array, on return contains the frequency grid
211 : * \param [in] dt the temporal sampling of the time series
212 : * \param [in] inverse [optional] if true
213 : *
214 : * \tparam eigenArr the Eigen-like array type
215 : */
216 : template<typename eigenArr>
217 : void frequency_grid1D( eigenArr & vec,
218 : typename eigenArr::Scalar dt,
219 : bool inverse = false )
220 : {
221 : typename eigenArr::Index dim, dim_1, dim_2;
222 : typename eigenArr::Scalar df;
223 :
224 : dim_1 = vec.rows();
225 : dim_2 = vec.cols();
226 :
227 : dim = std::max(dim_1, dim_2);
228 :
229 : df = freq_sampling(dim, 0.5/dt);
230 :
231 : if( !inverse )
232 : {
233 : for(int ii=0; ii < ceil(0.5*(dim-1) + 1); ++ii)
234 : {
235 : vec(ii) = ii*df;
236 : }
237 :
238 : for(int ii=ceil(0.5*(dim-1)+1); ii < dim_1; ++ii)
239 : {
240 : vec(ii) = (ii-dim)*df;
241 : }
242 : }
243 : else
244 : {
245 : for(int ii=0; ii < dim; ++ii)
246 : {
247 : vec(ii) = df * ii / dim;
248 : }
249 : }
250 : }
251 : #endif
252 :
253 : /// Create a 1-D frequency grid
254 : /**
255 : *
256 : * \tparam realT a real floating point type
257 : * \tparam realParamT a real floating point type, convenience to avoid double-float confusion.
258 : *
259 : * \test Verify creation of a 1D frequency grid \ref tests_sigproc_psdUtils_frequencyGrid_1D "[test doc]"
260 : */
261 : template <typename realT, typename realParamT>
262 4 : int frequencyGrid(
263 : std::vector<realT> &vec, ///< [out] vec the pre-allocated vector, on return contains the frequency grid
264 : realParamT dt, ///< [in] dt the temporal sampling of the time series
265 : bool fftOrder = true ///< [in] fftOrder [optional] if true the frequency grid is in FFT order
266 : )
267 : {
268 4 : realT dtTT = dt;
269 :
270 4 : if( fftOrder )
271 : {
272 2 : if( vec.size() % 2 == 1 )
273 : {
274 0 : internal::mxlib_error_report(error_t::invalidarg,"Frequency scale can't be odd-sized for FFT order" );
275 0 : return -1;
276 : }
277 :
278 2 : realT df = ( 1.0 / dtTT ) / ( (realT)vec.size() );
279 :
280 14 : for( ssize_t ii = 0; ii < ceil( 0.5 * ( vec.size() - 1 ) + 1 ); ++ii )
281 : {
282 12 : vec[ii] = ii * df;
283 : }
284 :
285 10 : for( ssize_t ii = ceil( 0.5 * ( vec.size() - 1 ) + 1 ); ii < vec.size(); ++ii )
286 : {
287 8 : vec[ii] = ( ii - (ssize_t)vec.size() ) * df;
288 : }
289 :
290 2 : return 0;
291 : }
292 : else
293 : {
294 2 : if( vec.size() % 2 == 0 )
295 : {
296 1 : realT df = ( 0.5 / dtTT ) / ( (realT)vec.size() - 1 );
297 7 : for( int ii = 0; ii < vec.size(); ++ii )
298 : {
299 6 : vec[ii] = df * ii;
300 : }
301 :
302 1 : return 0;
303 : }
304 : else
305 : {
306 1 : realT df = ( 0.5 / dt ) / ( (realT)vec.size() );
307 6 : for( int ii = 0; ii < vec.size(); ++ii )
308 : {
309 5 : vec[ii] = df * ( ii + 1 );
310 : }
311 :
312 1 : return 0;
313 : }
314 : }
315 : }
316 :
317 : /// Create a 2-D frequency grid
318 : template <typename eigenArr, typename realParamT>
319 4 : void frequencyGrid( eigenArr &arr, realParamT drT, eigenArr *k_x, eigenArr *k_y )
320 : {
321 4 : typename eigenArr::Scalar dr = drT;
322 :
323 : typename eigenArr::Index dim_1, dim_2;
324 : typename eigenArr::Scalar k_1, k_2, df;
325 :
326 4 : dim_1 = arr.rows();
327 4 : dim_2 = arr.cols();
328 :
329 4 : if( k_x )
330 0 : k_x->resize( dim_1, dim_2 );
331 4 : if( k_y )
332 0 : k_y->resize( dim_1, dim_2 );
333 :
334 4 : df = freq_sampling( std::max( dim_1, dim_2 ), 0.5 / dr );
335 :
336 104 : for( int ii = 0; ii < 0.5 * ( dim_1 - 1 ) + 1; ++ii )
337 : {
338 100 : k_1 = ii * df;
339 2856 : for( int jj = 0; jj < 0.5 * ( dim_2 - 1 ) + 1; ++jj )
340 : {
341 2756 : k_2 = jj * df;
342 :
343 2756 : arr( ii, jj ) = sqrt( k_1 * k_1 + k_2 * k_2 );
344 :
345 2756 : if( k_x )
346 0 : ( *k_x )( ii, jj ) = k_1;
347 2756 : if( k_x )
348 0 : ( *k_y )( ii, jj ) = k_2;
349 : }
350 :
351 2756 : for( int jj = 0.5 * ( dim_2 - 1 ) + 1; jj < dim_2; ++jj )
352 : {
353 2656 : k_2 = ( jj - dim_2 ) * df;
354 :
355 2656 : arr( ii, jj ) = sqrt( k_1 * k_1 + k_2 * k_2 );
356 :
357 2656 : if( k_x )
358 0 : ( *k_x )( ii, jj ) = k_1;
359 2656 : if( k_x )
360 0 : ( *k_y )( ii, jj ) = k_2;
361 : }
362 : }
363 :
364 100 : for( int ii = 0.5 * ( dim_1 - 1 ) + 1; ii < dim_1; ++ii )
365 : {
366 96 : k_1 = ( ii - dim_1 ) * df;
367 2752 : for( int jj = 0; jj < 0.5 * ( dim_2 - 1 ) + 1; ++jj )
368 : {
369 2656 : k_2 = jj * df;
370 :
371 2656 : arr( ii, jj ) = sqrt( k_1 * k_1 + k_2 * k_2 );
372 :
373 2656 : if( k_x )
374 0 : ( *k_x )( ii, jj ) = k_1;
375 2656 : if( k_x )
376 0 : ( *k_y )( ii, jj ) = k_2;
377 : }
378 :
379 2656 : for( int jj = 0.5 * ( dim_2 - 1 ) + 1; jj < dim_2; ++jj )
380 : {
381 2560 : k_2 = ( jj - dim_2 ) * df;
382 :
383 2560 : arr( ii, jj ) = sqrt( k_1 * k_1 + k_2 * k_2 );
384 :
385 2560 : if( k_x )
386 0 : ( *k_x )( ii, jj ) = k_1;
387 2560 : if( k_x )
388 0 : ( *k_y )( ii, jj ) = k_2;
389 : }
390 : }
391 4 : }
392 :
393 : /// Create a frequency grid
394 : template <typename eigenArr>
395 4 : void frequencyGrid( eigenArr &arr, typename eigenArr::Scalar dt )
396 : {
397 4 : frequencyGrid( arr, dt, (eigenArr *)0, (eigenArr *)0 );
398 4 : }
399 :
400 : /// Create a frequency grid
401 : template <typename eigenArr>
402 : void frequencyGrid( eigenArr &arr, typename eigenArr::Scalar dt, eigenArr &k_x, eigenArr &k_y )
403 : {
404 : frequencyGrid( arr, dt, &k_x, &k_y );
405 : }
406 :
407 : /// Calculate the normalization for a 1-D @f$ 1/|f|^\alpha @f$ PSD.
408 : /**
409 : * \param [in] fmin is the minimum non-zero absolute value of frequency
410 : * \param [in] fmax is the maximum absolute value of frequencey
411 : * \param [in] alpha is the power-law exponent, by convention @f$ \alpha > 0 @f$.
412 : *
413 : * \returns the normalization for a 2-sided power law PSD.
414 : *
415 : * \tparam realT is the real floating point type used for calculations.
416 : */
417 : template <typename realT>
418 : realT oneoverf_norm( realT fmin, realT fmax, realT alpha )
419 : {
420 : realT integ = 2 * ( pow( fmax, -1.0 * alpha + 1.0 ) - pow( fmin, -1.0 * alpha + 1.0 ) ) / ( -1.0 * alpha + 1.0 );
421 :
422 : return 1 / integ;
423 : }
424 :
425 : /// Calculate the normalization for a 2-D @f$ 1/|k|^\alpha @f$ PSD.
426 : /**
427 : * \param [in] kmin is the minimum non-zero absolute value of frequency
428 : * \param [in] kmax is the maximum absolute value of frequencey
429 : * \param [in] alpha is the power-law exponent, by convention @f$ \alpha > 0 @f$.
430 : *
431 : * \returns the normalization for a 2-D, 2-sided power law PSD.
432 : *
433 : * \tparam realT is the real floating point type used for calculations.
434 : */
435 : template <typename realT>
436 : realT oneoverk_norm( realT kmin, realT kmax, realT alpha )
437 : {
438 : realT integ = 2 * ( pow( kmax, -1 * alpha + 2.0 ) - pow( kmin, -1.0 * alpha + 2.0 ) ) / ( -1 * alpha + 2.0 );
439 :
440 : return 1 / integ;
441 : }
442 :
443 : /// Normalize a 1-D PSD to have a given variance
444 : /** A frequency range can be specified to calculate the norm, otherwise f[0] to f[f.size()-1] is the range. The entire
445 : * PSD is normalized regardless.
446 : *
447 : * \tparam floatT the floating point type of the PSD.
448 : * \tparam floatParamT a floating point type, convenience to avoid double-float cofusion.
449 : *
450 : * \test Verify scaling and normalization of augment1SidedPSD \ref tests_sigproc_psdUtils_augment1SidedPSD "[test doc]"
451 : */
452 : template <typename floatT, typename floatParamT>
453 2049 : int normPSD( std::vector<floatT> &psd, ///< [in.out] the PSD to normalize, will be altered.
454 : std::vector<floatT> &f, ///< [in] the frequency points for the PSD
455 : floatParamT normT, ///< [in] the desired total variance (or integral) of the PSD.
456 : floatT fmin = std::numeric_limits<floatT>::min(), ///< [in] [optiona] the minimum frequency of the range
457 : ///< over which to normalize.
458 : floatT fmax = std::numeric_limits<floatT>::max() ///< [in] [optiona] the maximum frequency of the range
459 : ///< over which to normalize.
460 : )
461 : {
462 2049 : floatT norm = normT;
463 :
464 2049 : floatT s = 0; // accumulate
465 :
466 : // Check if inside-the-loop branch is needed
467 2049 : if( fmin != std::numeric_limits<floatT>::min() || fmax != std::numeric_limits<floatT>::max() )
468 : {
469 137089 : for( size_t i = 0; i < psd.size(); ++i )
470 : {
471 135041 : if( fabs( f[i] ) < fmin || fabs( f[i] ) > fmax )
472 0 : continue;
473 135041 : s += psd[i];
474 : }
475 : }
476 : else
477 : {
478 6 : for( size_t i = 0; i < psd.size(); ++i )
479 : {
480 5 : s += psd[i];
481 : }
482 : }
483 :
484 2049 : s *= ( f[1] - f[0] );
485 :
486 137095 : for( size_t i = 0; i < psd.size(); ++i )
487 135046 : psd[i] *= norm / s;
488 :
489 2049 : return 0;
490 : }
491 :
492 : /// Normalize a 2-D PSD to have a given variance
493 : /** A frequency range can be specified for calculating the norm, otherwise the entire PSD is used. The entire PSD is
494 : * normalized regardless.
495 : *
496 : * \tparam floatT the floating point type of the PSD.
497 : * \tparam floatParamT a floating point type, convenience to avoid double-float cofusion.
498 : *
499 : * \test Verify scaling and normalization of augment1SidedPSD \ref tests_sigproc_psdUtils_augment1SidedPSD "[test doc]"
500 : */
501 : template <typename floatT, typename floatParamT>
502 : floatT
503 4 : normPSD( Eigen::Array<floatT, Eigen::Dynamic, Eigen::Dynamic> &psd, ///< [in.out] the PSD to normalize, will be altered.
504 : Eigen::Array<floatT, Eigen::Dynamic, Eigen::Dynamic> &k, ///< [in] the frequency grid for psd.
505 : floatParamT normT, ///< [in] the desired total variance (or integral) of the PSD.
506 : floatT kmin = std::numeric_limits<floatT>::min(), ///< [in] [optiona] the minimum frequency of the range over
507 : ///< which to normalize.
508 : floatT kmax = std::numeric_limits<floatT>::max() ///< [in] [optiona] the maximum frequency of the range over
509 : ///< which to normalize.
510 : )
511 : {
512 4 : floatT norm = normT;
513 :
514 : floatT dk1, dk2;
515 :
516 4 : if( k.rows() > 1 )
517 4 : dk1 = k( 1, 0 ) - k( 0, 0 );
518 : else
519 0 : dk1 = 1;
520 :
521 4 : if( k.cols() > 1 )
522 4 : dk2 = k( 0, 1 ) - k( 0, 0 );
523 : else
524 0 : dk2 = 1;
525 :
526 4 : floatT s = 0;
527 :
528 : // Check if inside-the-loop branch is needed
529 4 : if( kmin != std::numeric_limits<floatT>::min() || kmax != std::numeric_limits<floatT>::max() )
530 : {
531 0 : for( int c = 0; c < psd.cols(); ++c )
532 : {
533 0 : for( int r = 0; r < psd.rows(); ++r )
534 : {
535 0 : if( fabs( k( r, c ) ) < kmin || fabs( k( r, c ) ) > kmax )
536 0 : continue;
537 0 : s += psd( r, c );
538 : }
539 : }
540 : }
541 : else
542 : {
543 196 : for( int c = 0; c < psd.cols(); ++c )
544 : {
545 10432 : for( int r = 0; r < psd.rows(); ++r )
546 : {
547 10240 : s += psd( r, c );
548 : }
549 : }
550 : }
551 :
552 4 : s *= dk1 * dk2;
553 :
554 196 : for( int c = 0; c < psd.cols(); ++c )
555 : {
556 10432 : for( int r = 0; r < psd.rows(); ++r )
557 : {
558 10240 : psd( r, c ) *= norm / s;
559 : }
560 : }
561 :
562 4 : return 0;
563 : }
564 :
565 : /// Generates a @f$ 1/|f|^\alpha @f$ power spectrum
566 : /**
567 : * Populates an Eigen array with
568 : * \f[
569 : * P(|f| = 0) = 0
570 : * \f]
571 : * \f[
572 : * P(|f| > 0) = \frac{\beta}{|f|^{\alpha}}
573 : * \f]
574 : *
575 : *
576 : * \param [out] psd is the array to populate
577 : * \param [in] freq is a frequency grid, must be the same logical size as psd
578 : * \param [in] alpha is the power law exponent, by convention @f$ alpha > 0 @f$.
579 : * \param [in] beta [optional is a normalization constant to multiply the raw spectrum by. If beta==-1 (default) then
580 : * the PSD is normalized using \ref oneoverf_norm.
581 : *
582 : * \tparam eigenArrp is the Eigen-like array type of the psd
583 : * \tparam eigenArrf is the Eigen-like array type of the frequency grid
584 : */
585 : template <typename eigenArrp, typename eigenArrf>
586 : void oneoverf_psd( eigenArrp &psd,
587 : eigenArrf &freq,
588 : typename eigenArrp::Scalar alpha,
589 : typename eigenArrp::Scalar beta = -1 )
590 : {
591 : typedef typename eigenArrp::Scalar Scalar;
592 :
593 : typename eigenArrp::Index dim_1, dim_2;
594 : Scalar f_x, f_y, p;
595 :
596 : dim_1 = psd.rows();
597 : dim_2 = psd.cols();
598 :
599 : if( beta == -1 )
600 : {
601 : Scalar fmin;
602 : Scalar fmax;
603 :
604 : fmax = freq.abs().maxCoeff();
605 :
606 : // Find minimum non-zero Coeff.
607 : fmin = ( freq.abs() > 0 ).select( freq.abs(), freq.abs() + fmax ).minCoeff();
608 :
609 : beta = oneoverf_norm( fmin, fmax, alpha );
610 : }
611 :
612 : for( int ii = 0; ii < dim_1; ++ii )
613 : {
614 : for( int jj = 0; jj < dim_2; ++jj )
615 : {
616 : if( freq( ii, jj ) == 0 )
617 : {
618 : p = 0;
619 : }
620 : else
621 : {
622 : p = beta / std::pow( std::abs( freq( ii, jj ) ), alpha );
623 : }
624 : psd( ii, jj ) = p;
625 : }
626 : }
627 : }
628 :
629 : /// Generate a 1-D von Karman power spectrum
630 : /**
631 : * Populates an Eigen array with
632 : *
633 : * \f[
634 : * P(f) = \frac{\beta}{ (f^2 + (1/T_0)^2)^{\alpha/2}} e^{ - f^2 t_0^2}
635 : * \f]
636 : *
637 : * If you set \f$ T_0 \le 0 \f$ and \f$ t_0 = 0\f$ this reverts to a simple \f$ 1/f^\alpha \f$ law (i.e.
638 : * it treats this as infinite outer scale and inner scale).
639 : *
640 : * \returns error_t::noerror on success
641 : *
642 : * \tparam floatT a floating point
643 : */
644 : template <typename floatT,
645 : typename floatfT,
646 : typename alphaT,
647 : typename T0T = double,
648 : typename t0T = double,
649 : typename betaT = double>
650 : mx::error_t vonKarmanPSD( std::vector<floatT> &psd, ///< [out] the PSD vector, will be resized.
651 : std::vector<floatfT> &f, ///< [in] the frequency vector
652 : alphaT alpha, ///< [in] the exponent, by convention @f$ alpha > 0 @f$.
653 : T0T T0 = 0, ///< [in] the outer scale, default is 0 (not used).
654 : t0T t0 = 0, ///< [in] the inner scale, default is 0 (not used).
655 : betaT beta = 1 ///< [in] the scaling constant, default is 1
656 : )
657 : {
658 :
659 : floatT T02;
660 : if( T0 > 0 )
661 : {
662 : T02 = 1.0 / ( T0 * T0 );
663 : }
664 : else
665 : {
666 : T02 = 0;
667 : }
668 :
669 : floatT sqrt_alpha = 0.5 * alpha;
670 :
671 : floatT _beta;
672 : if( beta <= 0 )
673 : {
674 : _beta = 1;
675 : }
676 : else
677 : {
678 : _beta = beta;
679 : }
680 :
681 : psd.resize( f.size() );
682 :
683 : if( t0 > 0 )
684 : {
685 : floatT _t0 = static_cast<floatT>( t0 );
686 : for( size_t i = 0; i < f.size(); ++i )
687 : {
688 : psd[i] = _beta / pow( pow( f[i], 2 ) + T02, sqrt_alpha ) * exp( -1 * pow( f[i] * _t0, 2 ) );
689 : }
690 : }
691 : else
692 : {
693 : for( size_t i = 0; i < f.size(); ++i )
694 : {
695 : psd[i] = _beta / pow( pow( f[i], 2 ) + T02, sqrt_alpha );
696 : }
697 : }
698 :
699 : return mx::error_t::noerror;
700 : }
701 :
702 : /// Generate a 1-D "knee" PSD
703 : /**
704 : * Populates an Eigen array with
705 : *
706 : * \f[
707 : * P(f) = \frac{\beta}{ 1 + (f/f_n)^{\alpha}}
708 : * \f]
709 : *
710 : * If you set \f$ T_0 \le 0 \f$ and \f$ t_0 = 0\f$ this reverts to a simple \f$ 1/f^\alpha \f$ law (i.e.
711 : * it treats this as infinite outer scale and inner scale).
712 : *
713 : * \tparam floatT a floating point
714 : */
715 : template <typename floatT>
716 : int kneePSD( std::vector<floatT> &psd, ///< [out] the PSD vector, will be resized.
717 : std::vector<floatT> &f, ///< [in] the frequency vector
718 : floatT beta, ///< [in] the scaling constant
719 : floatT fn, ///< [in] the knee frequency
720 : floatT alpha ///< [in] the exponent, by convention @f$ alpha > 0 @f$.
721 : )
722 : {
723 :
724 : psd.resize( f.size() );
725 :
726 : for( int i = 0; i < f.size(); ++i )
727 : {
728 : floatT p = beta / ( 1 + pow( f[i] / fn, alpha ) );
729 : psd[i] = p;
730 : }
731 :
732 : return 0;
733 : }
734 :
735 : /// Generates a von Karman power spectrum
736 : /**
737 : * Populates an Eigen array with
738 : *
739 : * \f[
740 : * P(k) = \frac{\beta}{ (k^2 + (1/L_0)^2)^{\alpha/2}} e^{ - k^2 l_0^2}
741 : * \f]
742 : *
743 : * If you set \f$ L_0 \le 0 \f$ and \f$ l_0 = 0\f$ this reverts to a simple \f$ 1/f^\alpha \f$ law (i.e.
744 : * it treats this as infinite outer scale and inner scale).
745 : *
746 : * \param [out] psd is the array to populate, allocated.
747 : * \param [in] freq is a frequency grid, must be the same logical size as psd
748 : * \param [in] alpha is the power law exponent, by convention @f$ alpha > 0 @f$.
749 : * \param [in] L0 [optional] is the outer scale.
750 : * \param [in] l0 [optional] is the inner scale.
751 : * \param [in] beta [optional] is a normalization constant to multiply the raw spectrum by. If beta==-1 (default)
752 : * then the PSD is normalized using \ref oneoverf_norm.
753 : *
754 : * \tparam eigenArrp is the Eigen array type of the psd
755 : * \tparam eigenArrf is the Eigen array type of the frequency grid
756 : */
757 : template <typename eigenArrp, typename eigenArrf, typename alphaT, typename L0T, typename l0T, typename betaT>
758 : void vonKarmanPSD( eigenArrp &psd, eigenArrf &freq, alphaT alpha, L0T L0 = 0, l0T l0 = 0, betaT beta = -1 )
759 : {
760 : typedef typename eigenArrp::Scalar Scalar;
761 :
762 : typename eigenArrp::Index dim_1, dim_2;
763 : Scalar p;
764 :
765 : dim_1 = psd.rows();
766 : dim_2 = psd.cols();
767 :
768 : Scalar _beta;
769 :
770 : if( beta == -1 )
771 : {
772 : Scalar fmin;
773 : Scalar fmax;
774 :
775 : fmax = freq.abs().maxCoeff();
776 :
777 : // Find minimum non-zero Coeff.
778 : fmin = ( freq.abs() > 0 ).select( freq.abs(), freq.abs() + fmax ).minCoeff();
779 :
780 : _beta = beta = oneoverf_norm( fmin, fmax, static_cast<Scalar>( alpha ) );
781 : }
782 : else
783 : _beta = static_cast<Scalar>( beta );
784 :
785 : Scalar L02;
786 : if( L0 > 0 )
787 : L02 = 1.0 / ( L0 * L0 );
788 : else
789 : L02 = 0;
790 :
791 : Scalar sqrt_alpha = 0.5 * alpha; // std::sqrt(alpha);
792 :
793 : for( int ii = 0; ii < dim_1; ++ii )
794 : {
795 : for( int jj = 0; jj < dim_2; ++jj )
796 : {
797 : if( freq( ii, jj ) == 0 && L02 == 0 )
798 : {
799 : p = 0;
800 : }
801 : else
802 : {
803 : p = _beta / pow( pow( freq( ii, jj ), 2 ) + L02, sqrt_alpha );
804 : if( l0 > 0 )
805 : p *= exp( -1 * pow( freq( ii, jj ) * static_cast<Scalar>( l0 ), 2 ) );
806 : }
807 : psd( ii, jj ) = p;
808 : }
809 : }
810 : }
811 :
812 : /// Augment a 1-sided PSD to standard 2-sided FFT form.
813 : /** Allocates psdTwoSided to hold a flipped copy of psdOneSided.
814 : * Default assumes that psdOneSided[0] corresponds to 0 frequency,
815 : * but this can be changed by setting zeroFreq to a non-zero value.
816 : * In this case psdTwoSided[0] is set to 0, and the augmented psd
817 : * is shifted by 1.
818 : *
819 : * To illustrate, the bins are re-ordered as:
820 : * \verbatim
821 : * {1,2,3,4,5} --> {0,1,2,3,4,5,-4,-3,-2,-1}
822 : * \endverbatim
823 : *
824 : * The output is scaled so that the total power remains the same. The 0-freq and
825 : * Nyquist freq are not scaled.
826 : *
827 : *
828 : * Entries in psdOneSided are cast to the value_type of psdTwoSided,
829 : * for instance to allow for conversion to complex type.
830 : *
831 : * \test Verify scaling and normalization of augment1SidedPSD \ref tests_sigproc_psdUtils_augment1SidedPSD "[test doc]"
832 : */
833 : template <typename vectorTout, typename vectorTin>
834 3 : void augment1SidedPSD(
835 : vectorTout &psdTwoSided, ///< [out] on return contains the FFT storage order copy of psdOneSided.
836 : vectorTin &psdOneSided, ///< [in] the one-sided PSD to augment
837 : bool addZeroFreq =
838 : false, ///< [in] [optional] set to true if psdOneSided does not contain a zero frequency component.
839 : typename vectorTin::value_type scale = 0.5 ///< [in] [optional] value to scale the input by when copying to the
840 : ///< output. The default 0.5 re-normalizes for a 2-sided PSD.
841 : )
842 : {
843 : typedef typename vectorTout::value_type outT;
844 :
845 3 : bool needZero = 1;
846 :
847 : size_t N;
848 :
849 3 : if( addZeroFreq == 0 )
850 : {
851 3 : needZero = 0;
852 3 : N = 2 * psdOneSided.size() - 2;
853 : }
854 : else
855 : {
856 0 : N = 2 * psdOneSided.size();
857 : }
858 :
859 3 : psdTwoSided.resize( N );
860 :
861 : // First set the 0-freq point
862 3 : if( needZero )
863 : {
864 0 : psdTwoSided[0] = outT( 0.0 );
865 : }
866 : else
867 : {
868 3 : psdTwoSided[0] = outT( psdOneSided[0] );
869 : }
870 :
871 : // Now set all the rest.
872 : unsigned long i;
873 3076 : for( i = 0; i < psdOneSided.size() - 1 - ( 1 - needZero ); ++i )
874 : {
875 3073 : psdTwoSided[i + 1] = outT( psdOneSided[i + ( 1 - needZero )] * scale );
876 3073 : psdTwoSided[i + psdOneSided.size() + needZero] = outT( psdOneSided[psdOneSided.size() - 2 - i] * scale );
877 : }
878 3 : psdTwoSided[i + 1] = outT( psdOneSided[i + ( 1 - needZero )] );
879 3 : }
880 :
881 : /// Augment a 1-sided frequency scale to standard FFT form.
882 : /** Allocates freqTwoSided to hold a flipped copy of freqOneSided.
883 : * If freqOneSided[0] is not 0, freqTwoSided[0] is set to 0, and the augmented
884 : * frequency scale is shifted by 1.
885 : *
886 : * Example:
887 : *
888 : * {1,2,3,4,5} --> {0,1,2,3,4,5,-4,-3,-2,-1}
889 : *
890 : * \test Verify scaling and normalization of augment1SidedPSD \ref tests_sigproc_psdUtils_augment1SidedPSD "[test doc]"
891 : */
892 : template <typename T>
893 5 : void augment1SidedPSDFreq(
894 : std::vector<T> &freqTwoSided, ///< [out] on return contains the FFT storage order copy of freqOneSided.
895 : std::vector<T> &freqOneSided ///< [in] the one-sided frequency scale to augment
896 : )
897 : {
898 5 : int needZero = 1;
899 :
900 : size_t N;
901 :
902 5 : if( freqOneSided[0] != 0 )
903 : {
904 0 : N = 2 * freqOneSided.size();
905 : }
906 : else
907 : {
908 5 : needZero = 0;
909 5 : N = 2 * freqOneSided.size() - 2;
910 : }
911 :
912 5 : freqTwoSided.resize( N );
913 :
914 5 : if( needZero )
915 : {
916 0 : freqTwoSided[0] = 0.0;
917 : }
918 : else
919 : {
920 5 : freqTwoSided[0] = freqOneSided[0]; // 0
921 : }
922 :
923 : int i;
924 3140 : for( i = 0; i < freqOneSided.size() - 1 - ( 1 - needZero ); ++i )
925 : {
926 3135 : freqTwoSided[i + 1] = freqOneSided[i + ( 1 - needZero )];
927 3135 : freqTwoSided[i + freqOneSided.size() + needZero] = -freqOneSided[freqOneSided.size() - 2 - i];
928 : }
929 5 : freqTwoSided[i + 1] = freqOneSided[i + ( 1 - needZero )];
930 5 : }
931 :
932 : /// Rebin a PSD, including its frequency scale, to a larger frequency bin size (fewer bins)
933 : /** The rebinning uses trapezoid integration within bins to ensure minimum signal loss.
934 : *
935 : * Maintains DFT sampling. That is, if initial frequency grid is 0,0.1,0.2...
936 : * and the binSize is 1.0, the new grid will be 0,1,2 (as opposed to 0.5, 1.5, 2.5).
937 : *
938 : * This introduces a question of what to do with first half-bin, which includes 0. It can be
939 : * integrated (binAtZero = true, the default). This may cause inaccurate behavior if the value of the PSD when
940 : * f=0 is important (e.g. when analyzing correlated noise), so setting binAtZero=false causes the f=0 value to be
941 : * copied (using the nearest neighbor if no f=0 point is in the input.
942 : *
943 : * The last half bin is always integrated.
944 : *
945 : * The output is variance normalized to match the input variance.
946 : *
947 : * \tparam realT the real floating point type
948 : */
949 : template <typename realT>
950 : int rebin1SidedPSD( std::vector<realT> &binFreq, ///< [out] the binned frequency scale, resized.
951 : std::vector<realT> &binPSD, ///< [out] the binned PSD, resized.
952 : std::vector<realT> &freq, ///< [in] the frequency scale of the PSD to bin.
953 : std::vector<realT> &PSD, ///< [in] the PSD to bin.
954 : realT binSize, ///< [in] in same units as freq
955 : bool binAtZero = true ///< [in] [optional] controls whether the zero point is binned or copied.
956 : )
957 : {
958 : binFreq.clear();
959 : binPSD.clear();
960 :
961 : realT sumPSD = 0;
962 : realT startFreq = 0;
963 : realT sumFreq = 0;
964 : int nSum = 0;
965 :
966 : int i = 0;
967 :
968 : realT df = freq[1] - freq[0];
969 :
970 : // Now move to first bin
971 : while( freq[i] <= 0.5 * binSize + 0.5 * df )
972 : {
973 : sumPSD += PSD[i];
974 : ++nSum;
975 : ++i;
976 : if( i >= freq.size() )
977 : break;
978 : }
979 :
980 : if( !binAtZero )
981 : {
982 : binFreq.push_back( 0 );
983 : binPSD.push_back( PSD[0] );
984 : }
985 : else
986 : {
987 : binFreq.push_back( 0 );
988 : binPSD.push_back( sumPSD / nSum );
989 : }
990 :
991 : --i;
992 : startFreq = freq[i];
993 : nSum = 0;
994 : sumFreq = 0;
995 : sumPSD = 0;
996 :
997 : while( i < freq.size() )
998 : {
999 : realT sc = 0.5; // First one is multiplied by 1/2 for trapezoid rule.
1000 : while( freq[i] - startFreq + 0.5 * df < binSize )
1001 : {
1002 : sumFreq += freq[i];
1003 : sumPSD += sc * PSD[i];
1004 : sc = 1.0;
1005 : ++nSum;
1006 :
1007 : ++i;
1008 : if( i >= freq.size() - 1 )
1009 : break; // break 1 element early so last point is mult by 0.5
1010 : }
1011 :
1012 : if( i < freq.size() )
1013 : {
1014 : sumFreq += freq[i];
1015 : sumPSD += 0.5 * PSD[i]; // last one is multiplied by 1/2 for trapezoid rule
1016 : ++nSum;
1017 : ++i;
1018 : }
1019 :
1020 : // Check if this is last point
1021 : if( i < freq.size() )
1022 : {
1023 : binFreq.push_back( sumFreq / nSum );
1024 : }
1025 : else
1026 : {
1027 : // last point frequencyis not averaged.
1028 : binFreq.push_back( freq[freq.size() - 1] );
1029 : }
1030 :
1031 : binPSD.push_back( sumPSD / ( nSum - 1 ) );
1032 :
1033 : sumFreq = 0;
1034 : sumPSD = 0;
1035 : nSum = 0;
1036 : if( i >= freq.size() )
1037 : break;
1038 :
1039 : --i; // Step back one, so averages are edge to edge.
1040 :
1041 : startFreq = freq[i];
1042 : }
1043 :
1044 : // Now normalize variance
1045 : realT var = psdVar( freq, PSD );
1046 : realT binv = psdVar( binFreq, binPSD );
1047 :
1048 : for( int i = 0; i < binFreq.size(); ++i )
1049 : binPSD[i] *= var / binv;
1050 :
1051 : return 0;
1052 : }
1053 : ///@}
1054 :
1055 : } // namespace sigproc
1056 : } // namespace mx
1057 :
1058 : #endif // psdUtils_hpp
|
