30SCENARIO(
"compiling psdFilter",
"[sigproc::psdFilter]" )
32 GIVEN(
"a psdFilter, sqrt pointer" )
36 mx::sigproc::psdFilter<double, 1> psdF;
38 std::vector<double> psdSqrt( 1024, 1 );
40 int rv = psdF.psdSqrt( &psdSqrt, 1 );
43 REQUIRE( psdF.rows() == 1024 );
44 REQUIRE( psdF.cols() == 1 );
45 REQUIRE( psdF.planes() == 1 );
48 REQUIRE( psdF.rows() == 0 );
49 REQUIRE( psdF.cols() == 0 );
50 REQUIRE( psdF.planes() == 0 );
54 mx::sigproc::psdFilter<double, 2> psdF;
56 Eigen::Array<double, -1, -1> psdSqrt;
57 psdSqrt.resize( 256, 256 );
58 psdSqrt.setConstant( 1 );
60 int rv = psdF.psdSqrt( &psdSqrt, 1, 1 );
63 REQUIRE( psdF.rows() == 256 );
64 REQUIRE( psdF.cols() == 256 );
65 REQUIRE( psdF.planes() == 1 );
68 REQUIRE( psdF.rows() == 0 );
69 REQUIRE( psdF.cols() == 0 );
70 REQUIRE( psdF.planes() == 0 );
74 mx::sigproc::psdFilter<double, 3> psdF;
77 psdSqrt.resize( 128, 128, 256 );
80 int rv = psdF.psdSqrt( &psdSqrt, 1, 1, 1 );
83 REQUIRE( psdF.rows() == 128 );
84 REQUIRE( psdF.cols() == 128 );
85 REQUIRE( psdF.planes() == 256 );
88 REQUIRE( psdF.rows() == 0 );
89 REQUIRE( psdF.cols() == 0 );
90 REQUIRE( psdF.planes() == 0 );
94 GIVEN(
"a psdFilter, sqrt reference" )
98 mx::sigproc::psdFilter<double, 1> psdF;
100 std::vector<double> psdSqrt( 1024, 1 );
102 int rv = psdF.psdSqrt( psdSqrt, 1 );
105 REQUIRE( psdF.rows() == 1024 );
106 REQUIRE( psdF.cols() == 1 );
107 REQUIRE( psdF.planes() == 1 );
110 REQUIRE( psdF.rows() == 0 );
111 REQUIRE( psdF.cols() == 0 );
112 REQUIRE( psdF.planes() == 0 );
116 mx::sigproc::psdFilter<double, 2> psdF;
118 Eigen::Array<double, -1, -1> psdSqrt;
119 psdSqrt.resize( 256, 256 );
120 psdSqrt.setConstant( 1 );
122 int rv = psdF.psdSqrt( psdSqrt, 1, 1 );
125 REQUIRE( psdF.rows() == 256 );
126 REQUIRE( psdF.cols() == 256 );
127 REQUIRE( psdF.planes() == 1 );
130 REQUIRE( psdF.rows() == 0 );
131 REQUIRE( psdF.cols() == 0 );
132 REQUIRE( psdF.planes() == 0 );
136 mx::sigproc::psdFilter<double, 3> psdF;
139 psdSqrt.resize( 128, 128, 256 );
142 int rv = psdF.psdSqrt( psdSqrt, 1, 1, 1 );
145 REQUIRE( psdF.rows() == 128 );
146 REQUIRE( psdF.cols() == 128 );
147 REQUIRE( psdF.planes() == 256 );
150 REQUIRE( psdF.rows() == 0 );
151 REQUIRE( psdF.cols() == 0 );
152 REQUIRE( psdF.planes() == 0 );
156 GIVEN(
"a psdFilter, psd reference" )
160 mx::sigproc::psdFilter<double, 1> psdF;
162 std::vector<double> psd( 1024, 1 );
164 int rv = psdF.psd( psd, 1.0 );
167 REQUIRE( psdF.rows() == 1024 );
168 REQUIRE( psdF.cols() == 1 );
169 REQUIRE( psdF.planes() == 1 );
172 REQUIRE( psdF.rows() == 0 );
173 REQUIRE( psdF.cols() == 0 );
174 REQUIRE( psdF.planes() == 0 );
178 mx::sigproc::psdFilter<double, 2> psdF;
180 Eigen::Array<double, -1, -1> psd;
181 psd.resize( 256, 256 );
182 psd.setConstant( 1 );
184 int rv = psdF.psd( psd, 1.0, 1.0 );
187 REQUIRE( psdF.rows() == 256 );
188 REQUIRE( psdF.cols() == 256 );
189 REQUIRE( psdF.planes() == 1 );
192 REQUIRE( psdF.rows() == 0 );
193 REQUIRE( psdF.cols() == 0 );
194 REQUIRE( psdF.planes() == 0 );
198 mx::sigproc::psdFilter<double, 3> psdF;
201 psd.resize( 128, 128, 256 );
204 int rv = psdF.psd( psd, 1, 1, 1 );
207 REQUIRE( psdF.rows() == 128 );
208 REQUIRE( psdF.cols() == 128 );
209 REQUIRE( psdF.planes() == 256 );
212 REQUIRE( psdF.rows() == 0 );
213 REQUIRE( psdF.cols() == 0 );
214 REQUIRE( psdF.planes() == 0 );
225SCENARIO(
"filtering with psdFilter",
"[sigproc::psdFilter]" )
227 GIVEN(
"a rank 1 psd" )
229 WHEN(
"alpha=-2.5, df Nyquist matched to array size, var=1" )
231 mx::sigproc::psdFilter<double, 1> psdF;
233 std::vector<double> f( 2049 ), psd( 2049 );
235 for(
size_t n = 0; n < psd.size(); ++n )
236 f[n] = n * 1.0 / 4096.;
238 for(
size_t n = 1; n < psd.size(); ++n )
239 psd[n] = pow( f[n], -2.5 );
244 std::vector<double> f2s, psd2s;
248 double df = f2s[1] - f2s[0];
250 int rv = psdF.psd( psd2s, df );
253 REQUIRE( psdF.rows() == 2. * psd.size() - 2 );
254 REQUIRE( psdF.cols() == 1 );
255 REQUIRE( psdF.planes() == 1 );
257 std::vector<double> noise( psdF.rows() );
263 for(
int k = 0; k < psdFilterTrials; ++k )
265 for(
size_t n = 0; n < noise.size(); ++n )
271 avgRms = sqrt( avgRms / psdFilterTrials );
273 REQUIRE_THAT( avgRms, Catch::Matchers::WithinAbs( 1.0, psdFilterTol ) );
275 WHEN(
"alpha=-1.5, df arbitrary, var = 2.2" )
277 mx::sigproc::psdFilter<double, 1> psdF;
279 std::vector<double> f( 1025 ), psd( 1025 );
281 for(
size_t n = 0; n < psd.size(); ++n )
282 f[n] = n * 1.0 / 7000.;
284 for(
size_t n = 1; n < psd.size(); ++n )
285 psd[n] = pow( f[n], -1.5 );
288 std::vector<double> f2s, psd2s;
292 double df = f2s[1] - f2s[0];
295 int rv = psdF.psd( psd2s, df );
298 REQUIRE( psdF.rows() == 2. * psd.size() - 2 );
299 REQUIRE( psdF.cols() == 1 );
300 REQUIRE( psdF.planes() == 1 );
302 std::vector<double> noise( psdF.rows() );
308 for(
int k = 0; k < psdFilterTrials; ++k )
310 for(
size_t n = 0; n < noise.size(); ++n )
316 avgRms = sqrt( avgRms / psdFilterTrials );
318 REQUIRE_THAT( avgRms, Catch::Matchers::WithinAbs( sqrt( 2.2 ), psdFilterTol * sqrt( 2.2 ) ) );
321 GIVEN(
"a rank 2 psd" )
323 WHEN(
"alpha=-2.5, dk Nyquist matched to array size, var=1" )
325 mx::sigproc::psdFilter<double, 2> psdF;
327 Eigen::Array<double, -1, -1> k, psd;
330 psd.resize( 64, 64 );
333 for(
int cc = 0; cc < psd.cols(); ++cc )
335 for(
int rr = 0; rr < psd.rows(); ++rr )
337 if( k( rr, cc ) == 0 )
340 psd( rr, cc ) = pow( k( rr, cc ), -2.5 );
344 double dk = k( 0, 1 ) - k( 0, 0 );
348 int rv = psdF.psd( psd, dk, dk );
351 REQUIRE( psdF.rows() == psd.rows() );
352 REQUIRE( psdF.cols() == psd.cols() );
353 REQUIRE( psdF.planes() == 1 );
355 Eigen::Array<double, -1, -1> noise( psdF.rows(), psdF.cols() );
361 for(
int k = 0; k < psdFilterTrials; ++k )
363 for(
int cc = 0; cc < psd.cols(); ++cc )
365 for(
int rr = 0; rr < psd.rows(); ++rr )
367 noise( rr, cc ) = normVar;
372 avgRms += noise.square().sum();
375 avgRms = sqrt( avgRms / ( psd.rows() * psd.cols() ) / psdFilterTrials );
377 REQUIRE_THAT( avgRms, Catch::Matchers::WithinAbs( 1.0, psdFilterTol ) );
379 WHEN(
"alpha=-1.5, dk arb, var=2.2" )
381 mx::sigproc::psdFilter<double, 2> psdF;
383 Eigen::Array<double, -1, -1> k, psd;
386 psd.resize( 64, 64 );
389 for(
int cc = 0; cc < psd.cols(); ++cc )
391 for(
int rr = 0; rr < psd.rows(); ++rr )
393 if( k( rr, cc ) == 0 )
396 psd( rr, cc ) = pow( k( rr, cc ), -1.5 );
400 double dk = k( 0, 1 ) - k( 0, 0 );
404 int rv = psdF.psd( psd, dk, dk );
407 REQUIRE( psdF.rows() == psd.rows() );
408 REQUIRE( psdF.cols() == psd.cols() );
409 REQUIRE( psdF.planes() == 1 );
411 Eigen::Array<double, -1, -1> noise( psdF.rows(), psdF.cols() );
417 for(
int k = 0; k < psdFilterTrials; ++k )
419 for(
int cc = 0; cc < psd.cols(); ++cc )
421 for(
int rr = 0; rr < psd.rows(); ++rr )
423 noise( rr, cc ) = normVar;
428 avgRms += noise.square().sum();
431 avgRms = sqrt( avgRms / ( psd.rows() * psd.cols() ) / psdFilterTrials );
433 REQUIRE_THAT( avgRms, Catch::Matchers::WithinAbs( sqrt( 2.2 ), psdFilterTol * sqrt( 2.2 ) ) );
436 GIVEN(
"a rank 3 psd" )
438 WHEN(
"k-alpha=-2.5, f-alph=-2.5, dk Nyquist matched to array size, df Nyquist matched to array size, var=1" )
440 mx::sigproc::psdFilter<double, 3> psdF;
442 Eigen::Array<double, -1, -1> k, psdk;
443 std::vector<double> f, f2s, psd2s;
451 psdk.resize( k.rows(), k.cols() );
452 for(
int cc = 0; cc < psdk.cols(); ++cc )
454 for(
int rr = 0; rr < psdk.rows(); ++rr )
456 if( k( rr, cc ) == 0 )
459 psdk( rr, cc ) = pow( k( rr, cc ), -2.5 );
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 );
472 psd.resize( k.rows(), k.cols(), f2s.size() );
474 for(
int cc = 0; cc < psd.cols(); ++cc )
476 for(
int rr = 0; rr < psd.rows(); ++rr )
478 if( k( rr, cc ) == 0 )
479 psd.
pixel( rr, cc ).setZero();
482 double p = psdk( rr, cc );
485 for(
int pp = 0; pp < psd.planes(); ++pp )
486 psd.
image( pp )( rr, cc ) = psd2s[pp];
491 double dk = k( 0, 1 ) - k( 0, 0 );
492 double df = f[1] - f[0];
494 int rv = psdF.psd( psd, dk, dk, df );
497 REQUIRE( psdF.rows() == psd.rows() );
498 REQUIRE( psdF.cols() == psd.cols() );
499 REQUIRE( psdF.planes() == psd.planes() );
507 for(
int k = 0; k < psdFilterTrials; ++k )
509 for(
int pp = 0; pp < psd.planes(); ++pp )
511 for(
int cc = 0; cc < psd.cols(); ++cc )
513 for(
int rr = 0; rr < psd.rows(); ++rr )
515 noise.
image( pp )( rr, cc ) = normVar;
521 for(
int pp = 0; pp < noise.planes(); ++pp )
522 avgRms += noise.
image( pp ).square().sum();
525 avgRms = sqrt( avgRms / ( psd.rows() * psd.cols() * psd.planes() ) / psdFilterTrials );
527 REQUIRE_THAT( avgRms, Catch::Matchers::WithinAbs( 1.0, psdFilterTol ) );
529 WHEN(
"k-alpha=-3.5, f-alph=-1.5, dk arb, df arb, var=2" )
531 mx::sigproc::psdFilter<double, 3> psdF;
533 Eigen::Array<double, -1, -1> k, psdk;
534 std::vector<double> f, f2s, psd2s;
542 psdk.resize( k.rows(), k.cols() );
543 for(
int cc = 0; cc < psdk.cols(); ++cc )
545 for(
int rr = 0; rr < psdk.rows(); ++rr )
547 if( k( rr, cc ) == 0 )
550 psdk( rr, cc ) = pow( k( rr, cc ), -3.5 );
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 );
563 psd.resize( k.rows(), k.cols(), f2s.size() );
565 for(
int cc = 0; cc < psd.cols(); ++cc )
567 for(
int rr = 0; rr < psd.rows(); ++rr )
569 if( k( rr, cc ) == 0 )
570 psd.
pixel( rr, cc ).setZero();
573 double p = psdk( rr, cc );
576 for(
int pp = 0; pp < psd.planes(); ++pp )
577 psd.
image( pp )( rr, cc ) = psd2s[pp];
582 double dk = k( 0, 1 ) - k( 0, 0 );
583 double df = f[1] - f[0];
585 int rv = psdF.psd( psd, dk, dk, df );
588 REQUIRE( psdF.rows() == psd.rows() );
589 REQUIRE( psdF.cols() == psd.cols() );
590 REQUIRE( psdF.planes() == psd.planes() );
598 for(
int k = 0; k < psdFilterTrials; ++k )
600 for(
int pp = 0; pp < psd.planes(); ++pp )
602 for(
int cc = 0; cc < psd.cols(); ++cc )
604 for(
int rr = 0; rr < psd.rows(); ++rr )
606 noise.
image( pp )( rr, cc ) = normVar;
612 for(
int pp = 0; pp < noise.planes(); ++pp )
613 avgRms += noise.
image( pp ).square().sum();
616 avgRms = sqrt( avgRms / ( psd.rows() * psd.cols() * psd.planes() ) / psdFilterTrials );
618 REQUIRE_THAT( avgRms, Catch::Matchers::WithinAbs( sqrt( 2.0 ), psdFilterTol * sqrt( 2 ) ) );