1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
|
#include "utils.h"
#include "fitting.h" // for unit test
#include "complexdata.h" // for unit test
#include <tjutils/tjtest.h>
Data<float,1> unwrap_phase(const Data<float,1>& phase, int startindex) {
Log<OdinData> odinlog("","unwrap_phase");
int n=phase.size();
Data<float,1> result(n);
if(startindex<0 || startindex>=n) {
ODINLOG(odinlog,errorLog) << "startindex=" << startindex << " out of range (0," << (n-1) << ")" << STD_endl;
return result;
}
for(int j=0; j<n; j++) {
if(phase(j)<float(-PII) || phase(j)>float(PII)) {
ODINLOG(odinlog,errorLog) << "phase(" << j << "/" << n << ")=" << phase(j) << " out of range (" << -PII << "," << PII << ")" << STD_endl;
return result;
}
}
result(startindex)=phase(startindex);
// unwrap going to higher indices
int i=startindex;
int mod=0;
while( (i+1) < n ) {
float diff=phase(i+1)-phase(i);
if(diff>PII) mod--;
if(diff<-PII) mod++;
result(i+1)=phase(i+1)+mod*2.0*PII;
i++;
}
// unwrap going to lower indices
i=startindex;
mod=0;
while( (i-1) >= 0 ) {
float diff=phase(i-1)-phase(i);
if(diff>PII) mod--;
if(diff<-PII) mod++;
result(i-1)=phase(i-1)+mod*2.0*PII;
i--;
}
return result;
}
//////////////////////////////////////////////////
// Unit Test
#ifndef NO_UNIT_TEST
class DataUtilsTest : public UnitTest {
public:
DataUtilsTest() : UnitTest("DataUtils") {}
private:
bool check() const {
Log<UnitTest> odinlog(this,"check");
PolynomialFunction<3> poly;
poly.a[0].val=-20.0;
poly.a[1].val=-10.0;
poly.a[2].val=2.0;
poly.a[3].val=0.5;
int nx=1000;
float xrange=10.0;
Data<float,1> xvals(nx);
for(int i=0; i<nx; i++) xvals(i)=(float(i)/float(nx)-0.5)*xrange;
Data<float,1> yvals(poly.get_function(xvals));
// yvals.autowrite("yvals.asc");
ComplexData<1> cplxdata(expc(float2imag(yvals)));
Data<float,1> pha(phase(cplxdata));
// pha.autowrite("pha.asc");
int nstart=5;
int startindex[nstart];
startindex[0]=0;
startindex[1]=nx/3;
startindex[2]=nx/2;
startindex[3]=3*nx/4;
startindex[4]=nx-1;
for(int istart=0; istart<nstart; istart++) {
Data<float,1> pha_unwrap(unwrap_phase(pha,startindex[istart]));
pha_unwrap-=pha_unwrap(startindex[istart])-yvals(startindex[istart]);
// pha_unwrap.autowrite("pha_unwrap.asc");
float diff=sum(fabs(pha_unwrap-yvals))/nx;
if(diff>1.0e-5) {
ODINLOG(odinlog,errorLog) << "unwrap_phase(...," << startindex[istart] << "), diff=" << diff << STD_endl;
return false;
}
}
return true;
}
};
void alloc_DataUtilsTest() {new DataUtilsTest();} // create test instance
#endif
//////////////////////////////////////////////////
// Test Instantiations
#ifdef ODIN_DEBUG
template Array<float,1> matrix_product(const Array<float,2>&, const Array<float,1>&);
template Array<float,1> vector_product(const Array<float,1>&, const Array<float,1>&);
template bool same_shape(const Array<float,1>&, const Array<float,1>&,const TinyVector<int,1>&);
template bool check_range(float&, float, float);
#endif
|
