1
2
3
4
5
6
7
8
9 """Simple preprocessing utilities"""
10
11 __docformat__ = 'restructuredtext'
12
13
14 import numpy as N
15 from scipy import signal
16 from scipy.linalg import lstsq
17 from scipy.special import legendre
18
19 from operator import isSequenceType
20
21 from mvpa.misc.support import getBreakPoints
22
23 -def detrend(data, perchunk=False, model='linear',
24 polyord=None, opt_reg=None):
25 """
26 Given a dataset, detrend the data inplace either entirely
27 or per each chunk
28
29 :Parameters:
30 `data` : `Dataset`
31 dataset to operate on
32 `perchunk` : bool
33 either to operate on whole dataset at once or on each chunk
34 separately
35 `model`
36 Type of detrending model to run. If 'linear' or 'constant',
37 scipy.signal.detrend is used to perform a linear or demeaning
38 detrend. If 'regress', then you specify the polyord and opt_reg
39 arguments to define regressors to regress out of the dataset.
40 `polyord` : int or list
41 Order of the Legendre polynomial to remove from the data. This
42 will remove every polynomial up to and including the provided
43 value. For example, 3 will remove 0th, 1st, 2nd, and 3rd order
44 polynomials from the data. N.B.: The 0th polynomial is the
45 baseline shift, the 1st is the linear trend.
46 If you specify a single int and perchunk is True, then this value
47 is used for each chunk. You can also specify a differnt polyord
48 value for each chunk by providing a list or ndarray of polyord
49 values the length of the number of chunks.
50 `opt_reg` : ndarray
51 Optional ndarray of additional information to regress out from the
52 dataset. One example would be to regress out motion parameters.
53 As with the data, time is on the first axis.
54
55 """
56
57 if model in ['linear', 'constant']:
58
59 bp = 0
60
61 if perchunk:
62 bp = getBreakPoints(data.chunks)
63
64 data.samples[:] = signal.detrend(data.samples, axis=0,
65 type=model, bp=bp)
66 elif model in ['regress']:
67
68 return __detrend_regress(data, perchunk=perchunk,
69 polyord=polyord, opt_reg=opt_reg)
70 else:
71
72 raise ValueError('Specified model type (%s) is unknown.'
73 % (model))
74
76 """
77 Given a dataset, perform a detrend inplace, regressing out polynomial
78 terms as well as optional regressors, such as motion parameters.
79
80 :Parameters:
81 `data`: `Dataset`
82 Dataset to operate on
83 `perchunk` : bool
84 Either to operate on whole dataset at once or on each chunk
85 separately. If perchunk is True, all the samples within a
86 chunk should be contiguous and the chunks should be sorted in
87 order from low to high.
88 `polyord` : int
89 Order of the Legendre polynomial to remove from the data. This
90 will remove every polynomial up to and including the provided
91 value. For example, 3 will remove 0th, 1st, 2nd, and 3rd order
92 polynomials from the data. N.B.: The 0th polynomial is the
93 baseline shift, the 1st is the linear trend.
94 If you specify a single int and perchunk is True, then this value
95 is used for each chunk. You can also specify a differnt polyord
96 value for each chunk by providing a list or ndarray of polyord
97 values the length of the number of chunks.
98 `opt_reg` : ndarray
99 Optional ndarray of additional information to regress out from the
100 dataset. One example would be to regress out motion parameters.
101 As with the data, time is on the first axis.
102 """
103
104
105 if not polyord is None:
106 if not isSequenceType(polyord):
107
108 polyord = [polyord]*len(data.uniquechunks)
109 elif perchunk and len(polyord) != len(data.uniquechunks):
110 raise ValueError("If you specify a sequence of polyord values " + \
111 "they sequence length must match the " + \
112 "number of unique chunks in the data.")
113
114
115 if perchunk:
116
117 uchunks = data.uniquechunks
118
119
120 reg = []
121 for n, chunk in enumerate(uchunks):
122
123 cinds = data.chunks == chunk
124
125
126 if not polyord is None:
127
128 x = N.linspace(-1, 1, cinds.sum())
129
130 for n in range(polyord[n] + 1):
131 newreg = N.zeros((data.nsamples, 1))
132 newreg[cinds, 0] = legendre(n)(x)
133 reg.append(newreg)
134 else:
135
136 reg = []
137
138 if not polyord is None:
139
140 x = N.linspace(-1, 1, data.nsamples)
141 for n in range(polyord[0] + 1):
142 reg.append(legendre(n)(x)[:, N.newaxis])
143
144
145 if not opt_reg is None:
146
147 reg.append(opt_reg)
148
149
150 if len(reg) > 0:
151 if len(reg) > 1:
152 regs = N.hstack(reg)
153 else:
154 regs = reg[0]
155 else:
156
157 raise ValueError("You must specify at least one " + \
158 "regressor to regress out.")
159
160
161 res = lstsq(regs, data.samples)
162
163
164 yhat = N.dot(regs, res[0])
165 data.samples -= yhat
166
167
168 return res, regs
169