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