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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
|
#include "complexdata.h"
#include <tjutils/tjtest.h>
#ifdef HAVE_LIBGSL
#include <gsl/gsl_fft_complex.h>
struct GslFftWorkSpace {
int length_cache;
gsl_fft_complex_wavetable* wavetable;
gsl_fft_complex_workspace* workspace;
};
GslFft::GslFft(int length) : ws(new GslFftWorkSpace) {
ws->length_cache=length;
ws->wavetable=gsl_fft_complex_wavetable_alloc(length);
ws->workspace=gsl_fft_complex_workspace_alloc(length);
}
GslFft::~GslFft() {
gsl_fft_complex_wavetable_free(ws->wavetable);
gsl_fft_complex_workspace_free(ws->workspace);
delete ws;
}
void GslFft::fft1d(double* complex_data, bool forward_fft) {
if(forward_fft) gsl_fft_complex_forward (complex_data, 1, ws->length_cache, ws->wavetable, ws->workspace);
else gsl_fft_complex_backward(complex_data, 1, ws->length_cache, ws->wavetable, ws->workspace);
}
#else
#error "GNU Scientific library is missing!"
#endif
//////////////////////////////////////////////////
// Unit Test
#ifndef NO_UNIT_TEST
class ComplexDataTest : public UnitTest {
public:
ComplexDataTest() : UnitTest("ComplexData") {}
private:
bool check() const {
Log<UnitTest> odinlog(this,"check");
int testsize=11; // Use odd number for more thourough testing
ComplexData<2> testarray(testsize,testsize);
TinyVector<int,2> index;
for(unsigned int i=0; i<testarray.numElements(); i++) {
index=testarray.create_index(i);
float radius=norm(index(0)-testsize/2-1,index(1)-testsize/2+2); // off-center disk
if(radius<4.0) testarray(index)=STD_complex(2.0,0.0);
else testarray(index)=STD_complex(0.0,1.0);
}
ComplexData<2> original(testarray);
original.makeUnique();
// test FFT by back and forth transform
testarray.fft(true);
testarray.fft(false);
float diff=sum(cabs(testarray-original));
if(diff>1e-4) {
// Data<float,2>(cabs(testarray)).autowrite("testarray.jdx");
// Data<float,2>(cabs(original)).autowrite("original.jdx");
ODINLOG(odinlog,errorLog) << "FFT test failed, diff=" << diff << STD_endl;
return false;
}
// Test Data::shift against ComplexData::modulate_offset
ComplexData<2> shifttest1(original.shape()); shifttest1=original;
ComplexData<2> shifttest2(original.shape()); shifttest2=original;
int shiftpix=3;
shifttest1.shift(0, shiftpix);
shifttest2.fft(true);
float relshift=float(shiftpix)/float(testsize);
shifttest2.modulate_offset(TinyVector<float,2>(relshift, 0.0));
shifttest2.fft(false);
diff=sum(cabs(shifttest1)-cabs(shifttest2)); // take magnitude difference as modulate_offset adds extra phase
if(diff>0) {
// Data<float,2>(cabs(shifttest1)).autowrite("shifttest1.jdx");
// Data<float,2>(cabs(shifttest2)).autowrite("shifttest2.jdx");
ODINLOG(odinlog,errorLog) << "modulate_offset failed, diff=" << diff << STD_endl;
return false;
}
// Test Data::convert_to complex->complex
ComplexData<3> data3d;
original.convert_to(data3d);
ComplexData<2> convtest;
data3d.convert_to(convtest);
diff=sum(cabs(original-convtest));
if(diff>0.0) {
ODINLOG(odinlog,errorLog) << "convert_to(complex->complex) failed, diff=" << diff << STD_endl;
ODINLOG(odinlog,errorLog) << "original " << original << STD_endl;
ODINLOG(odinlog,errorLog) << "convtest " << convtest << STD_endl;
return false;
}
// Testing convert_to with byte->complex->float
int nbytes=4;
Data<s8bit,1> bytedata(nbytes);
for(int i=0; i<nbytes; i++) bytedata(i)=i;
ComplexData<1> cplxdst;
bytedata.convert_to(cplxdst);
Data<float,1> floatdst;
cplxdst.convert_to(floatdst);
for(int i=0; i<nbytes; i++) {
if(float(bytedata(i))!=floatdst(i)) {
ODINLOG(odinlog,errorLog) << "bytedata=" << bytedata << STD_endl;
ODINLOG(odinlog,errorLog) << "floatdst=" << floatdst << STD_endl;
return false;
}
}
// Testing phasemap()
testsize=1000;
Data<float,1> expected_phase(testsize);
for(int i=0; i<testsize; i++) {
float x=10.0*(float(i)/float(testsize)-0.5); // phase is unwrapped starting at center so we need a function which zero at the center
expected_phase(i)=pow(x,3);
}
ComplexData<1> cplxdata(expc(float2imag(expected_phase)));
Data<float,1> unwrapped_phase(cplxdata.phasemap());
diff=sum(abs(expected_phase-unwrapped_phase));
if(diff>0.03) {
ODINLOG(odinlog,errorLog) << "phasemap failed, diff=" << diff << STD_endl;
return false;
}
return true;
}
};
void alloc_ComplexDataTest() {new ComplexDataTest();} // create test instance
#endif
|
