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