Package mvpa :: Package misc :: Module signal
[hide private]
[frames] | no frames]

Source Code for Module mvpa.misc.signal

  1  #emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*- 
  2  #ex: set sts=4 ts=4 sw=4 et: 
  3  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 
  4  # 
  5  #   See COPYING file distributed along with the PyMVPA package for the 
  6  #   copyright and license terms. 
  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 # perform scipy detrend 59 bp = 0 # no break points by default 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 # perform regression-based detrend 68 return __detrend_regress(data, perchunk=perchunk, 69 polyord=polyord, opt_reg=opt_reg) 70 else: 71 # raise exception because not found 72 raise ValueError('Specified model type (%s) is unknown.' 73 % (model))
74
75 -def __detrend_regress(data, perchunk=True, polyord=None, opt_reg=None):
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 # Process the polyord to be a list with length of the number of chunks 105 if not polyord is None: 106 if not isSequenceType(polyord): 107 # repeat to be proper length 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 # loop over chunks if necessary 115 if perchunk: 116 # get the unique chunks 117 uchunks = data.uniquechunks 118 119 # loop over each chunk 120 reg = [] 121 for n, chunk in enumerate(uchunks): 122 # get the indices for that chunk 123 cinds = data.chunks == chunk 124 125 # see if add in polyord values 126 if not polyord is None: 127 # create the timespan 128 x = N.linspace(-1, 1, cinds.sum()) 129 # create each polyord with the value for that chunk 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 # take out mean over entire dataset 136 reg = [] 137 # see if add in polyord values 138 if not polyord is None: 139 # create the timespan 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 # see if add in optional regs 145 if not opt_reg is None: 146 # add in the optional regressors, too 147 reg.append(opt_reg) 148 149 # combine the regs 150 if len(reg) > 0: 151 if len(reg) > 1: 152 regs = N.hstack(reg) 153 else: 154 regs = reg[0] 155 else: 156 # no regs to remove 157 raise ValueError("You must specify at least one " + \ 158 "regressor to regress out.") 159 160 # perform the regression 161 res = lstsq(regs, data.samples) 162 163 # remove all but the residuals 164 yhat = N.dot(regs, res[0]) 165 data.samples -= yhat 166 167 # return the results 168 return res, regs
169