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