mxlib
c++ tools for analyzing astronomical data and other tasks by Jared R. Males. [git repo]
Loading...
Searching...
No Matches
psdUtils_test.cpp
Go to the documentation of this file.
1/** \file psdUtils_test.cpp
2 * \brief Tests power-spectral-density utilities.
3 */
4#include "../../catch2/catch.hpp"
5
6#include <vector>
7#include <Eigen/Dense>
8
9#define MX_NO_ERROR_REPORTS
10
11#include "../../../include/sigproc/psdFilter.hpp"
12#include "../../../include/sigproc/psdUtils.hpp"
13#include "../../../include/improc/eigenCube.hpp"
14#include "../../../include/math/randomT.hpp"
15#include "../../../include/math/vectorUtils.hpp"
16
17/** calculating variance from a 1D PSD
18 *
19 * Verify calculations of psdVar1sided, psdVar2sided, and psdVar.
20 *
21 */
22/**
23 * \ingroup psdUtils_unit_tests
24 */
25TEST_CASE( "calculating variance from a 1D PSD", "[sigproc::psdUtils]" )
26{
27 GIVEN( "a 1 sided PSD," )
28 {
29 WHEN( "a flat PSD, midpoint rule" )
30 {
31 std::vector<double> f( 5 ), psd( 5 );
32 for( size_t n = 0; n < f.size(); ++n )
33 {
34 f[n] = n;
35 psd[n] = 1;
36 }
37
38 REQUIRE( mx::sigproc::psdVar( f, psd, 1.0 ) == 5 );
39 }
40 WHEN( "a flat PSD, trapezoid rule by default" )
41 {
42 std::vector<double> f( 5 ), psd( 5 );
43 for( size_t n = 0; n < f.size(); ++n )
44 {
45 f[n] = n;
46 psd[n] = 1;
47 }
48
49 REQUIRE( mx::sigproc::psdVar( f, psd ) == 4 );
50 }
51 }
52 GIVEN( "a 2 sided PSD," )
53 {
54 WHEN( "a flat PSD, midpoint rule" )
55 {
56 std::vector<double> f( 6 ), psd( 6 );
57 f[0] = 0;
58 f[1] = 1;
59 f[2] = 2;
60 f[3] = 3;
61 f[4] = -2;
62 f[5] = -1;
63
64 for( size_t n = 0; n < f.size(); ++n )
65 {
66 psd[n] = 1;
67 }
68 psd[3] = 2; // This one gets twice the power
69
70 REQUIRE( mx::sigproc::psdVar( f, psd, 1.0 ) == 7 );
71 }
72 WHEN( "a flat PSD, trapezoid rule by default" )
73 {
74 std::vector<double> f( 6 ), psd( 6 );
75 f[0] = 0;
76 f[1] = 1;
77 f[2] = 2;
78 f[3] = 3;
79 f[4] = -2;
80 f[5] = -1;
81
82 for( size_t n = 0; n < f.size(); ++n )
83 {
84 psd[n] = 1;
85 }
86 psd[3] = 2; // This one gets twice the power
87
88 REQUIRE( mx::sigproc::psdVar( f, psd ) == 6 );
89 }
90 }
91}
92
93/** Verify scaling and normalization of augment1SidedPSD
94 *
95 */
96/**
97 * \ingroup psdUtils_unit_tests
98 */
99TEST_CASE( "augmenting a 1 sided PSD", "[sigproc::psdUtils]" )
100{
101 GIVEN( "a 1 sided PSD, with a 0 freq value" )
102 {
103 WHEN( "1/f^2" )
104 {
105 std::vector<double> f( 5 ), psd( 5 );
106
107 for( size_t n = 0; n < psd.size(); ++n )
108 {
109 f[n] = n;
110 psd[n] = pow( f[n], -2. );
111 }
112 psd[0] = psd[1];
113
114 mx::sigproc::normPSD( psd, f, 1.0 );
115 // Now have a 1/f^2 PSD with total 1-sided variance of 1.0.
116 REQUIRE_THAT( mx::sigproc::psdVar( f, psd, 1.0 ),
117 Catch::Matchers::WithinAbs( 1.0, 1e-10 ) ); // handles epsilon
118
119 // proceed to augment:
120 std::vector<double> f2s, psd2s;
121 mx::sigproc::augment1SidedPSDFreq( f2s, f );
122 mx::sigproc::augment1SidedPSD( psd2s, psd );
123
124 REQUIRE( f2s.size() == 8 );
125 REQUIRE( psd2s.size() == 8 );
126
127 REQUIRE( f2s[0] == 0 );
128 REQUIRE( f2s[1] == 1 );
129 REQUIRE( f2s[2] == 2 );
130 REQUIRE( f2s[3] == 3 );
131 REQUIRE( f2s[4] == 4 );
132 REQUIRE( f2s[5] == -3 );
133 REQUIRE( f2s[6] == -2 );
134 REQUIRE( f2s[7] == -1 );
135
136 // Should now have 1.0 in bin 0, 0.5 in all other bins.
137 REQUIRE_THAT( psd2s[0], Catch::Matchers::WithinAbs( psd[0], 1e-12 ) );
138 REQUIRE_THAT( psd2s[1], Catch::Matchers::WithinAbs( 0.5 * psd[1], 1e-12 ) );
139 REQUIRE_THAT( psd2s[2], Catch::Matchers::WithinAbs( 0.5 * psd[2], 1e-12 ) );
140 REQUIRE_THAT( psd2s[3], Catch::Matchers::WithinAbs( 0.5 * psd[3], 1e-12 ) );
141 REQUIRE_THAT( psd2s[4], Catch::Matchers::WithinAbs( psd[4], 1e-12 ) );
142 REQUIRE_THAT( psd2s[5], Catch::Matchers::WithinAbs( 0.5 * psd[3], 1e-12 ) );
143 REQUIRE_THAT( psd2s[6], Catch::Matchers::WithinAbs( 0.5 * psd[2], 1e-12 ) );
144 REQUIRE_THAT( psd2s[7], Catch::Matchers::WithinAbs( 0.5 * psd[1], 1e-12 ) );
145
146 // handle machine precision
147 REQUIRE_THAT( mx::sigproc::psdVar( f2s, psd2s, 1.0 ), Catch::Matchers::WithinAbs( 1.0, 1e-10 ) );
148 }
149 }
150}
151
152/** Verify creation of a 1D frequency grid
153 *
154 */
155/**
156 * \ingroup psdUtils_unit_tests
157 */
158TEST_CASE( "creating a 1D frequency grid", "[sigproc::psdUtils]" )
159{
160 GIVEN( "2 sided FFT-order frequency grid" )
161 {
162 WHEN( "dt = 1" )
163 {
164 std::vector<double> f( 10 );
165
166 REQUIRE( mx::sigproc::frequencyGrid( f, 1.0 ) == 0 );
167
168 REQUIRE_THAT( f[0], Catch::Matchers::WithinAbs( 0, 1e-10 ) );
169 REQUIRE_THAT( f[1], Catch::Matchers::WithinAbs( 0.1, 1e-10 ) );
170 REQUIRE_THAT( f[2], Catch::Matchers::WithinAbs( 0.2, 1e-10 ) );
171 REQUIRE_THAT( f[3], Catch::Matchers::WithinAbs( 0.3, 1e-10 ) );
172 REQUIRE_THAT( f[4], Catch::Matchers::WithinAbs( 0.4, 1e-10 ) );
173 REQUIRE_THAT( f[5], Catch::Matchers::WithinAbs( 0.5, 1e-10 ) );
174 REQUIRE_THAT( f[6], Catch::Matchers::WithinAbs( -0.4, 1e-10 ) );
175 REQUIRE_THAT( f[7], Catch::Matchers::WithinAbs( -0.3, 1e-10 ) );
176 REQUIRE_THAT( f[8], Catch::Matchers::WithinAbs( -0.2, 1e-10 ) );
177 REQUIRE_THAT( f[9], Catch::Matchers::WithinAbs( -0.1, 1e-10 ) );
178 }
179
180 WHEN( "dt = 2" )
181 {
182 std::vector<double> f( 10 );
183
184 REQUIRE( mx::sigproc::frequencyGrid( f, 2.5 ) == 0 );
185
186 REQUIRE_THAT( f[0], Catch::Matchers::WithinAbs( 0, 1e-10 ) );
187 REQUIRE_THAT( f[1], Catch::Matchers::WithinAbs( 0.04, 1e-10 ) );
188 REQUIRE_THAT( f[2], Catch::Matchers::WithinAbs( 0.08, 1e-10 ) );
189 REQUIRE_THAT( f[3], Catch::Matchers::WithinAbs( 0.12, 1e-10 ) );
190 REQUIRE_THAT( f[4], Catch::Matchers::WithinAbs( 0.16, 1e-10 ) );
191 REQUIRE_THAT( f[5], Catch::Matchers::WithinAbs( 0.2, 1e-10 ) );
192 REQUIRE_THAT( f[6], Catch::Matchers::WithinAbs( -0.16, 1e-10 ) );
193 REQUIRE_THAT( f[7], Catch::Matchers::WithinAbs( -0.12, 1e-10 ) );
194 REQUIRE_THAT( f[8], Catch::Matchers::WithinAbs( -0.08, 1e-10 ) );
195 REQUIRE_THAT( f[9], Catch::Matchers::WithinAbs( -0.04, 1e-10 ) );
196 }
197 }
198
199 GIVEN( "1 sided frequency grid" )
200 {
201 WHEN( "dt = 1, odd size" )
202 {
203 std::vector<double> f( 5 );
204
205 REQUIRE( mx::sigproc::frequencyGrid( f, 1.0, false ) == 0 );
206
207 REQUIRE_THAT( f[0], Catch::Matchers::WithinAbs( 0.1, 1e-10 ) );
208 REQUIRE_THAT( f[1], Catch::Matchers::WithinAbs( 0.2, 1e-10 ) );
209 REQUIRE_THAT( f[2], Catch::Matchers::WithinAbs( 0.3, 1e-10 ) );
210 REQUIRE_THAT( f[3], Catch::Matchers::WithinAbs( 0.4, 1e-10 ) );
211 REQUIRE_THAT( f[4], Catch::Matchers::WithinAbs( 0.5, 1e-10 ) );
212 }
213
214 WHEN( "dt = 1, even size" )
215 {
216 std::vector<double> f( 6 );
217
218 REQUIRE( mx::sigproc::frequencyGrid( f, 1.0, false ) == 0 );
219
220 REQUIRE_THAT( f[0], Catch::Matchers::WithinAbs( 0.0, 1e-10 ) );
221 REQUIRE_THAT( f[1], Catch::Matchers::WithinAbs( 0.1, 1e-10 ) );
222 REQUIRE_THAT( f[2], Catch::Matchers::WithinAbs( 0.2, 1e-10 ) );
223 REQUIRE_THAT( f[3], Catch::Matchers::WithinAbs( 0.3, 1e-10 ) );
224 REQUIRE_THAT( f[4], Catch::Matchers::WithinAbs( 0.4, 1e-10 ) );
225 REQUIRE_THAT( f[5], Catch::Matchers::WithinAbs( 0.5, 1e-10 ) );
226 }
227 }
228}
eigenCube.hpp
An image cube with an Eigen API.
TEST_CASE
TEST_CASE("calculating variance from a 1D PSD", "[sigproc::psdUtils]")
Definition psdUtils_test.cpp:25
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:827
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:885
mx::sigproc::psdVar
realT psdVar(const std::vector< realT > &f, const std::vector< realT > &PSD, realT half=0.5)
Calculate the variance of a 1-D PSD.
Definition psdUtils.hpp:135
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:448
mx::sigproc::frequencyGrid
int frequencyGrid(std::vector< realT > &vec, realParamT dt, bool fftOrder=true)
Create a 1-D frequency grid.
Definition psdUtils.hpp:258
psdFilter.hpp
Declares and defines a class for filtering with PSDs.
psdUtils.hpp
Tools for working with PSDs.
randomT.hpp
Defines a random number type.
vectorUtils.hpp
Header for the std::vector utilities.