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 : * \test Scenario: compiling psdFilter. \ref tests_sigproc_psdFilter_compile "[test doc]"
142 : */
143 : template <typename _realT, size_t _rank>
144 : class psdFilter<_realT, _rank, 0>
145 : {
146 : public:
147 : typedef _realT realT; ///< Real floating point type
148 : typedef std::complex<_realT> complexT; ///< Complex floating point type.
149 :
150 : static const size_t rank = _rank;
151 :
152 : typedef typename psdFilterTypes::arrayT<realT, rank>::realArrayT
153 : realArrayT; ///< std::vector for rank==1, Eigen::Array for rank==2, eigenCube for rank==3.
154 : typedef typename psdFilterTypes::arrayT<realT, rank>::realArrayMapT
155 : realArrayMapT; ///< std::vector for rank==1, Eigen::Map for rank==2, eigenCube for rank==3.
156 : typedef typename psdFilterTypes::arrayT<realT, rank>::complexArrayT
157 : complexArrayT; ///< std::vector for rank==1, Eigen::Array for rank==2, eigenCube for rank==3.
158 :
159 : protected:
160 : int m_rows{ 0 }; ///< The number of rows in the filter, and the required number of rows in the noise array.
161 : int m_cols{ 0 }; ///< The number of columns in the filter, and the required number of columns in the noise array.
162 : int m_planes{ 0 }; ///< Then number of planes in the filter.
163 :
164 : realT m_dFreq1{ 1.0 }; ///< The frequency scaling of the x-dimension. Used to scale the output.
165 : realT m_dFreq2{ 1.0 }; ///< The frequency scaling of the y-dimension. Used to scale the output.
166 : realT m_dFreq3{ 1.0 }; ///< The frequency scaling of the z-dimension. Used to scale the output.
167 :
168 : realArrayT *m_psdSqrt{ nullptr }; ///< Pointer to the real array containing the square root of the PSD.
169 : bool m_owner{ false }; ///< Flag indicates whether or not m_psdSqrt was allocated by this instance, and so must be
170 : ///< deallocated.
171 :
172 : mutable complexArrayT
173 : m_ftWork; ///< Working memory for the FFT. Declared mutable so it can be accessed in the const filter method.
174 :
175 : math::ft::fftT<complexT, complexT, rank, 0> m_fft_fwd; ///< FFT object for the forward transform.
176 : math::ft::fftT<complexT, complexT, rank, 0> m_fft_bwd; ///< FFT object for the backward transfsorm.
177 :
178 : public:
179 : /// C'tor.
180 : /**
181 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
182 : * "[test doc]"
183 : */
184 : psdFilter();
185 :
186 : /// Destructor
187 : ~psdFilter();
188 :
189 : protected:
190 : /// Set the sqaure-root of the PSD to be a pointer to an array containing the square root of the properly normalized
191 : /// PSD.
192 : /** This does not allocate _npsdSqrt, it merely points to the specified array, which remains your responsibility for
193 : * deallocation, etc.
194 : *
195 : * See the discussion of PSD normalization above.
196 : *
197 : * This private version handles the actual setting of m_psdSqrt, which is rank independent.
198 : *
199 : * \returns 0 on success
200 : * \returns -1 on error
201 : *
202 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
203 : * "[test doc]"
204 : */
205 : int psdSqrt( realArrayT *npsdSqrt /**< [in] a pointer to an array containing the square root of the PSD. */ );
206 :
207 : /// Set the sqaure-root of the PSD.
208 : /** This allocates _npsdSqrt and fills it with the values in the array.
209 : *
210 : * See the discussion of PSD normalization above.
211 : *
212 : * This private version handles the actual setting of m_psdSqrt, which is rank independent.
213 : *
214 : * \returns 0 on success
215 : * \returns -1 on error
216 : *
217 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
218 : * "[test doc]"
219 : */
220 : int psdSqrt( const realArrayT &npsdSqrt /**< [in] an array containing the square root of the PSD. */ );
221 :
222 : /// Set the size of the filter.
223 : /** Handles allocation of the m_ftWork array and fftw planning.
224 : *
225 : * Requires m_psdSqrt to be set first. This is called by the psdSqrt() and psd() methods.
226 : *
227 : * This version compiles when rank==1
228 : *
229 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
230 : * "[test doc]"
231 : */
232 : template <size_t crank = rank>
233 : int setSize( typename std::enable_if<crank == 1>::type * = 0 );
234 :
235 : /// Set the size of the filter.
236 : /** Handles allocation of the m_ftWork array and fftw planning.
237 : *
238 : * Requires m_psdSqrt to be set first. This is called by the psdSqrt() and psd() methods.
239 : *
240 : * This version compiles when rank==2
241 : *
242 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
243 : * "[test doc]"
244 : */
245 : template <size_t crank = rank>
246 : int setSize( typename std::enable_if<crank == 2>::type * = 0 );
247 :
248 : /// Set the size of the filter.
249 : /** Handles allocation of the m_ftWork array and fftw planning.
250 : *
251 : * Requires m_psdSqrt to be set first. This is called by the psdSqrt() and psd() methods.
252 : *
253 : * This version compiles when rank==3
254 : *
255 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
256 : * "[test doc]"
257 : */
258 : template <size_t crank = rank>
259 : int setSize( typename std::enable_if<crank == 3>::type * = 0 );
260 :
261 : public:
262 : /// Get the number of rows in the filter
263 : /**
264 : * \returns the current value of m_rows.
265 : *
266 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
267 : * "[test doc]"
268 : */
269 : int rows();
270 :
271 : /// Get the number of columns in the filter
272 : /**
273 : * \returns the current value of m_cols.
274 : *
275 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
276 : * "[test doc]"
277 : */
278 : int cols();
279 :
280 : /// Get the number of planes in the filter
281 : /**
282 : * \returns the current value of m_planes.
283 : *
284 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
285 : * "[test doc]"
286 : */
287 : int planes();
288 :
289 : /// Set the sqaure-root of the PSD to be a pointer to an array containing the square root of the properly normalized
290 : /// PSD.
291 : /** This does not allocate _npsdSqrt, it merely points to the specified array, which remains your responsibility for
292 : * deallocation, etc.
293 : *
294 : * See the discussion of PSD normalization above.
295 : *
296 : * This version compiles when rank==1
297 : *
298 : * \returns 0 on success
299 : * \returns -1 on error
300 : *
301 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
302 : * "[test doc]"
303 : */
304 : template <size_t crank = rank>
305 : int psdSqrt( realArrayT *npsdSqrt, ///< [in] a pointer to an array containing the square root of the PSD.
306 : realT df, ///< [in] the frequency spacing
307 : typename std::enable_if<crank == 1>::type * = 0 );
308 :
309 : /// Set the square-root of the PSD to be a pointer to an array containing the square root of the properly normalized
310 : /// PSD.
311 : /** This does not allocate _npsdSqrt, it merely points to the specified array, which remains your responsibility for
312 : * deallocation, etc.
313 : *
314 : * See the discussion of PSD normalization above.
315 : *
316 : * This version compiles when rank==2
317 : *
318 : * \returns 0 on success
319 : * \returns -1 on error
320 : *
321 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
322 : * "[test doc]"
323 : */
324 : template <size_t crank = rank>
325 : int psdSqrt( realArrayT *npsdSqrt, ///< [in] a pointer to an array containing the square root of the PSD.
326 : realT dk1, ///< [in] the frequency spacing along dimension 1
327 : realT dk2, ///< [in] the frequency spacing along dimension 2
328 : typename std::enable_if<crank == 2>::type * = 0 );
329 :
330 : /// Set the sqaure-root of the PSD to be a pointer to an array containing the square root of the properly normalized
331 : /// PSD.
332 : /** This does not allocate _npsdSqrt, it merely points to the specified array, which remains your responsibility for
333 : * deallocation, etc.
334 : *
335 : * See the discussion of PSD normalization above.
336 : *
337 : * This version compiles when rank==3
338 : *
339 : * \returns 0 on success
340 : * \returns -1 on error
341 : *
342 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
343 : * "[test doc]"
344 : */
345 : template <size_t crank = rank>
346 : int psdSqrt( realArrayT *npsdSqrt, ///< [in] a pointer to an array containing the square root of the PSD.
347 : realT dk1, ///< [in] the frequency spacing along dimension 1
348 : realT dk2, ///< [in] the frequency spacing along dimension 2
349 : realT df, ///< [in] the frequency spacing along dimension 3
350 : typename std::enable_if<crank == 3>::type * = 0 );
351 :
352 : /// Set the sqaure-root of the PSD.
353 : /** This allocates _npsdSqrt and fills it with th evalues in the array.
354 : *
355 : * See the discussion of PSD normalization above.
356 : *
357 : * This version compiles when rank==1
358 : *
359 : * \returns 0 on success
360 : * \returns -1 on error
361 : *
362 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
363 : * "[test doc]"
364 : */
365 : template <size_t crank = rank>
366 : int psdSqrt( const realArrayT &npsdSqrt, ///< [in] an array containing the square root of the PSD.
367 : realT df, ///< [in] the frequency spacing
368 : typename std::enable_if<crank == 1>::type * = 0 );
369 :
370 : /// Set the sqaure-root of the PSD.
371 : /** This allocates _npsdSqrt and fills it with th evalues in the array.
372 : *
373 : * See the discussion of PSD normalization above.
374 : *
375 : * This version compiles when rank==2
376 : *
377 : * \returns 0 on success
378 : * \returns -1 on error
379 : *
380 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
381 : * "[test doc]"
382 : */
383 : template <size_t crank = rank>
384 : int psdSqrt( const realArrayT &npsdSqrt, ///< [in] an array containing the square root of the PSD.
385 : realT dk1, ///< [in] the frequency spacing along dimension 1
386 : realT dk2, ///< [in] the frequency spacing along dimension 2
387 : typename std::enable_if<crank == 2>::type * = 0 );
388 :
389 : /// Set the sqaure-root of the PSD.
390 : /** This allocates _npsdSqrt and fills it with th evalues in the array.
391 : *
392 : * See the discussion of PSD normalization above.
393 : *
394 : * This version compiles when rank==3
395 : *
396 : * \returns 0 on success
397 : * \returns -1 on error
398 : *
399 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
400 : * "[test doc]"
401 : */
402 : template <size_t crank = rank>
403 : int psdSqrt( const realArrayT &npsdSqrt, ///< [in] an array containing the square root of the PSD.
404 : realT dk1, ///< [in] the frequency spacing along dimension 1
405 : realT dk2, ///< [in] the frequency spacing along dimension 2
406 : realT df, ///< [in] the frequency spacing along dimension 3
407 : typename std::enable_if<crank == 3>::type * = 0 );
408 :
409 : /// Set the sqaure-root of the PSD from the PSD.
410 : /** This allocates _npsdSqrt and fills it with the square root of the values in the array.
411 : *
412 : * See the discussion of PSD normalization above.
413 : *
414 : * This version compiles when rank==1
415 : *
416 : * \returns 0 on success
417 : * \returns -1 on error
418 : *
419 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
420 : * "[test doc]"
421 : */
422 : template <size_t crank = rank>
423 : int psd( const realArrayT &npsd, ///< [in] an array containing the PSD
424 : const realT df, ///< [in] the frequency spacing
425 : typename std::enable_if<crank == 1>::type * = 0 );
426 :
427 : /// Set the sqaure-root of the PSD from the PSD.
428 : /** This allocates _npsdSqrt and fills it with the square root of the values in the array.
429 : *
430 : * See the discussion of PSD normalization above.
431 : *
432 : * This version compiles when rank==2
433 : *
434 : * \returns 0 on success
435 : * \returns -1 on error
436 : *
437 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
438 : * "[test doc]"
439 : */
440 : template <size_t crank = rank>
441 : int psd( const realArrayT &npsd, ///< [in] an array containing the PSD
442 : const realT dk1, ///< [in] the frequency spacing along dimension 1
443 : const realT dk2, ///< [in] the frequency spacing along dimension 2
444 : typename std::enable_if<crank == 2>::type * = 0 );
445 :
446 : /// Set the sqaure-root of the PSD from the PSD.
447 : /** This allocates _npsdSqrt and fills it with the square root of the values in the array.
448 : *
449 : * See the discussion of PSD normalization above.
450 : *
451 : * This version compiles when rank==3
452 : *
453 : * \returns 0 on success
454 : * \returns -1 on error
455 : *
456 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
457 : * "[test doc]"
458 : */
459 : template <size_t crank = rank>
460 : int psd( const realArrayT &npsd, ///< [in] an array containing the PSD
461 : const realT dk1, ///< [in] the frequency spacing along dimension 1
462 : const realT dk2, ///< [in] the frequency spacing along dimension 2
463 : const realT df, ///< [in] the frequency spacing along dimension 3
464 : typename std::enable_if<crank == 3>::type * = 0 );
465 :
466 : /// De-allocate all working memory and reset to initial state.
467 : /**
468 : *
469 : * \test Verify compilation and initialization of the 3 ranks for psdFilter. \ref tests_sigproc_psdFilter_compile
470 : * "[test doc]"
471 : */
472 : void clear();
473 :
474 : /// Apply the filter.
475 : /**
476 : * This version compiles when rank==1
477 : *
478 : * \returns 0 on success
479 : * \returns -1 on error
480 : *
481 : * \test Verify filtering and noise normalization. \ref tests_sigproc_psdFilter_filter "[test doc]"
482 : */
483 : template <size_t crank = rank>
484 : int filter( realArrayT &noise, ///< [in.out] the noise field of size rows() X cols(), which is filtered in-place.
485 : realArrayT *noiseIm = nullptr, ///< [out] [optional] an array to fill with the imaginary output of the
486 : ///< filter, allowing 2-for-1 calculation.
487 : typename std::enable_if<crank == 1>::type * = 0 ) const;
488 :
489 : /// Apply the filter.
490 : /**
491 : * This version compiles when rank==2
492 : *
493 : * \returns 0 on success
494 : * \returns -1 on error
495 : *
496 : * \test Verify filtering and noise normalization. \ref tests_sigproc_psdFilter_filter "[test doc]"
497 : */
498 : template <size_t crank = rank>
499 : int filter( realArrayT &noise, ///< [in.out] the noise field of size rows() X cols(), which is filtered in-place.
500 : realArrayT *noiseIm = nullptr, ///< [out] [optional] an array to fill with the imaginary output of the
501 : ///< filter, allowing 2-for-1 calculation.
502 : typename std::enable_if<crank == 2>::type * = 0 ) const;
503 :
504 : /// Apply the filter.
505 : /**
506 : * This version compiles when rank==2
507 : *
508 : * \returns 0 on success
509 : * \returns -1 on error
510 : *
511 : */
512 : template <size_t crank = rank>
513 : int filter( realArrayMapT noise, ///< [in.out] the noise field of size rows() X cols(), which is filtered in-place.
514 : realArrayT *noiseIm = nullptr, ///< [out] [optional] an array to fill with the imaginary output of the
515 : ///< filter, allowing 2-for-1 calculation.
516 : typename std::enable_if<crank == 2>::type * = 0 ) const;
517 :
518 : /// Apply the filter.
519 : /**
520 : * This version compiles when rank==3
521 : *
522 : * \returns 0 on success
523 : * \returns -1 on error
524 : *
525 : * \test Verify filtering and noise normalization. \ref tests_sigproc_psdFilter_filter "[test doc]"
526 : */
527 : template <size_t crank = rank>
528 : int filter( realArrayT &noise, ///< [in.out] the noise field of size rows() X cols(), which is filtered in-place.
529 : realArrayT *noiseIm = nullptr, ///< [out] [optional] an array to fill with the imaginary output of the
530 : ///< filter, allowing 2-for-1 calculation.
531 : typename std::enable_if<crank == 3>::type * = 0 ) const;
532 :
533 : /// Apply the filter.
534 : /**
535 : * \returns 0 on success
536 : * \returns -1 on error
537 : *
538 : * \test Verify filtering and noise normalization. \ref tests_sigproc_psdFilter_filter "[test doc]"
539 : */
540 : int operator()(
541 : realArrayT &noise /**< [in.out] the noise field of size rows() X cols(), which is filtered in-place. */ ) const;
542 :
543 : /// Apply the filter.
544 : /**
545 : * \overload
546 : *
547 : * \returns 0 on success
548 : * \returns -1 on error
549 : *
550 : * \test Verify filtering and noise normalization. \ref tests_sigproc_psdFilter_filter "[test doc]"
551 : */
552 : int operator()(
553 : realArrayMapT noise /**< [in.out] the noise field of size rows() X cols(), which is filtered in-place. */ )
554 : const;
555 :
556 : /// Apply the filter.
557 : /**
558 : * \returns 0 on success
559 : * \returns -1 on error
560 : *
561 : * \test Verify filtering and noise normalization. \ref tests_sigproc_psdFilter_filter "[test doc]"
562 : */
563 : int
564 : operator()( realArrayT &noise, ///< [in.out] the noise field of size rows() X cols(), which is filtered in-place.
565 : realArrayT &noiseIm ///< [out] [optional] an array to fill with the imaginary output of the filter,
566 : ///< allowing 2-for-1 calculation.
567 : ) const;
568 : };
569 :
570 : template <typename realT, size_t rank>
571 15 : psdFilter<realT, rank>::psdFilter()
572 : {
573 15 : }
574 :
575 : template <typename realT, size_t rank>
576 15 : psdFilter<realT, rank>::~psdFilter()
577 : {
578 15 : if( m_psdSqrt && m_owner )
579 : {
580 6 : delete m_psdSqrt;
581 : }
582 15 : }
583 :
584 : template <typename realT, size_t rank>
585 3 : int psdFilter<realT, rank>::psdSqrt( realArrayT *npsdSqrt )
586 : {
587 3 : if( m_psdSqrt && m_owner )
588 : {
589 0 : delete m_psdSqrt;
590 : }
591 :
592 3 : m_psdSqrt = npsdSqrt;
593 3 : m_owner = false;
594 :
595 3 : setSize();
596 :
597 3 : return 0;
598 : }
599 :
600 : template <typename realT, size_t rank>
601 3 : int psdFilter<realT, rank>::psdSqrt( const realArrayT &npsdSqrt )
602 : {
603 3 : if( m_psdSqrt && m_owner )
604 : {
605 0 : delete m_psdSqrt;
606 : }
607 :
608 3 : m_psdSqrt = new realArrayT;
609 :
610 3 : ( *m_psdSqrt ) = npsdSqrt;
611 3 : m_owner = true;
612 :
613 3 : setSize();
614 :
615 3 : return 0;
616 : }
617 :
618 : template <typename realT, size_t rank>
619 : template <size_t crank>
620 5 : int psdFilter<realT, rank>::setSize( typename std::enable_if<crank == 1>::type * )
621 : {
622 5 : if( m_psdSqrt == 0 )
623 : {
624 0 : internal::mxlib_error_report(error_t::paramnotset, "m_psdSqrt has not been set yet, is still NULL." );
625 0 : return -1;
626 : }
627 :
628 5 : if( m_rows == m_psdSqrt->size() )
629 : {
630 0 : return 0;
631 : }
632 :
633 5 : m_rows = m_psdSqrt->size();
634 5 : m_cols = 1;
635 5 : m_planes = 1;
636 :
637 5 : m_ftWork.resize( m_rows );
638 :
639 5 : m_fft_fwd.plan( m_rows, math::ft::dir::forward, true );
640 :
641 5 : m_fft_bwd.plan( m_rows, math::ft::dir::backward, true );
642 :
643 5 : return 0;
644 : }
645 :
646 : template <typename realT, size_t rank>
647 : template <size_t crank>
648 5 : int psdFilter<realT, rank>::setSize( typename std::enable_if<crank == 2>::type * )
649 : {
650 5 : if( m_psdSqrt == 0 )
651 : {
652 0 : internal::mxlib_error_report(error_t::paramnotset, "m_psdSqrt has not been set yet, is still NULL." );
653 0 : return -1;
654 : }
655 :
656 5 : if( m_rows == m_psdSqrt->rows() && m_cols == m_psdSqrt->cols() )
657 : {
658 0 : return 0;
659 : }
660 :
661 5 : m_rows = m_psdSqrt->rows();
662 5 : m_cols = m_psdSqrt->cols();
663 5 : m_planes = 1;
664 :
665 5 : m_ftWork.resize( m_rows, m_cols );
666 :
667 5 : m_fft_fwd.plan( m_rows, m_cols, math::ft::dir::forward, true );
668 :
669 5 : m_fft_bwd.plan( m_rows, m_cols, math::ft::dir::backward, true );
670 :
671 5 : return 0;
672 : }
673 :
674 : template <typename realT, size_t rank>
675 : template <size_t crank>
676 5 : int psdFilter<realT, rank>::setSize( typename std::enable_if<crank == 3>::type * )
677 : {
678 5 : if( m_psdSqrt == 0 )
679 : {
680 0 : internal::mxlib_error_report(error_t::paramnotset, "m_psdSqrt has not been set yet, is still NULL." );
681 0 : return -1;
682 : }
683 :
684 5 : if( m_rows == m_psdSqrt->rows() && m_cols == m_psdSqrt->cols() && m_planes == m_psdSqrt->planes() )
685 : {
686 0 : return 0;
687 : }
688 :
689 5 : m_rows = m_psdSqrt->rows();
690 5 : m_cols = m_psdSqrt->cols();
691 5 : m_planes = m_psdSqrt->planes();
692 :
693 5 : m_ftWork.resize( m_rows, m_cols, m_planes );
694 :
695 5 : m_fft_fwd.plan( m_planes, m_rows, m_cols, math::ft::dir::forward, true );
696 :
697 5 : m_fft_bwd.plan( m_planes, m_rows, m_cols, math::ft::dir::backward, true );
698 :
699 5 : return 0;
700 : }
701 :
702 : template <typename realT, size_t rank>
703 30 : int psdFilter<realT, rank>::rows()
704 : {
705 30 : return m_rows;
706 : }
707 :
708 : template <typename realT, size_t rank>
709 28 : int psdFilter<realT, rank>::cols()
710 : {
711 28 : return m_cols;
712 : }
713 :
714 : template <typename realT, size_t rank>
715 26 : int psdFilter<realT, rank>::planes()
716 : {
717 26 : return m_planes;
718 : }
719 :
720 : template <typename realT, size_t rank>
721 : template <size_t crank>
722 1 : int psdFilter<realT, rank>::psdSqrt( realArrayT *npsdSqrt, realT df, typename std::enable_if<crank == 1>::type * )
723 : {
724 1 : m_dFreq1 = df;
725 1 : return psdSqrt( npsdSqrt );
726 : }
727 :
728 : template <typename realT, size_t rank>
729 : template <size_t crank>
730 1 : int psdFilter<realT, rank>::psdSqrt( realArrayT *npsdSqrt,
731 : realT dk1,
732 : realT dk2,
733 : typename std::enable_if<crank == 2>::type * )
734 : {
735 1 : m_dFreq1 = dk1;
736 1 : m_dFreq2 = dk2;
737 1 : return psdSqrt( npsdSqrt );
738 : }
739 :
740 : template <typename realT, size_t rank>
741 : template <size_t crank>
742 1 : int psdFilter<realT, rank>::psdSqrt(
743 : realArrayT *npsdSqrt, realT dk1, realT dk2, realT df, typename std::enable_if<crank == 3>::type * )
744 : {
745 1 : m_dFreq1 = dk1;
746 1 : m_dFreq2 = dk2;
747 1 : m_dFreq3 = df;
748 1 : return psdSqrt( npsdSqrt );
749 : }
750 :
751 : template <typename realT, size_t rank>
752 : template <size_t crank>
753 1 : int psdFilter<realT, rank>::psdSqrt( const realArrayT &npsdSqrt, realT df, typename std::enable_if<crank == 1>::type * )
754 : {
755 1 : m_dFreq1 = df;
756 1 : return psdSqrt( npsdSqrt );
757 : }
758 :
759 : template <typename realT, size_t rank>
760 : template <size_t crank>
761 1 : int psdFilter<realT, rank>::psdSqrt( const realArrayT &npsdSqrt,
762 : realT dk1,
763 : realT dk2,
764 : typename std::enable_if<crank == 2>::type * )
765 : {
766 1 : m_dFreq1 = dk1;
767 1 : m_dFreq2 = dk2;
768 1 : return psdSqrt( npsdSqrt );
769 : }
770 :
771 : template <typename realT, size_t rank>
772 : template <size_t crank>
773 1 : int psdFilter<realT, rank>::psdSqrt(
774 : const realArrayT &npsdSqrt, realT dk1, realT dk2, realT df, typename std::enable_if<crank == 3>::type * )
775 : {
776 1 : m_dFreq1 = dk1;
777 1 : m_dFreq2 = dk2;
778 1 : m_dFreq3 = df;
779 1 : return psdSqrt( npsdSqrt );
780 : }
781 :
782 : template <typename realT, size_t rank>
783 : template <size_t crank>
784 3 : int psdFilter<realT, rank>::psd( const realArrayT &npsd, const realT df1, typename std::enable_if<crank == 1>::type * )
785 : {
786 3 : if( m_psdSqrt && m_owner )
787 : {
788 0 : delete m_psdSqrt;
789 : }
790 :
791 3 : m_psdSqrt = new realArrayT;
792 :
793 : // Vector
794 3 : m_psdSqrt->resize( npsd.size() );
795 7171 : for( size_t n = 0; n < npsd.size(); ++n )
796 7168 : ( *m_psdSqrt )[n] = sqrt( npsd[n] );
797 :
798 3 : m_owner = true;
799 :
800 3 : m_dFreq1 = df1;
801 :
802 3 : setSize();
803 :
804 3 : return 0;
805 : }
806 :
807 : template <typename realT, size_t rank>
808 : template <size_t crank>
809 3 : int psdFilter<realT, rank>::psd( const realArrayT &npsd,
810 : const realT dk1,
811 : const realT dk2,
812 : typename std::enable_if<crank == 2>::type * )
813 : {
814 3 : if( m_psdSqrt && m_owner )
815 : {
816 0 : delete m_psdSqrt;
817 : }
818 :
819 3 : m_psdSqrt = new realArrayT;
820 :
821 3 : ( *m_psdSqrt ) = npsd.sqrt();
822 3 : m_owner = true;
823 :
824 3 : m_dFreq1 = dk1;
825 3 : m_dFreq2 = dk2;
826 :
827 3 : setSize();
828 :
829 3 : return 0;
830 : }
831 :
832 : template <typename realT, size_t rank>
833 : template <size_t crank>
834 3 : int psdFilter<realT, rank>::psd( const realArrayT &npsd,
835 : const realT dk1,
836 : const realT dk2,
837 : const realT df,
838 : typename std::enable_if<crank == 3>::type * )
839 : {
840 3 : if( m_psdSqrt && m_owner )
841 : {
842 0 : delete m_psdSqrt;
843 : }
844 :
845 3 : m_psdSqrt = new realArrayT;
846 :
847 : // Cube
848 3 : m_psdSqrt->resize( npsd.rows(), npsd.cols(), npsd.planes() );
849 387 : for( int pp = 0; pp < npsd.planes(); ++pp )
850 : {
851 37248 : for( int cc = 0; cc < npsd.cols(); ++cc )
852 : {
853 4362240 : for( int rr = 0; rr < npsd.rows(); ++rr )
854 : {
855 4325376 : m_psdSqrt->image( pp )( rr, cc ) = sqrt( npsd.image( pp )( rr, cc ) );
856 : }
857 : }
858 : }
859 :
860 3 : m_owner = true;
861 :
862 3 : m_dFreq1 = dk1;
863 3 : m_dFreq2 = dk2;
864 3 : m_dFreq3 = df;
865 :
866 3 : setSize();
867 :
868 3 : return 0;
869 : }
870 :
871 : template <typename realT, size_t rank>
872 9 : void psdFilter<realT, rank>::clear()
873 : {
874 : // m_ftWork.resize(0,0);
875 9 : psdFilterTypes::arrayT<realT, rank>::clear( m_ftWork );
876 :
877 9 : m_rows = 0;
878 9 : m_cols = 0;
879 9 : m_planes = 0;
880 :
881 9 : if( m_psdSqrt && m_owner )
882 : {
883 6 : delete m_psdSqrt;
884 6 : m_psdSqrt = 0;
885 : }
886 9 : }
887 :
888 : template <typename realT, size_t rank>
889 : template <size_t crank>
890 200 : int psdFilter<realT, rank>::filter( realArrayT &noise,
891 : realArrayT *noiseIm,
892 : typename std::enable_if<crank == 1>::type * ) const
893 : {
894 614600 : for( int nn = 0; nn < noise.size(); ++nn )
895 614400 : m_ftWork[nn] = complexT( noise[nn], 0 );
896 :
897 : // Transform complex noise to Fourier domain.
898 200 : m_fft_fwd( m_ftWork.data(), m_ftWork.data() );
899 :
900 : // Apply the filter.
901 614600 : for( int nn = 0; nn < m_ftWork.size(); ++nn )
902 614400 : m_ftWork[nn] *= ( *m_psdSqrt )[nn];
903 :
904 200 : m_fft_bwd( m_ftWork.data(), m_ftWork.data() );
905 :
906 : // Now take the real part, and normalize.
907 200 : realT norm = sqrt( noise.size() / m_dFreq1 );
908 614600 : for( int nn = 0; nn < m_ftWork.size(); ++nn )
909 614400 : noise[nn] = m_ftWork[nn].real() / norm;
910 :
911 200 : if( noiseIm != nullptr )
912 : {
913 0 : for( int nn = 0; nn < m_ftWork.size(); ++nn )
914 0 : ( *noiseIm )[nn] = m_ftWork[nn].imag() / norm;
915 : }
916 :
917 200 : return 0;
918 : }
919 :
920 : template <typename realT, size_t rank>
921 : template <size_t crank>
922 200 : int psdFilter<realT, rank>::filter( realArrayT &noise,
923 : realArrayT *noiseIm,
924 : typename std::enable_if<crank == 2>::type * ) const
925 : {
926 : // Make noise a complex number
927 13000 : for( int ii = 0; ii < noise.rows(); ++ii )
928 : {
929 832000 : for( int jj = 0; jj < noise.cols(); ++jj )
930 : {
931 819200 : m_ftWork( ii, jj ) = complexT( noise( ii, jj ), 0 );
932 : }
933 : }
934 :
935 : // Transform complex noise to Fourier domain.
936 200 : m_fft_fwd( m_ftWork.data(), m_ftWork.data() );
937 :
938 : // Apply the filter.
939 200 : m_ftWork *= *m_psdSqrt;
940 :
941 200 : m_fft_bwd( m_ftWork.data(), m_ftWork.data() );
942 :
943 200 : realT norm = sqrt( noise.rows() * noise.cols() / ( m_dFreq1 * m_dFreq2 ) );
944 :
945 : // Now take the real part, and normalize.
946 200 : noise = m_ftWork.real() / norm;
947 :
948 200 : if( noiseIm != nullptr )
949 : {
950 0 : *noiseIm = m_ftWork.imag() / norm;
951 : }
952 :
953 200 : return 0;
954 : }
955 :
956 : template <typename realT, size_t rank>
957 : template <size_t crank>
958 : int psdFilter<realT, rank>::filter( realArrayMapT noise,
959 : realArrayT *noiseIm,
960 : typename std::enable_if<crank == 2>::type * ) const
961 : {
962 : // Make noise a complex number
963 : for( int ii = 0; ii < noise.rows(); ++ii )
964 : {
965 : for( int jj = 0; jj < noise.cols(); ++jj )
966 : {
967 : m_ftWork( ii, jj ) = complexT( noise( ii, jj ), 0 );
968 : }
969 : }
970 :
971 : // Transform complex noise to Fourier domain.
972 : m_fft_fwd( m_ftWork.data(), m_ftWork.data() );
973 :
974 : // Apply the filter.
975 : m_ftWork *= *m_psdSqrt;
976 :
977 : m_fft_bwd( m_ftWork.data(), m_ftWork.data() );
978 :
979 : realT norm = sqrt( noise.rows() * noise.cols() / ( m_dFreq1 * m_dFreq2 ) );
980 :
981 : // Now take the real part, and normalize.
982 : noise = m_ftWork.real() / norm;
983 :
984 : if( noiseIm != nullptr )
985 : {
986 : *noiseIm = m_ftWork.imag() / norm;
987 : }
988 :
989 : return 0;
990 : }
991 :
992 : template <typename realT, size_t rank>
993 : template <size_t crank>
994 200 : int psdFilter<realT, rank>::filter( realArrayT &noise,
995 : realArrayT *noiseIm,
996 : typename std::enable_if<crank == 3>::type * ) const
997 : {
998 : // Make noise a complex number
999 13000 : for( int pp = 0; pp < noise.planes(); ++pp )
1000 : {
1001 422400 : for( int ii = 0; ii < noise.rows(); ++ii )
1002 : {
1003 13516800 : for( int jj = 0; jj < noise.cols(); ++jj )
1004 : {
1005 13107200 : m_ftWork.image( pp )( ii, jj ) = complexT( noise.image( pp )( ii, jj ), 0 );
1006 : }
1007 : }
1008 : }
1009 :
1010 : // Transform complex noise to Fourier domain.
1011 200 : m_fft_fwd( m_ftWork.data(), m_ftWork.data() );
1012 :
1013 : // Apply the filter.
1014 13000 : for( int pp = 0; pp < noise.planes(); ++pp )
1015 12800 : m_ftWork.image( pp ) *= m_psdSqrt->image( pp );
1016 :
1017 200 : m_fft_bwd( m_ftWork.data(), m_ftWork.data() );
1018 :
1019 : // Now take the real part, and normalize.
1020 :
1021 200 : realT norm = sqrt( m_rows * m_cols * m_planes / ( m_dFreq1 * m_dFreq2 * m_dFreq3 ) );
1022 13000 : for( int pp = 0; pp < noise.planes(); ++pp )
1023 12800 : noise.image( pp ) = m_ftWork.image( pp ).real() / norm;
1024 :
1025 200 : if( noiseIm != nullptr )
1026 : {
1027 0 : for( int pp = 0; pp < noise.planes(); ++pp )
1028 0 : noiseIm->image( pp ) = m_ftWork.image( pp ).imag() / norm;
1029 : }
1030 :
1031 200 : return 0;
1032 : }
1033 :
1034 : template <typename realT, size_t rank>
1035 600 : int psdFilter<realT, rank>::operator()( realArrayT &noise ) const
1036 : {
1037 600 : return filter( noise );
1038 : }
1039 :
1040 : template <typename realT, size_t rank>
1041 : int psdFilter<realT, rank>::operator()( realArrayMapT noise ) const
1042 : {
1043 : return filter( noise );
1044 : }
1045 :
1046 : template <typename realT, size_t rank>
1047 : int psdFilter<realT, rank>::operator()( realArrayT &noise, realArrayT &noiseIm ) const
1048 : {
1049 : return filter( noise, &noiseIm );
1050 : }
1051 :
1052 : } // namespace sigproc
1053 : } // namespace mx
1054 :
1055 : #endif // psdFilter_hpp
|
