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