Line data Source code
1 : /** \file psdFilter.hpp
2 : * \brief Declares and defines a class for filtering with PSDs
3 : * \ingroup signal_processing_files
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2015, 2016, 2017, 2018, 2019, 2020 Jared R. Males (jaredmales@gmail.com)
10 : //
11 : // This file is part of mxlib.
12 : //
13 : // mxlib is free software: you can redistribute it and/or modify
14 : // it under the terms of the GNU General Public License as published by
15 : // the Free Software Foundation, either version 3 of the License, or
16 : // (at your option) any later version.
17 : //
18 : // mxlib is distributed in the hope that it will be useful,
19 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 : // GNU General Public License for more details.
22 : //
23 : // You should have received a copy of the GNU General Public License
24 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
25 : //***********************************************************************//
26 :
27 : #ifndef psdFilter_hpp
28 : #define psdFilter_hpp
29 :
30 : #include <vector>
31 : #include <complex>
32 : #include <Eigen/Dense>
33 :
34 : #include "../mxlib.hpp"
35 : #include "../math/ft/fftT.hpp"
36 : #include "../improc/eigenCube.hpp"
37 :
38 : namespace mx
39 : {
40 : namespace sigproc
41 : {
42 :
43 : namespace psdFilterTypes
44 : {
45 :
46 : /// Types for different ranks in psdFilter
47 : template <typename realT, size_t rank>
48 : struct arrayT;
49 :
50 : template <typename realT>
51 : struct arrayT<realT, 1>
52 : {
53 : typedef std::vector<realT> realArrayT;
54 : typedef std::vector<realT> *realArrayMapT;
55 : typedef std::vector<std::complex<realT>> complexArrayT;
56 :
57 : static void clear( realArrayT &arr )
58 : {
59 : arr.clear();
60 : }
61 :
62 3 : static void clear( complexArrayT &arr )
63 : {
64 3 : arr.clear();
65 3 : }
66 : };
67 :
68 : template <typename realT>
69 : struct arrayT<realT, 2>
70 : {
71 : typedef Eigen::Array<realT, Eigen::Dynamic, Eigen::Dynamic> realArrayT;
72 : typedef Eigen::Map<Eigen::Array<realT, Eigen::Dynamic, Eigen::Dynamic>> realArrayMapT;
73 :
74 : typedef Eigen::Array<std::complex<realT>, Eigen::Dynamic, Eigen::Dynamic> complexArrayT;
75 :
76 : static void clear( realArrayT &arr )
77 : {
78 : arr.resize( 0, 0 );
79 : }
80 :
81 3 : static void clear( complexArrayT &arr )
82 : {
83 3 : arr.resize( 0, 0 );
84 3 : }
85 : };
86 :
87 : template <typename realT>
88 : struct arrayT<realT, 3>
89 : {
90 : typedef improc::eigenCube<realT> realArrayT;
91 : typedef improc::eigenCube<realT> *realArrayMapT;
92 : typedef improc::eigenCube<std::complex<realT>> complexArrayT;
93 :
94 : static void clear( realArrayT &arr )
95 : {
96 : arr.resize( 0, 0, 0 );
97 : }
98 :
99 3 : static void clear( complexArrayT &arr )
100 : {
101 3 : arr.resize( 0, 0, 0 );
102 3 : }
103 : };
104 :
105 : } // namespace psdFilterTypes
106 :
107 : // Forward declaration
108 : template <typename _realT, size_t rank, int cuda = 0>
109 : class psdFilter;
110 :
111 : /** \defgroup psd_filter PSD Filter
112 : * \brief Filtering with a PSD to generate correlated noise.
113 : *
114 : * \ingroup psds
115 : */
116 :
117 : /// A class for filtering noise with PSDs
118 : /** The square-root of the PSD is maintained by this class, either as a pointer to an external array or using internally
119 : * allocated memory (which will be de-allocated on destruction).
120 : *
121 : * PSD Requirements:
122 : * - the PSD must be in FFT storage order form. That means including negative frequencies reversed from the end of the
123 : * array.
124 : * - the PSD used for this needs to be normalized properly, \ref psds "according to the mxlib standard", to produce
125 : * filtered noise with the correct statistics.
126 : *
127 : *
128 : * Array type varies based on rank.
129 : * - For rank==1, the array type is std::vector<realT>
130 : * - for rank==2, the array type is Eigen::Array<realT, -1, -1>
131 : * - for rank==3, the array type is mx::improc::eigenCube<realT>
132 : * and likewise for the complex types.
133 : *
134 : * \tparam _realT real floating type
135 : * \tparam _rank the rank, or dimension, of the PSD
136 : *
137 : * \ingroup psd_filter
138 : *
139 : * \todo once fftT has a plan interface with pointers for working memory, use it.
140 : *
141 : */
142 : template <typename _realT, size_t _rank>
143 : class psdFilter<_realT, _rank, 0>
144 : {
145 : public:
146 : typedef _realT realT; ///< Real floating point type
147 : typedef std::complex<_realT> complexT; ///< Complex floating point type.
148 :
149 : static const size_t rank = _rank;
150 :
151 : typedef typename psdFilterTypes::arrayT<realT, rank>::realArrayT
152 : realArrayT; ///< std::vector for rank==1, Eigen::Array for rank==2, eigenCube for rank==3.
153 : typedef typename psdFilterTypes::arrayT<realT, rank>::realArrayMapT
154 : realArrayMapT; ///< std::vector for rank==1, Eigen::Map for rank==2, eigenCube for rank==3.
155 : typedef typename psdFilterTypes::arrayT<realT, rank>::complexArrayT
156 : complexArrayT; ///< std::vector for rank==1, Eigen::Array for rank==2, eigenCube for rank==3.
157 :
158 : protected:
159 : int m_rows{ 0 }; ///< The number of rows in the filter, and the required number of rows in the noise array.
160 : int m_cols{ 0 }; ///< The number of columns in the filter, and the required number of columns in the noise array.
161 : int m_planes{ 0 }; ///< Then number of planes in the filter.
162 :
163 : realT m_dFreq1{ 1.0 }; ///< The frequency scaling of the x-dimension. Used to scale the output.
164 : realT m_dFreq2{ 1.0 }; ///< The frequency scaling of the y-dimension. Used to scale the output.
165 : realT m_dFreq3{ 1.0 }; ///< The frequency scaling of the z-dimension. Used to scale the output.
166 :
167 : realArrayT *m_psdSqrt{ nullptr }; ///< Pointer to the real array containing the square root of the PSD.
168 : bool m_owner{ false }; ///< Flag indicates whether or not m_psdSqrt was allocated by this instance, and so must be
169 : ///< deallocated.
170 :
171 : mutable complexArrayT
172 : m_ftWork; ///< Working memory for the FFT. Declared mutable so it can be accessed in the const filter method.
173 :
174 : math::ft::fftT<complexT, complexT, rank, 0> m_fft_fwd; ///< FFT object for the forward transform.
175 : math::ft::fftT<complexT, complexT, rank, 0> m_fft_bwd; ///< FFT object for the backward transfsorm.
176 :
177 : public:
178 : /// C'tor.
179 : /**
180 : */
181 : psdFilter();
182 :
183 : /// Destructor
184 : ~psdFilter();
185 :
186 : protected:
187 : /// Set the sqaure-root of the PSD to be a pointer to an array containing the square root of the properly normalized
188 : /// PSD.
189 : /** This does not allocate _npsdSqrt, it merely points to the specified array, which remains your responsibility for
190 : * deallocation, etc.
191 : *
192 : * See the discussion of PSD normalization above.
193 : *
194 : * This private version handles the actual setting of m_psdSqrt, which is rank independent.
195 : *
196 : * \returns 0 on success
197 : * \returns -1 on error
198 : *
199 : */
200 : int psdSqrt( realArrayT *npsdSqrt /**< [in] a pointer to an array containing the square root of the PSD. */ );
201 :
202 : /// Set the sqaure-root of the PSD.
203 : /** This allocates _npsdSqrt and fills it with the values in the array.
204 : *
205 : * See the discussion of PSD normalization above.
206 : *
207 : * This private version handles the actual setting of m_psdSqrt, which is rank independent.
208 : *
209 : * \returns 0 on success
210 : * \returns -1 on error
211 : *
212 : */
213 : int psdSqrt( const realArrayT &npsdSqrt /**< [in] an array containing the square root of the PSD. */ );
214 :
215 : /// Set the size of the filter.
216 : /** Handles allocation of the m_ftWork array and fftw planning.
217 : *
218 : * Requires m_psdSqrt to be set first. This is called by the psdSqrt() and psd() methods.
219 : *
220 : * This version compiles when rank==1
221 : *
222 : */
223 : template <size_t crank = rank>
224 : int setSize( typename std::enable_if<crank == 1>::type * = 0 );
225 :
226 : /// Set the size of the filter.
227 : /** Handles allocation of the m_ftWork array and fftw planning.
228 : *
229 : * Requires m_psdSqrt to be set first. This is called by the psdSqrt() and psd() methods.
230 : *
231 : * This version compiles when rank==2
232 : *
233 : */
234 : template <size_t crank = rank>
235 : int setSize( typename std::enable_if<crank == 2>::type * = 0 );
236 :
237 : /// Set the size of the filter.
238 : /** Handles allocation of the m_ftWork array and fftw planning.
239 : *
240 : * Requires m_psdSqrt to be set first. This is called by the psdSqrt() and psd() methods.
241 : *
242 : * This version compiles when rank==3
243 : *
244 : */
245 : template <size_t crank = rank>
246 : int setSize( typename std::enable_if<crank == 3>::type * = 0 );
247 :
248 : public:
249 : /// Get the number of rows in the filter
250 : /**
251 : * \returns the current value of m_rows.
252 : *
253 : */
254 : int rows();
255 :
256 : /// Get the number of columns in the filter
257 : /**
258 : * \returns the current value of m_cols.
259 : *
260 : */
261 : int cols();
262 :
263 : /// Get the number of planes in the filter
264 : /**
265 : * \returns the current value of m_planes.
266 : *
267 : */
268 : int planes();
269 :
270 : /// Set the sqaure-root of the PSD to be a pointer to an array containing the square root of the properly normalized
271 : /// PSD.
272 : /** This does not allocate _npsdSqrt, it merely points to the specified array, which remains your responsibility for
273 : * deallocation, etc.
274 : *
275 : * See the discussion of PSD normalization above.
276 : *
277 : * This version compiles when rank==1
278 : *
279 : * \returns 0 on success
280 : * \returns -1 on error
281 : *
282 : */
283 : template <size_t crank = rank>
284 : int psdSqrt( realArrayT *npsdSqrt, ///< [in] a pointer to an array containing the square root of the PSD.
285 : realT df, ///< [in] the frequency spacing
286 : typename std::enable_if<crank == 1>::type * = 0 );
287 :
288 : /// Set the square-root of the PSD to be a pointer to an array containing the square root of the properly normalized
289 : /// PSD.
290 : /** This does not allocate _npsdSqrt, it merely points to the specified array, which remains your responsibility for
291 : * deallocation, etc.
292 : *
293 : * See the discussion of PSD normalization above.
294 : *
295 : * This version compiles when rank==2
296 : *
297 : * \returns 0 on success
298 : * \returns -1 on error
299 : *
300 : */
301 : template <size_t crank = rank>
302 : int psdSqrt( realArrayT *npsdSqrt, ///< [in] a pointer to an array containing the square root of the PSD.
303 : realT dk1, ///< [in] the frequency spacing along dimension 1
304 : realT dk2, ///< [in] the frequency spacing along dimension 2
305 : typename std::enable_if<crank == 2>::type * = 0 );
306 :
307 : /// Set the sqaure-root of the PSD to be a pointer to an array containing the square root of the properly normalized
308 : /// PSD.
309 : /** This does not allocate _npsdSqrt, it merely points to the specified array, which remains your responsibility for
310 : * deallocation, etc.
311 : *
312 : * See the discussion of PSD normalization above.
313 : *
314 : * This version compiles when rank==3
315 : *
316 : * \returns 0 on success
317 : * \returns -1 on error
318 : *
319 : */
320 : template <size_t crank = rank>
321 : int psdSqrt( realArrayT *npsdSqrt, ///< [in] a pointer to an array containing the square root of the PSD.
322 : realT dk1, ///< [in] the frequency spacing along dimension 1
323 : realT dk2, ///< [in] the frequency spacing along dimension 2
324 : realT df, ///< [in] the frequency spacing along dimension 3
325 : typename std::enable_if<crank == 3>::type * = 0 );
326 :
327 : /// Set the sqaure-root of the PSD.
328 : /** This allocates _npsdSqrt and fills it with th evalues in the array.
329 : *
330 : * See the discussion of PSD normalization above.
331 : *
332 : * This version compiles when rank==1
333 : *
334 : * \returns 0 on success
335 : * \returns -1 on error
336 : *
337 : */
338 : template <size_t crank = rank>
339 : int psdSqrt( const realArrayT &npsdSqrt, ///< [in] an array containing the square root of the PSD.
340 : realT df, ///< [in] the frequency spacing
341 : typename std::enable_if<crank == 1>::type * = 0 );
342 :
343 : /// Set the sqaure-root of the PSD.
344 : /** This allocates _npsdSqrt and fills it with th evalues in the array.
345 : *
346 : * See the discussion of PSD normalization above.
347 : *
348 : * This version compiles when rank==2
349 : *
350 : * \returns 0 on success
351 : * \returns -1 on error
352 : *
353 : */
354 : template <size_t crank = rank>
355 : int psdSqrt( const realArrayT &npsdSqrt, ///< [in] an array containing the square root of the PSD.
356 : realT dk1, ///< [in] the frequency spacing along dimension 1
357 : realT dk2, ///< [in] the frequency spacing along dimension 2
358 : typename std::enable_if<crank == 2>::type * = 0 );
359 :
360 : /// Set the sqaure-root of the PSD.
361 : /** This allocates _npsdSqrt and fills it with th evalues in the array.
362 : *
363 : * See the discussion of PSD normalization above.
364 : *
365 : * This version compiles when rank==3
366 : *
367 : * \returns 0 on success
368 : * \returns -1 on error
369 : *
370 : */
371 : template <size_t crank = rank>
372 : int psdSqrt( const realArrayT &npsdSqrt, ///< [in] an array containing the square root of the PSD.
373 : realT dk1, ///< [in] the frequency spacing along dimension 1
374 : realT dk2, ///< [in] the frequency spacing along dimension 2
375 : realT df, ///< [in] the frequency spacing along dimension 3
376 : typename std::enable_if<crank == 3>::type * = 0 );
377 :
378 : /// Set the sqaure-root of the PSD from the PSD.
379 : /** This allocates _npsdSqrt and fills it with the square root of the values in the array.
380 : *
381 : * See the discussion of PSD normalization above.
382 : *
383 : * This version compiles when rank==1
384 : *
385 : * \returns 0 on success
386 : * \returns -1 on error
387 : *
388 : */
389 : template <size_t crank = rank>
390 : int psd( const realArrayT &npsd, ///< [in] an array containing the PSD
391 : const realT df, ///< [in] the frequency spacing
392 : typename std::enable_if<crank == 1>::type * = 0 );
393 :
394 : /// Set the sqaure-root of the PSD from the PSD.
395 : /** This allocates _npsdSqrt and fills it with the square root of the values in the array.
396 : *
397 : * See the discussion of PSD normalization above.
398 : *
399 : * This version compiles when rank==2
400 : *
401 : * \returns 0 on success
402 : * \returns -1 on error
403 : *
404 : */
405 : template <size_t crank = rank>
406 : int psd( const realArrayT &npsd, ///< [in] an array containing the PSD
407 : const realT dk1, ///< [in] the frequency spacing along dimension 1
408 : const realT dk2, ///< [in] the frequency spacing along dimension 2
409 : typename std::enable_if<crank == 2>::type * = 0 );
410 :
411 : /// Set the sqaure-root of the PSD from the PSD.
412 : /** This allocates _npsdSqrt and fills it with the square root of the values in the array.
413 : *
414 : * See the discussion of PSD normalization above.
415 : *
416 : * This version compiles when rank==3
417 : *
418 : * \returns 0 on success
419 : * \returns -1 on error
420 : *
421 : */
422 : template <size_t crank = rank>
423 : int psd( const realArrayT &npsd, ///< [in] an array containing the PSD
424 : const realT dk1, ///< [in] the frequency spacing along dimension 1
425 : const realT dk2, ///< [in] the frequency spacing along dimension 2
426 : const realT df, ///< [in] the frequency spacing along dimension 3
427 : typename std::enable_if<crank == 3>::type * = 0 );
428 :
429 : /// De-allocate all working memory and reset to initial state.
430 : /**
431 : *
432 : */
433 : void clear();
434 :
435 : /// Apply the filter.
436 : /**
437 : * This version compiles when rank==1
438 : *
439 : * \returns 0 on success
440 : * \returns -1 on error
441 : *
442 : */
443 : template <size_t crank = rank>
444 : int filter( realArrayT &noise, ///< [in.out] the noise field of size rows() X cols(), which is filtered in-place.
445 : realArrayT *noiseIm = nullptr, ///< [out] [optional] an array to fill with the imaginary output of the
446 : ///< filter, allowing 2-for-1 calculation.
447 : typename std::enable_if<crank == 1>::type * = 0 ) const;
448 :
449 : /// Apply the filter.
450 : /**
451 : * This version compiles when rank==2
452 : *
453 : * \returns 0 on success
454 : * \returns -1 on error
455 : *
456 : */
457 : template <size_t crank = rank>
458 : int filter( realArrayT &noise, ///< [in.out] the noise field of size rows() X cols(), which is filtered in-place.
459 : realArrayT *noiseIm = nullptr, ///< [out] [optional] an array to fill with the imaginary output of the
460 : ///< filter, allowing 2-for-1 calculation.
461 : typename std::enable_if<crank == 2>::type * = 0 ) const;
462 :
463 : /// Apply the filter.
464 : /**
465 : * This version compiles when rank==2
466 : *
467 : * \returns 0 on success
468 : * \returns -1 on error
469 : *
470 : */
471 : template <size_t crank = rank>
472 : int filter( realArrayMapT noise, ///< [in.out] the noise field of size rows() X cols(), which is filtered in-place.
473 : realArrayT *noiseIm = nullptr, ///< [out] [optional] an array to fill with the imaginary output of the
474 : ///< filter, allowing 2-for-1 calculation.
475 : typename std::enable_if<crank == 2>::type * = 0 ) const;
476 :
477 : /// Apply the filter.
478 : /**
479 : * This version compiles when rank==3
480 : *
481 : * \returns 0 on success
482 : * \returns -1 on error
483 : *
484 : */
485 : template <size_t crank = rank>
486 : int filter( realArrayT &noise, ///< [in.out] the noise field of size rows() X cols(), which is filtered in-place.
487 : realArrayT *noiseIm = nullptr, ///< [out] [optional] an array to fill with the imaginary output of the
488 : ///< filter, allowing 2-for-1 calculation.
489 : typename std::enable_if<crank == 3>::type * = 0 ) const;
490 :
491 : /// Apply the filter.
492 : /**
493 : * \returns 0 on success
494 : * \returns -1 on error
495 : *
496 : */
497 : int operator()(
498 : realArrayT &noise /**< [in.out] the noise field of size rows() X cols(), which is filtered in-place. */ ) const;
499 :
500 : /// Apply the filter.
501 : /**
502 : * \overload
503 : *
504 : * \returns 0 on success
505 : * \returns -1 on error
506 : *
507 : */
508 : int operator()(
509 : realArrayMapT noise /**< [in.out] the noise field of size rows() X cols(), which is filtered in-place. */ )
510 : const;
511 :
512 : /// Apply the filter.
513 : /**
514 : * \returns 0 on success
515 : * \returns -1 on error
516 : *
517 : */
518 : int
519 : operator()( realArrayT &noise, ///< [in.out] the noise field of size rows() X cols(), which is filtered in-place.
520 : realArrayT &noiseIm ///< [out] [optional] an array to fill with the imaginary output of the filter,
521 : ///< allowing 2-for-1 calculation.
522 : ) const;
523 : };
524 :
525 : template <typename realT, size_t rank>
526 17 : psdFilter<realT, rank>::psdFilter()
527 : {
528 17 : }
529 :
530 : template <typename realT, size_t rank>
531 17 : psdFilter<realT, rank>::~psdFilter()
532 : {
533 17 : if( m_psdSqrt && m_owner )
534 : {
535 7 : delete m_psdSqrt;
536 : }
537 17 : }
538 :
539 : template <typename realT, size_t rank>
540 3 : int psdFilter<realT, rank>::psdSqrt( realArrayT *npsdSqrt )
541 : {
542 3 : if( m_psdSqrt && m_owner )
543 : {
544 0 : delete m_psdSqrt;
545 : }
546 :
547 3 : m_psdSqrt = npsdSqrt;
548 3 : m_owner = false;
549 :
550 3 : setSize();
551 :
552 3 : return 0;
553 : }
554 :
555 : template <typename realT, size_t rank>
556 4 : int psdFilter<realT, rank>::psdSqrt( const realArrayT &npsdSqrt )
557 : {
558 4 : if( m_psdSqrt && m_owner )
559 : {
560 0 : delete m_psdSqrt;
561 : }
562 :
563 4 : m_psdSqrt = new realArrayT;
564 :
565 4 : ( *m_psdSqrt ) = npsdSqrt;
566 4 : m_owner = true;
567 :
568 4 : setSize();
569 :
570 4 : return 0;
571 : }
572 :
573 : template <typename realT, size_t rank>
574 : template <size_t crank>
575 5 : int psdFilter<realT, rank>::setSize( typename std::enable_if<crank == 1>::type * )
576 : {
577 5 : if( m_psdSqrt == 0 )
578 : {
579 0 : internal::mxlib_error_report(error_t::paramnotset, "m_psdSqrt has not been set yet, is still NULL." );
580 0 : return -1;
581 : }
582 :
583 5 : if( m_rows == m_psdSqrt->size() )
584 : {
585 0 : return 0;
586 : }
587 :
588 5 : m_rows = m_psdSqrt->size();
589 5 : m_cols = 1;
590 5 : m_planes = 1;
591 :
592 5 : m_ftWork.resize( m_rows );
593 :
594 5 : m_fft_fwd.plan( m_rows, math::ft::dir::forward, true );
595 :
596 5 : m_fft_bwd.plan( m_rows, math::ft::dir::backward, true );
597 :
598 5 : return 0;
599 : }
600 :
601 : template <typename realT, size_t rank>
602 : template <size_t crank>
603 6 : int psdFilter<realT, rank>::setSize( typename std::enable_if<crank == 2>::type * )
604 : {
605 6 : if( m_psdSqrt == 0 )
606 : {
607 0 : internal::mxlib_error_report(error_t::paramnotset, "m_psdSqrt has not been set yet, is still NULL." );
608 0 : return -1;
609 : }
610 :
611 6 : if( m_rows == m_psdSqrt->rows() && m_cols == m_psdSqrt->cols() )
612 : {
613 0 : return 0;
614 : }
615 :
616 6 : m_rows = m_psdSqrt->rows();
617 6 : m_cols = m_psdSqrt->cols();
618 6 : m_planes = 1;
619 :
620 6 : m_ftWork.resize( m_rows, m_cols );
621 :
622 6 : m_fft_fwd.plan( m_rows, m_cols, math::ft::dir::forward, true );
623 :
624 6 : m_fft_bwd.plan( m_rows, m_cols, math::ft::dir::backward, true );
625 :
626 6 : return 0;
627 : }
628 :
629 : template <typename realT, size_t rank>
630 : template <size_t crank>
631 5 : int psdFilter<realT, rank>::setSize( typename std::enable_if<crank == 3>::type * )
632 : {
633 5 : if( m_psdSqrt == 0 )
634 : {
635 0 : internal::mxlib_error_report(error_t::paramnotset, "m_psdSqrt has not been set yet, is still NULL." );
636 0 : return -1;
637 : }
638 :
639 5 : if( m_rows == m_psdSqrt->rows() && m_cols == m_psdSqrt->cols() && m_planes == m_psdSqrt->planes() )
640 : {
641 0 : return 0;
642 : }
643 :
644 5 : m_rows = m_psdSqrt->rows();
645 5 : m_cols = m_psdSqrt->cols();
646 5 : m_planes = m_psdSqrt->planes();
647 :
648 5 : m_ftWork.resize( m_rows, m_cols, m_planes );
649 :
650 5 : m_fft_fwd.plan( m_planes, m_rows, m_cols, math::ft::dir::forward, true );
651 :
652 5 : m_fft_bwd.plan( m_planes, m_rows, m_cols, math::ft::dir::backward, true );
653 :
654 5 : return 0;
655 : }
656 :
657 : template <typename realT, size_t rank>
658 30 : int psdFilter<realT, rank>::rows()
659 : {
660 30 : return m_rows;
661 : }
662 :
663 : template <typename realT, size_t rank>
664 28 : int psdFilter<realT, rank>::cols()
665 : {
666 28 : return m_cols;
667 : }
668 :
669 : template <typename realT, size_t rank>
670 26 : int psdFilter<realT, rank>::planes()
671 : {
672 26 : return m_planes;
673 : }
674 :
675 : template <typename realT, size_t rank>
676 : template <size_t crank>
677 1 : int psdFilter<realT, rank>::psdSqrt( realArrayT *npsdSqrt, realT df, typename std::enable_if<crank == 1>::type * )
678 : {
679 1 : m_dFreq1 = df;
680 1 : return psdSqrt( npsdSqrt );
681 : }
682 :
683 : template <typename realT, size_t rank>
684 : template <size_t crank>
685 1 : int psdFilter<realT, rank>::psdSqrt( realArrayT *npsdSqrt,
686 : realT dk1,
687 : realT dk2,
688 : typename std::enable_if<crank == 2>::type * )
689 : {
690 1 : m_dFreq1 = dk1;
691 1 : m_dFreq2 = dk2;
692 1 : return psdSqrt( npsdSqrt );
693 : }
694 :
695 : template <typename realT, size_t rank>
696 : template <size_t crank>
697 1 : int psdFilter<realT, rank>::psdSqrt(
698 : realArrayT *npsdSqrt, realT dk1, realT dk2, realT df, typename std::enable_if<crank == 3>::type * )
699 : {
700 1 : m_dFreq1 = dk1;
701 1 : m_dFreq2 = dk2;
702 1 : m_dFreq3 = df;
703 1 : return psdSqrt( npsdSqrt );
704 : }
705 :
706 : template <typename realT, size_t rank>
707 : template <size_t crank>
708 1 : int psdFilter<realT, rank>::psdSqrt( const realArrayT &npsdSqrt, realT df, typename std::enable_if<crank == 1>::type * )
709 : {
710 1 : m_dFreq1 = df;
711 1 : return psdSqrt( npsdSqrt );
712 : }
713 :
714 : template <typename realT, size_t rank>
715 : template <size_t crank>
716 2 : int psdFilter<realT, rank>::psdSqrt( const realArrayT &npsdSqrt,
717 : realT dk1,
718 : realT dk2,
719 : typename std::enable_if<crank == 2>::type * )
720 : {
721 2 : m_dFreq1 = dk1;
722 2 : m_dFreq2 = dk2;
723 2 : return psdSqrt( npsdSqrt );
724 : }
725 :
726 : template <typename realT, size_t rank>
727 : template <size_t crank>
728 1 : int psdFilter<realT, rank>::psdSqrt(
729 : const realArrayT &npsdSqrt, realT dk1, realT dk2, realT df, typename std::enable_if<crank == 3>::type * )
730 : {
731 1 : m_dFreq1 = dk1;
732 1 : m_dFreq2 = dk2;
733 1 : m_dFreq3 = df;
734 1 : return psdSqrt( npsdSqrt );
735 : }
736 :
737 : template <typename realT, size_t rank>
738 : template <size_t crank>
739 3 : int psdFilter<realT, rank>::psd( const realArrayT &npsd, const realT df1, typename std::enable_if<crank == 1>::type * )
740 : {
741 3 : if( m_psdSqrt && m_owner )
742 : {
743 0 : delete m_psdSqrt;
744 : }
745 :
746 3 : m_psdSqrt = new realArrayT;
747 :
748 : // Vector
749 3 : m_psdSqrt->resize( npsd.size() );
750 7171 : for( size_t n = 0; n < npsd.size(); ++n )
751 7168 : ( *m_psdSqrt )[n] = sqrt( npsd[n] );
752 :
753 3 : m_owner = true;
754 :
755 3 : m_dFreq1 = df1;
756 :
757 3 : setSize();
758 :
759 3 : return 0;
760 : }
761 :
762 : template <typename realT, size_t rank>
763 : template <size_t crank>
764 3 : int psdFilter<realT, rank>::psd( const realArrayT &npsd,
765 : const realT dk1,
766 : const realT dk2,
767 : typename std::enable_if<crank == 2>::type * )
768 : {
769 3 : if( m_psdSqrt && m_owner )
770 : {
771 0 : delete m_psdSqrt;
772 : }
773 :
774 3 : m_psdSqrt = new realArrayT;
775 :
776 3 : ( *m_psdSqrt ) = npsd.sqrt();
777 3 : m_owner = true;
778 :
779 3 : m_dFreq1 = dk1;
780 3 : m_dFreq2 = dk2;
781 :
782 3 : setSize();
783 :
784 3 : return 0;
785 : }
786 :
787 : template <typename realT, size_t rank>
788 : template <size_t crank>
789 3 : int psdFilter<realT, rank>::psd( const realArrayT &npsd,
790 : const realT dk1,
791 : const realT dk2,
792 : const realT df,
793 : typename std::enable_if<crank == 3>::type * )
794 : {
795 3 : if( m_psdSqrt && m_owner )
796 : {
797 0 : delete m_psdSqrt;
798 : }
799 :
800 3 : m_psdSqrt = new realArrayT;
801 :
802 : // Cube
803 3 : m_psdSqrt->resize( npsd.rows(), npsd.cols(), npsd.planes() );
804 387 : for( int pp = 0; pp < npsd.planes(); ++pp )
805 : {
806 37248 : for( int cc = 0; cc < npsd.cols(); ++cc )
807 : {
808 4362240 : for( int rr = 0; rr < npsd.rows(); ++rr )
809 : {
810 4325376 : m_psdSqrt->image( pp )( rr, cc ) = sqrt( npsd.image( pp )( rr, cc ) );
811 : }
812 : }
813 : }
814 :
815 3 : m_owner = true;
816 :
817 3 : m_dFreq1 = dk1;
818 3 : m_dFreq2 = dk2;
819 3 : m_dFreq3 = df;
820 :
821 3 : setSize();
822 :
823 3 : return 0;
824 : }
825 :
826 : template <typename realT, size_t rank>
827 9 : void psdFilter<realT, rank>::clear()
828 : {
829 : // m_ftWork.resize(0,0);
830 9 : psdFilterTypes::arrayT<realT, rank>::clear( m_ftWork );
831 :
832 9 : m_rows = 0;
833 9 : m_cols = 0;
834 9 : m_planes = 0;
835 :
836 9 : if( m_psdSqrt && m_owner )
837 : {
838 6 : delete m_psdSqrt;
839 6 : m_psdSqrt = 0;
840 : }
841 9 : }
842 :
843 : template <typename realT, size_t rank>
844 : template <size_t crank>
845 200 : int psdFilter<realT, rank>::filter( realArrayT &noise,
846 : realArrayT *noiseIm,
847 : typename std::enable_if<crank == 1>::type * ) const
848 : {
849 614600 : for( int nn = 0; nn < noise.size(); ++nn )
850 614400 : m_ftWork[nn] = complexT( noise[nn], 0 );
851 :
852 : // Transform complex noise to Fourier domain.
853 200 : m_fft_fwd( m_ftWork.data(), m_ftWork.data() );
854 :
855 : // Apply the filter.
856 614600 : for( int nn = 0; nn < m_ftWork.size(); ++nn )
857 614400 : m_ftWork[nn] *= ( *m_psdSqrt )[nn];
858 :
859 200 : m_fft_bwd( m_ftWork.data(), m_ftWork.data() );
860 :
861 : // Now take the real part, and normalize.
862 200 : realT norm = sqrt( noise.size() / m_dFreq1 );
863 614600 : for( int nn = 0; nn < m_ftWork.size(); ++nn )
864 614400 : noise[nn] = m_ftWork[nn].real() / norm;
865 :
866 200 : if( noiseIm != nullptr )
867 : {
868 0 : for( int nn = 0; nn < m_ftWork.size(); ++nn )
869 0 : ( *noiseIm )[nn] = m_ftWork[nn].imag() / norm;
870 : }
871 :
872 200 : return 0;
873 : }
874 :
875 : template <typename realT, size_t rank>
876 : template <size_t crank>
877 201 : int psdFilter<realT, rank>::filter( realArrayT &noise,
878 : realArrayT *noiseIm,
879 : typename std::enable_if<crank == 2>::type * ) const
880 : {
881 : // Make noise a complex number
882 13193 : for( int ii = 0; ii < noise.rows(); ++ii )
883 : {
884 869056 : for( int jj = 0; jj < noise.cols(); ++jj )
885 : {
886 856064 : m_ftWork( ii, jj ) = complexT( noise( ii, jj ), 0 );
887 : }
888 : }
889 :
890 : // Transform complex noise to Fourier domain.
891 201 : m_fft_fwd( m_ftWork.data(), m_ftWork.data() );
892 :
893 : // Apply the filter.
894 201 : m_ftWork *= *m_psdSqrt;
895 :
896 201 : m_fft_bwd( m_ftWork.data(), m_ftWork.data() );
897 :
898 201 : realT norm = sqrt( noise.rows() * noise.cols() / ( m_dFreq1 * m_dFreq2 ) );
899 :
900 : // Now take the real part, and normalize.
901 201 : noise = m_ftWork.real() / norm;
902 :
903 201 : if( noiseIm != nullptr )
904 : {
905 0 : *noiseIm = m_ftWork.imag() / norm;
906 : }
907 :
908 201 : return 0;
909 : }
910 :
911 : template <typename realT, size_t rank>
912 : template <size_t crank>
913 : int psdFilter<realT, rank>::filter( realArrayMapT noise,
914 : realArrayT *noiseIm,
915 : typename std::enable_if<crank == 2>::type * ) const
916 : {
917 : // Make noise a complex number
918 : for( int ii = 0; ii < noise.rows(); ++ii )
919 : {
920 : for( int jj = 0; jj < noise.cols(); ++jj )
921 : {
922 : m_ftWork( ii, jj ) = complexT( noise( ii, jj ), 0 );
923 : }
924 : }
925 :
926 : // Transform complex noise to Fourier domain.
927 : m_fft_fwd( m_ftWork.data(), m_ftWork.data() );
928 :
929 : // Apply the filter.
930 : m_ftWork *= *m_psdSqrt;
931 :
932 : m_fft_bwd( m_ftWork.data(), m_ftWork.data() );
933 :
934 : realT norm = sqrt( noise.rows() * noise.cols() / ( m_dFreq1 * m_dFreq2 ) );
935 :
936 : // Now take the real part, and normalize.
937 : noise = m_ftWork.real() / norm;
938 :
939 : if( noiseIm != nullptr )
940 : {
941 : *noiseIm = m_ftWork.imag() / norm;
942 : }
943 :
944 : return 0;
945 : }
946 :
947 : template <typename realT, size_t rank>
948 : template <size_t crank>
949 200 : int psdFilter<realT, rank>::filter( realArrayT &noise,
950 : realArrayT *noiseIm,
951 : typename std::enable_if<crank == 3>::type * ) const
952 : {
953 : // Make noise a complex number
954 13000 : for( int pp = 0; pp < noise.planes(); ++pp )
955 : {
956 422400 : for( int ii = 0; ii < noise.rows(); ++ii )
957 : {
958 13516800 : for( int jj = 0; jj < noise.cols(); ++jj )
959 : {
960 13107200 : m_ftWork.image( pp )( ii, jj ) = complexT( noise.image( pp )( ii, jj ), 0 );
961 : }
962 : }
963 : }
964 :
965 : // Transform complex noise to Fourier domain.
966 200 : m_fft_fwd( m_ftWork.data(), m_ftWork.data() );
967 :
968 : // Apply the filter.
969 13000 : for( int pp = 0; pp < noise.planes(); ++pp )
970 12800 : m_ftWork.image( pp ) *= m_psdSqrt->image( pp );
971 :
972 200 : m_fft_bwd( m_ftWork.data(), m_ftWork.data() );
973 :
974 : // Now take the real part, and normalize.
975 :
976 200 : realT norm = sqrt( m_rows * m_cols * m_planes / ( m_dFreq1 * m_dFreq2 * m_dFreq3 ) );
977 13000 : for( int pp = 0; pp < noise.planes(); ++pp )
978 12800 : noise.image( pp ) = m_ftWork.image( pp ).real() / norm;
979 :
980 200 : if( noiseIm != nullptr )
981 : {
982 0 : for( int pp = 0; pp < noise.planes(); ++pp )
983 0 : noiseIm->image( pp ) = m_ftWork.image( pp ).imag() / norm;
984 : }
985 :
986 200 : return 0;
987 : }
988 :
989 : template <typename realT, size_t rank>
990 601 : int psdFilter<realT, rank>::operator()( realArrayT &noise ) const
991 : {
992 601 : return filter( noise );
993 : }
994 :
995 : template <typename realT, size_t rank>
996 : int psdFilter<realT, rank>::operator()( realArrayMapT noise ) const
997 : {
998 : return filter( noise );
999 : }
1000 :
1001 : template <typename realT, size_t rank>
1002 : int psdFilter<realT, rank>::operator()( realArrayT &noise, realArrayT &noiseIm ) const
1003 : {
1004 : return filter( noise, &noiseIm );
1005 : }
1006 :
1007 : } // namespace sigproc
1008 : } // namespace mx
1009 :
1010 : #endif // psdFilter_hpp
|
