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( fabs( mx::sigproc::psdVar( f, psd, 1.0 ) - 1.0 ) < 1e-10 ); // handles epsilon
112
113 // proceed to augment:
114 std::vector<double> f2s, psd2s;
115 mx::sigproc::augment1SidedPSDFreq( f2s, f );
116 mx::sigproc::augment1SidedPSD( psd2s, psd );
117
118 REQUIRE( f2s.size() == 8 );
119 REQUIRE( psd2s.size() == 8 );
120
121 REQUIRE( f2s[0] == 0 );
122 REQUIRE( f2s[1] == 1 );
123 REQUIRE( f2s[2] == 2 );
124 REQUIRE( f2s[3] == 3 );
125 REQUIRE( f2s[4] == 4 );
126 REQUIRE( f2s[5] == -3 );
127 REQUIRE( f2s[6] == -2 );
128 REQUIRE( f2s[7] == -1 );
129
130 // Should now have 1.0 in bin 0, 0.5 in all other bins.
131 REQUIRE( psd2s[0] == psd[0] );
132 REQUIRE( psd2s[1] == 0.5 * psd[1] );
133 REQUIRE( psd2s[2] == 0.5 * psd[2] );
134 REQUIRE( psd2s[3] == 0.5 * psd[3] );
135 REQUIRE( psd2s[4] == psd[4] );
136 REQUIRE( psd2s[5] == 0.5 * psd[3] );
137 REQUIRE( psd2s[6] == 0.5 * psd[2] );
138 REQUIRE( psd2s[7] == 0.5 * psd[1] );
139
140 // handle machine precision
141 REQUIRE( fabs( mx::sigproc::psdVar( f2s, psd2s, 1.0 ) - 1.0 ) < 1e-10 );
142 }
143 }
144}
145
146/** Verify creation of a 1D frequency grid
147 *
148 * \anchor tests_sigproc_psdUtils_frequencyGrid_1D
149 */
150SCENARIO( "creating a 1D frequency grid", "[sigproc::psdUtils]" )
151{
152 GIVEN( "2 sided FFT-order frequency grid" )
153 {
154 WHEN( "dt = 1" )
155 {
156 std::vector<double> f( 10 );
157
158 REQUIRE( mx::sigproc::frequencyGrid( f, 1.0 ) == 0 );
159
160 REQUIRE( fabs( f[0] - 0 ) < 1e-10 );
161 REQUIRE( fabs( f[1] - 0.1 ) < 1e-10 );
162 REQUIRE( fabs( f[2] - 0.2 ) < 1e-10 );
163 REQUIRE( fabs( f[3] - 0.3 ) < 1e-10 );
164 REQUIRE( fabs( f[4] - 0.4 ) < 1e-10 );
165 REQUIRE( fabs( f[5] - 0.5 ) < 1e-10 );
166 REQUIRE( fabs( f[6] - -0.4 ) < 1e-10 );
167 REQUIRE( fabs( f[7] - -0.3 ) < 1e-10 );
168 REQUIRE( fabs( f[8] - -0.2 ) < 1e-10 );
169 REQUIRE( fabs( f[9] - -0.1 ) < 1e-10 );
170 }
171
172 WHEN( "dt = 2" )
173 {
174 std::vector<double> f( 10 );
175
176 REQUIRE( mx::sigproc::frequencyGrid( f, 2.5 ) == 0 );
177
178 REQUIRE( fabs( f[0] - 0 ) < 1e-10 );
179 REQUIRE( fabs( f[1] - 0.04 ) < 1e-10 );
180 REQUIRE( fabs( f[2] - 0.08 ) < 1e-10 );
181 REQUIRE( fabs( f[3] - 0.12 ) < 1e-10 );
182 REQUIRE( fabs( f[4] - 0.16 ) < 1e-10 );
183 REQUIRE( fabs( f[5] - 0.2 ) < 1e-10 );
184 REQUIRE( fabs( f[6] - -0.16 ) < 1e-10 );
185 REQUIRE( fabs( f[7] - -0.12 ) < 1e-10 );
186 REQUIRE( fabs( f[8] - -0.08 ) < 1e-10 );
187 REQUIRE( fabs( f[9] - -0.04 ) < 1e-10 );
188 }
189 }
190
191 GIVEN( "1 sided frequency grid" )
192 {
193 WHEN( "dt = 1, odd size" )
194 {
195 std::vector<double> f( 5 );
196
197 REQUIRE( mx::sigproc::frequencyGrid( f, 1.0, false ) == 0 );
198
199 REQUIRE( fabs( f[0] - 0.1 ) < 1e-10 );
200 REQUIRE( fabs( f[1] - 0.2 ) < 1e-10 );
201 REQUIRE( fabs( f[2] - 0.3 ) < 1e-10 );
202 REQUIRE( fabs( f[3] - 0.4 ) < 1e-10 );
203 REQUIRE( fabs( f[4] - 0.5 ) < 1e-10 );
204 }
205
206 WHEN( "dt = 1, even size" )
207 {
208 std::vector<double> f( 6 );
209
210 REQUIRE( mx::sigproc::frequencyGrid( f, 1.0, false ) == 0 );
211
212 REQUIRE( fabs( f[0] - 0.0 ) < 1e-10 );
213 REQUIRE( fabs( f[1] - 0.1 ) < 1e-10 );
214 REQUIRE( fabs( f[2] - 0.2 ) < 1e-10 );
215 REQUIRE( fabs( f[3] - 0.3 ) < 1e-10 );
216 REQUIRE( fabs( f[4] - 0.4 ) < 1e-10 );
217 REQUIRE( fabs( f[5] - 0.5 ) < 1e-10 );
218 }
219 }
220}
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::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:136
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
SCENARIO
SCENARIO("calculating variance from a 1D PSD", "[sigproc::psdUtils]")
Definition psdUtils_test.cpp:22