mxlib
c++ tools for analyzing astronomical data and other tasks by Jared R. Males. [git repo]
Loading...
Searching...
No Matches
psdFilter_test.cpp
Go to the documentation of this file.
1/** \file psdFilter_test.cpp
2 */
3#include "../../catch2/catch.hpp"
4
5#include <vector>
6#include <Eigen/Dense>
7
8#define MX_NO_ERROR_REPORTS
9
10#include "../../../include/sigproc/psdFilter.hpp"
11#include "../../../include/sigproc/psdUtils.hpp"
12#include "../../../include/improc/eigenCube.hpp"
13#include "../../../include/math/randomT.hpp"
14#include "../../../include/math/vectorUtils.hpp"
15
16#ifdef MXLIB_BUILD_COVERAGE
17constexpr int psdFilterTrials = 100;
18constexpr double psdFilterTol = 0.09;
19#else
20constexpr int psdFilterTrials = 10000;
21constexpr double psdFilterTol = 0.02;
22#endif
23
24/** Scenario: compiling psdFilter
25 *
26 * Verify compilation and initilization of the 3 ranks for psdFilter.
27 *
28 * \anchor tests_sigproc_psdFilter_compile
29 */
30SCENARIO( "compiling psdFilter", "[sigproc::psdFilter]" )
31{
32 GIVEN( "a psdFilter, sqrt pointer" )
33 {
34 WHEN( "rank==1" )
35 {
36 mx::sigproc::psdFilter<double, 1> psdF;
37
38 std::vector<double> psdSqrt( 1024, 1 );
39
40 int rv = psdF.psdSqrt( &psdSqrt, 1 );
41
42 REQUIRE( rv == 0 );
43 REQUIRE( psdF.rows() == 1024 );
44 REQUIRE( psdF.cols() == 1 );
45 REQUIRE( psdF.planes() == 1 );
46
47 psdF.clear();
48 REQUIRE( psdF.rows() == 0 );
49 REQUIRE( psdF.cols() == 0 );
50 REQUIRE( psdF.planes() == 0 );
51 }
52 WHEN( "rank==2" )
53 {
54 mx::sigproc::psdFilter<double, 2> psdF;
55
56 Eigen::Array<double, -1, -1> psdSqrt;
57 psdSqrt.resize( 256, 256 );
58 psdSqrt.setConstant( 1 );
59
60 int rv = psdF.psdSqrt( &psdSqrt, 1, 1 );
61
62 REQUIRE( rv == 0 );
63 REQUIRE( psdF.rows() == 256 );
64 REQUIRE( psdF.cols() == 256 );
65 REQUIRE( psdF.planes() == 1 );
66
67 psdF.clear();
68 REQUIRE( psdF.rows() == 0 );
69 REQUIRE( psdF.cols() == 0 );
70 REQUIRE( psdF.planes() == 0 );
71 }
72 WHEN( "rank==3" )
73 {
74 mx::sigproc::psdFilter<double, 3> psdF;
75
77 psdSqrt.resize( 128, 128, 256 );
78 // psdSqrt.setConstant(1);
79
80 int rv = psdF.psdSqrt( &psdSqrt, 1, 1, 1 );
81
82 REQUIRE( rv == 0 );
83 REQUIRE( psdF.rows() == 128 );
84 REQUIRE( psdF.cols() == 128 );
85 REQUIRE( psdF.planes() == 256 );
86
87 psdF.clear();
88 REQUIRE( psdF.rows() == 0 );
89 REQUIRE( psdF.cols() == 0 );
90 REQUIRE( psdF.planes() == 0 );
91 }
92 }
93
94 GIVEN( "a psdFilter, sqrt reference" )
95 {
96 WHEN( "rank==1" )
97 {
98 mx::sigproc::psdFilter<double, 1> psdF;
99
100 std::vector<double> psdSqrt( 1024, 1 );
101
102 int rv = psdF.psdSqrt( psdSqrt, 1 );
103
104 REQUIRE( rv == 0 );
105 REQUIRE( psdF.rows() == 1024 );
106 REQUIRE( psdF.cols() == 1 );
107 REQUIRE( psdF.planes() == 1 );
108
109 psdF.clear();
110 REQUIRE( psdF.rows() == 0 );
111 REQUIRE( psdF.cols() == 0 );
112 REQUIRE( psdF.planes() == 0 );
113 }
114 WHEN( "rank==2" )
115 {
116 mx::sigproc::psdFilter<double, 2> psdF;
117
118 Eigen::Array<double, -1, -1> psdSqrt;
119 psdSqrt.resize( 256, 256 );
120 psdSqrt.setConstant( 1 );
121
122 int rv = psdF.psdSqrt( psdSqrt, 1, 1 );
123
124 REQUIRE( rv == 0 );
125 REQUIRE( psdF.rows() == 256 );
126 REQUIRE( psdF.cols() == 256 );
127 REQUIRE( psdF.planes() == 1 );
128
129 psdF.clear();
130 REQUIRE( psdF.rows() == 0 );
131 REQUIRE( psdF.cols() == 0 );
132 REQUIRE( psdF.planes() == 0 );
133 }
134 WHEN( "rank==3" )
135 {
136 mx::sigproc::psdFilter<double, 3> psdF;
137
139 psdSqrt.resize( 128, 128, 256 );
140 // psdSqrt.setConstant(1);
141
142 int rv = psdF.psdSqrt( psdSqrt, 1, 1, 1 );
143
144 REQUIRE( rv == 0 );
145 REQUIRE( psdF.rows() == 128 );
146 REQUIRE( psdF.cols() == 128 );
147 REQUIRE( psdF.planes() == 256 );
148
149 psdF.clear();
150 REQUIRE( psdF.rows() == 0 );
151 REQUIRE( psdF.cols() == 0 );
152 REQUIRE( psdF.planes() == 0 );
153 }
154 }
155
156 GIVEN( "a psdFilter, psd reference" )
157 {
158 WHEN( "rank==1" )
159 {
160 mx::sigproc::psdFilter<double, 1> psdF;
161
162 std::vector<double> psd( 1024, 1 );
163
164 int rv = psdF.psd( psd, 1.0 );
165
166 REQUIRE( rv == 0 );
167 REQUIRE( psdF.rows() == 1024 );
168 REQUIRE( psdF.cols() == 1 );
169 REQUIRE( psdF.planes() == 1 );
170
171 psdF.clear();
172 REQUIRE( psdF.rows() == 0 );
173 REQUIRE( psdF.cols() == 0 );
174 REQUIRE( psdF.planes() == 0 );
175 }
176 WHEN( "rank==2" )
177 {
178 mx::sigproc::psdFilter<double, 2> psdF;
179
180 Eigen::Array<double, -1, -1> psd;
181 psd.resize( 256, 256 );
182 psd.setConstant( 1 );
183
184 int rv = psdF.psd( psd, 1.0, 1.0 );
185
186 REQUIRE( rv == 0 );
187 REQUIRE( psdF.rows() == 256 );
188 REQUIRE( psdF.cols() == 256 );
189 REQUIRE( psdF.planes() == 1 );
190
191 psdF.clear();
192 REQUIRE( psdF.rows() == 0 );
193 REQUIRE( psdF.cols() == 0 );
194 REQUIRE( psdF.planes() == 0 );
195 }
196 WHEN( "rank==3" )
197 {
198 mx::sigproc::psdFilter<double, 3> psdF;
199
201 psd.resize( 128, 128, 256 );
202 // psdSqrt.setConstant(1);
203
204 int rv = psdF.psd( psd, 1, 1, 1 );
205
206 REQUIRE( rv == 0 );
207 REQUIRE( psdF.rows() == 128 );
208 REQUIRE( psdF.cols() == 128 );
209 REQUIRE( psdF.planes() == 256 );
210
211 psdF.clear();
212 REQUIRE( psdF.rows() == 0 );
213 REQUIRE( psdF.cols() == 0 );
214 REQUIRE( psdF.planes() == 0 );
215 }
216 }
217}
218
219/** Verify filtering and noise normalization
220 * Conducts random noise tests, verifying that the resultant rms is within 2% of expected value on average over many
221 * trials. Results are usually better than 1%, but 2% makes sure we don't get false failures.
222 *
223 * \anchor tests_sigproc_psdFilter_filter
224 */
225SCENARIO( "filtering with psdFilter", "[sigproc::psdFilter]" )
226{
227 GIVEN( "a rank 1 psd" )
228 {
229 WHEN( "alpha=-2.5, df Nyquist matched to array size, var=1" )
230 {
231 mx::sigproc::psdFilter<double, 1> psdF;
232
233 std::vector<double> f( 2049 ), psd( 2049 );
234
235 for( size_t n = 0; n < psd.size(); ++n )
236 f[n] = n * 1.0 / 4096.;
237
238 for( size_t n = 1; n < psd.size(); ++n )
239 psd[n] = pow( f[n], -2.5 );
240 psd[0] = psd[1];
241
242 mx::sigproc::normPSD( psd, f, 1.0, -1e5, 1e5 );
243
244 std::vector<double> f2s, psd2s;
246 mx::sigproc::augment1SidedPSD( psd2s, psd );
247
248 double df = f2s[1] - f2s[0];
249 // mx::sigproc::normPSD(psd2s, f2s, 1.0, -1e5, 1e5);
250 int rv = psdF.psd( psd2s, df );
251
252 REQUIRE( rv == 0 );
253 REQUIRE( psdF.rows() == 2. * psd.size() - 2 );
254 REQUIRE( psdF.cols() == 1 );
255 REQUIRE( psdF.planes() == 1 );
256
257 std::vector<double> noise( psdF.rows() );
258
260
261 double avgRms = 0;
262
263 for( int k = 0; k < psdFilterTrials; ++k )
264 {
265 for( size_t n = 0; n < noise.size(); ++n )
266 noise[n] = normVar;
267 psdF( noise );
268 avgRms += ( mx::math::vectorVariance( noise, 0.0 ) );
269 }
270
271 avgRms = sqrt( avgRms / psdFilterTrials );
272
273 REQUIRE_THAT( avgRms, Catch::Matchers::WithinAbs( 1.0, psdFilterTol ) );
274 }
275 WHEN( "alpha=-1.5, df arbitrary, var = 2.2" )
276 {
277 mx::sigproc::psdFilter<double, 1> psdF;
278
279 std::vector<double> f( 1025 ), psd( 1025 );
280
281 for( size_t n = 0; n < psd.size(); ++n )
282 f[n] = n * 1.0 / 7000.;
283
284 for( size_t n = 1; n < psd.size(); ++n )
285 psd[n] = pow( f[n], -1.5 );
286 psd[0] = psd[1];
287
288 std::vector<double> f2s, psd2s;
290 mx::sigproc::augment1SidedPSD( psd2s, psd );
291
292 double df = f2s[1] - f2s[0];
293 mx::sigproc::normPSD( psd2s, f2s, 2.2, -1e5, 1e5 );
294
295 int rv = psdF.psd( psd2s, df );
296
297 REQUIRE( rv == 0 );
298 REQUIRE( psdF.rows() == 2. * psd.size() - 2 );
299 REQUIRE( psdF.cols() == 1 );
300 REQUIRE( psdF.planes() == 1 );
301
302 std::vector<double> noise( psdF.rows() );
303
305
306 double avgRms = 0;
307
308 for( int k = 0; k < psdFilterTrials; ++k )
309 {
310 for( size_t n = 0; n < noise.size(); ++n )
311 noise[n] = normVar;
312 psdF( noise );
313 avgRms += ( mx::math::vectorVariance( noise, 0.0 ) );
314 }
315
316 avgRms = sqrt( avgRms / psdFilterTrials );
317
318 REQUIRE_THAT( avgRms, Catch::Matchers::WithinAbs( sqrt( 2.2 ), psdFilterTol * sqrt( 2.2 ) ) );
319 }
320 }
321 GIVEN( "a rank 2 psd" )
322 {
323 WHEN( "alpha=-2.5, dk Nyquist matched to array size, var=1" )
324 {
325 mx::sigproc::psdFilter<double, 2> psdF;
326
327 Eigen::Array<double, -1, -1> k, psd;
328
329 k.resize( 64, 64 );
330 psd.resize( 64, 64 );
331
332 mx::sigproc::frequencyGrid( k, 1. / 128. );
333 for( int cc = 0; cc < psd.cols(); ++cc )
334 {
335 for( int rr = 0; rr < psd.rows(); ++rr )
336 {
337 if( k( rr, cc ) == 0 )
338 psd( rr, cc ) = 0;
339 else
340 psd( rr, cc ) = pow( k( rr, cc ), -2.5 );
341 }
342 }
343
344 double dk = k( 0, 1 ) - k( 0, 0 );
345
346 mx::sigproc::normPSD( psd, k, 1.0 );
347
348 int rv = psdF.psd( psd, dk, dk );
349
350 REQUIRE( rv == 0 );
351 REQUIRE( psdF.rows() == psd.rows() );
352 REQUIRE( psdF.cols() == psd.cols() );
353 REQUIRE( psdF.planes() == 1 );
354
355 Eigen::Array<double, -1, -1> noise( psdF.rows(), psdF.cols() );
356
358
359 double avgRms = 0;
360
361 for( int k = 0; k < psdFilterTrials; ++k )
362 {
363 for( int cc = 0; cc < psd.cols(); ++cc )
364 {
365 for( int rr = 0; rr < psd.rows(); ++rr )
366 {
367 noise( rr, cc ) = normVar;
368 }
369 }
370
371 psdF( noise );
372 avgRms += noise.square().sum(); //(mx::math::vectorVariance(noise,0.0));
373 }
374
375 avgRms = sqrt( avgRms / ( psd.rows() * psd.cols() ) / psdFilterTrials );
376
377 REQUIRE_THAT( avgRms, Catch::Matchers::WithinAbs( 1.0, psdFilterTol ) );
378 }
379 WHEN( "alpha=-1.5, dk arb, var=2.2" )
380 {
381 mx::sigproc::psdFilter<double, 2> psdF;
382
383 Eigen::Array<double, -1, -1> k, psd;
384
385 k.resize( 64, 64 );
386 psd.resize( 64, 64 );
387
388 mx::sigproc::frequencyGrid( k, 1. / 302. );
389 for( int cc = 0; cc < psd.cols(); ++cc )
390 {
391 for( int rr = 0; rr < psd.rows(); ++rr )
392 {
393 if( k( rr, cc ) == 0 )
394 psd( rr, cc ) = 0;
395 else
396 psd( rr, cc ) = pow( k( rr, cc ), -1.5 );
397 }
398 }
399
400 double dk = k( 0, 1 ) - k( 0, 0 );
401
402 mx::sigproc::normPSD( psd, k, 2.2 );
403
404 int rv = psdF.psd( psd, dk, dk );
405
406 REQUIRE( rv == 0 );
407 REQUIRE( psdF.rows() == psd.rows() );
408 REQUIRE( psdF.cols() == psd.cols() );
409 REQUIRE( psdF.planes() == 1 );
410
411 Eigen::Array<double, -1, -1> noise( psdF.rows(), psdF.cols() );
412
414
415 double avgRms = 0;
416
417 for( int k = 0; k < psdFilterTrials; ++k )
418 {
419 for( int cc = 0; cc < psd.cols(); ++cc )
420 {
421 for( int rr = 0; rr < psd.rows(); ++rr )
422 {
423 noise( rr, cc ) = normVar;
424 }
425 }
426
427 psdF( noise );
428 avgRms += noise.square().sum(); //(mx::math::vectorVariance(noise,0.0));
429 }
430
431 avgRms = sqrt( avgRms / ( psd.rows() * psd.cols() ) / psdFilterTrials );
432
433 REQUIRE_THAT( avgRms, Catch::Matchers::WithinAbs( sqrt( 2.2 ), psdFilterTol * sqrt( 2.2 ) ) );
434 }
435 }
436 GIVEN( "a rank 3 psd" )
437 {
438 WHEN( "k-alpha=-2.5, f-alph=-2.5, dk Nyquist matched to array size, df Nyquist matched to array size, var=1" )
439 {
440 mx::sigproc::psdFilter<double, 3> psdF;
441
442 Eigen::Array<double, -1, -1> k, psdk;
443 std::vector<double> f, f2s, psd2s;
444
446
447 k.resize( 32, 32 );
448 f.resize( 33 );
449
450 mx::sigproc::frequencyGrid( k, 1. / 64. );
451 psdk.resize( k.rows(), k.cols() );
452 for( int cc = 0; cc < psdk.cols(); ++cc )
453 {
454 for( int rr = 0; rr < psdk.rows(); ++rr )
455 {
456 if( k( rr, cc ) == 0 )
457 psdk( rr, cc ) = 0;
458 else
459 psdk( rr, cc ) = pow( k( rr, cc ), -2.5 );
460 }
461 }
462 mx::sigproc::normPSD( psdk, k, 1.0 );
463
464 for( size_t n = 0; n < f.size(); ++n )
465 f[n] = n * 1.0 / 64.;
467 psd2s.resize( f2s.size() );
468 for( size_t n = 0; n < psd2s.size(); ++n )
469 psd2s[n] = pow( fabs( f2s[n] ), -2.5 );
470 psd2s[0] = psd2s[1];
471
472 psd.resize( k.rows(), k.cols(), f2s.size() );
473
474 for( int cc = 0; cc < psd.cols(); ++cc )
475 {
476 for( int rr = 0; rr < psd.rows(); ++rr )
477 {
478 if( k( rr, cc ) == 0 )
479 psd.pixel( rr, cc ).setZero();
480 else
481 {
482 double p = psdk( rr, cc );
483 mx::sigproc::normPSD( psd2s, f2s, p, -1e5, 1e5 );
484
485 for( int pp = 0; pp < psd.planes(); ++pp )
486 psd.image( pp )( rr, cc ) = psd2s[pp];
487 }
488 }
489 }
490
491 double dk = k( 0, 1 ) - k( 0, 0 );
492 double df = f[1] - f[0];
493
494 int rv = psdF.psd( psd, dk, dk, df );
495
496 REQUIRE( rv == 0 );
497 REQUIRE( psdF.rows() == psd.rows() );
498 REQUIRE( psdF.cols() == psd.cols() );
499 REQUIRE( psdF.planes() == psd.planes() );
500
501 mx::improc::eigenCube<double> noise( psdF.rows(), psdF.cols(), psdF.planes() );
502
504
505 double avgRms = 0;
506
507 for( int k = 0; k < psdFilterTrials; ++k )
508 {
509 for( int pp = 0; pp < psd.planes(); ++pp )
510 {
511 for( int cc = 0; cc < psd.cols(); ++cc )
512 {
513 for( int rr = 0; rr < psd.rows(); ++rr )
514 {
515 noise.image( pp )( rr, cc ) = normVar;
516 }
517 }
518 }
519
520 psdF( noise );
521 for( int pp = 0; pp < noise.planes(); ++pp )
522 avgRms += noise.image( pp ).square().sum();
523 }
524
525 avgRms = sqrt( avgRms / ( psd.rows() * psd.cols() * psd.planes() ) / psdFilterTrials );
526
527 REQUIRE_THAT( avgRms, Catch::Matchers::WithinAbs( 1.0, psdFilterTol ) );
528 }
529 WHEN( "k-alpha=-3.5, f-alph=-1.5, dk arb, df arb, var=2" )
530 {
531 mx::sigproc::psdFilter<double, 3> psdF;
532
533 Eigen::Array<double, -1, -1> k, psdk;
534 std::vector<double> f, f2s, psd2s;
535
537
538 k.resize( 32, 32 );
539 f.resize( 33 );
540
541 mx::sigproc::frequencyGrid( k, 1. / 640. );
542 psdk.resize( k.rows(), k.cols() );
543 for( int cc = 0; cc < psdk.cols(); ++cc )
544 {
545 for( int rr = 0; rr < psdk.rows(); ++rr )
546 {
547 if( k( rr, cc ) == 0 )
548 psdk( rr, cc ) = 0;
549 else
550 psdk( rr, cc ) = pow( k( rr, cc ), -3.5 );
551 }
552 }
553 mx::sigproc::normPSD( psdk, k, 2.0 );
554
555 for( size_t n = 0; n < f.size(); ++n )
556 f[n] = n * 1.0 / 78.;
558 psd2s.resize( f2s.size() );
559 for( size_t n = 0; n < psd2s.size(); ++n )
560 psd2s[n] = pow( fabs( f2s[n] ), -1.5 );
561 psd2s[0] = psd2s[1];
562
563 psd.resize( k.rows(), k.cols(), f2s.size() );
564
565 for( int cc = 0; cc < psd.cols(); ++cc )
566 {
567 for( int rr = 0; rr < psd.rows(); ++rr )
568 {
569 if( k( rr, cc ) == 0 )
570 psd.pixel( rr, cc ).setZero();
571 else
572 {
573 double p = psdk( rr, cc );
574 mx::sigproc::normPSD( psd2s, f2s, p, -1e5, 1e5 );
575
576 for( int pp = 0; pp < psd.planes(); ++pp )
577 psd.image( pp )( rr, cc ) = psd2s[pp];
578 }
579 }
580 }
581
582 double dk = k( 0, 1 ) - k( 0, 0 );
583 double df = f[1] - f[0];
584
585 int rv = psdF.psd( psd, dk, dk, df );
586
587 REQUIRE( rv == 0 );
588 REQUIRE( psdF.rows() == psd.rows() );
589 REQUIRE( psdF.cols() == psd.cols() );
590 REQUIRE( psdF.planes() == psd.planes() );
591
592 mx::improc::eigenCube<double> noise( psdF.rows(), psdF.cols(), psdF.planes() );
593
595
596 double avgRms = 0;
597
598 for( int k = 0; k < psdFilterTrials; ++k )
599 {
600 for( int pp = 0; pp < psd.planes(); ++pp )
601 {
602 for( int cc = 0; cc < psd.cols(); ++cc )
603 {
604 for( int rr = 0; rr < psd.rows(); ++rr )
605 {
606 noise.image( pp )( rr, cc ) = normVar;
607 }
608 }
609 }
610
611 psdF( noise );
612 for( int pp = 0; pp < noise.planes(); ++pp )
613 avgRms += noise.image( pp ).square().sum();
614 }
615
616 avgRms = sqrt( avgRms / ( psd.rows() * psd.cols() * psd.planes() ) / psdFilterTrials );
617
618 REQUIRE_THAT( avgRms, Catch::Matchers::WithinAbs( sqrt( 2.0 ), psdFilterTol * sqrt( 2 ) ) );
619 }
620 }
621}
An image cube with an Eigen-like API.
Definition eigenCube.hpp:30
Eigen::Map< Eigen::Array< dataT, Eigen::Dynamic, Eigen::Dynamic >, Eigen::Unaligned, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > pixel(Index i, Index j)
Returns an Eigen::Eigen::Map-ed vector of the pixels at the given coordinate.
Eigen::Map< Eigen::Array< dataT, Eigen::Dynamic, Eigen::Dynamic > > image(Index n)
Returns a 2D Eigen::Eigen::Map pointed at the specified image.
void augment1SidedPSD(vectorTout &psdTwoSided, vectorTin &psdOneSided, bool addZeroFreq=false, typename vectorTin::value_type scale=0.5)
Augment a 1-sided PSD to standard 2-sided FFT form.
Definition psdUtils.hpp:834
void augment1SidedPSDFreq(std::vector< T > &freqTwoSided, std::vector< T > &freqOneSided)
Augment a 1-sided frequency scale to standard FFT form.
Definition psdUtils.hpp:893
int normPSD(std::vector< floatT > &psd, std::vector< floatT > &f, floatParamT normT, floatT fmin=std::numeric_limits< floatT >::min(), floatT fmax=std::numeric_limits< floatT >::max())
Normalize a 1-D PSD to have a given variance.
Definition psdUtils.hpp:453
int frequencyGrid(std::vector< realT > &vec, realParamT dt, bool fftOrder=true)
Create a 1-D frequency grid.
Definition psdUtils.hpp:262
valueT vectorVariance(const valueT *vec, size_t sz, valueT mean)
Calculate the variance of a vector relative to a supplied mean value.
SCENARIO("compiling psdFilter", "[sigproc::psdFilter]")