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
|
#include "filter_detrend.h"
#include "fitting.h"
void FilterDeTrend::init(){
nlow=5;
nlow.set_description("Number of low frequency components to be removed");
append_arg(nlow,"nlow");
zeromean=false;
zeromean.set_description("Zero mean of resulting timecourse");
append_arg(zeromean,"zeromean");
}
bool FilterDeTrend::process(Data<float,4>& data, Protocol& prot) const {
Log<Filter> odinlog(c_label(),"process");
Range all=Range::all();
TinyVector<int,4> shape=data.shape();
if(nlow<2 || shape(0)<2) {
ODINLOG(odinlog,warningLog) << "Too few time points: nlow=" << int(nlow) << ", shape=" << shape << STD_endl;
return true;
}
if(nlow==2) { // Remove linear trend
LinearFunction linf;
for(int islice=0; islice<shape(1); islice++) {
for(int iphase=0; iphase<shape(2); iphase++) {
for(int iread=0; iread<shape(3); iread++) {
linf.fit(data(all,islice,iphase,iread));
if(zeromean) {
for(int irep=0; irep<shape(0); irep++) data(irep,islice,iphase,iread)-=(irep*linf.m.val+linf.c.val);
} else {
float meanval=0.5*(shape(0)-1)*linf.m.val;
for(int irep=0; irep<shape(0); irep++) data(irep,islice,iphase,iread)-=(irep*linf.m.val-meanval);
}
}
}
}
} else {
TinyVector<int,4> lowshape=shape;
lowshape(0)=nlow;
Data<float,4> lowdata(data.copy());
lowdata.congrid(lowshape);
lowdata.congrid(shape);
if(!zeromean) {
for(int islice=0; islice<lowshape(1); islice++) {
for(int iphase=0; iphase<lowshape(2); iphase++) {
for(int iread=0; iread<lowshape(3); iread++) {
lowdata(all,islice,iphase,iread)-=mean(lowdata(all,islice,iphase,iread));
}
}
}
}
data=data-lowdata;
}
return true;
}
|
