ivs.sigproc.filtering

1 """ 2 Filter non-equidistant signals (spectra, timeseries...). 3 4 Currently implemented: 5 - Gaussian 6 - Sinc 7 - Box 8 - Iterative Nonlinear filter 9 10 Example usage: 11 12 Generate some data: 13 14 >>> x = np.linspace(0,150,10000) 15 >>> x_= np.linspace(0,150,100) 16 >>> y = np.sin(2*pi/20.*x)+np.random.normal(size=len(x)) 17 18 Apply some filter: 19 20 >>> x1,y1,pnts = filter_signal(x,y,"gauss",sigma=1) 21 >>> x2,y2,pnts = filter_signal(x,y,"pijpers",delta=1) 22 >>> x3,y3,pnts = filter_signal(x,y,"box",window_width=10) 23 >>> x3_,y3_,pnts = filter_signal(x,y,"box",window_width=10,x_template=x_,threads=2) 24 25 >>> import pylab as pl 26 >>> p = pl.plot(x,y,'ko') 27 >>> p = pl.plot(x1,y1,'ro',mec='r',label='Gauss') 28 >>> p = pl.plot(x2,y2,'bo',mec='b',label='Pijpers') 29 >>> p = pl.plot(x3,y3,'go',mec='g',label='Box') 30 >>> p = pl.plot(x3_,y3_,'gs',label='Box with template') 31 >>> p = pl.legend() 32 33 """ 34 import copy 35 import logging 36 import numpy as np 37 from numpy import sqrt,exp,pi,cos,sinc,trapz,average 38 39 from scipy.integrate import trapz 40 from multiprocessing import Process, Manager,cpu_count 41 42 from ivs.aux.decorators import make_parallel 43 from ivs.timeseries.decorators import parallel_pergram,defaults_filtering 44 45 logger = logging.getLogger("TS.FILTERING")

46 47 #{ Main filter program for convolution-based filters 48 49 @defaults_filtering 50 @parallel_pergram 51 @make_parallel 52 -def filter_signal(x,y,ftype,f0=None,fn=None,step=1,x_template=None,**kwargs):

53 """ 54 Filter a signal. 55 56 Handles equidistant and non-equidistant data. C{x} coordinate array 57 must be monotonically increasing, but does not need to be strict monotonically 58 increasing (e.g. their can be duplicate C{x} coordinates). 59 60 It is possible to define a new template C{x} coordinate array on which to 61 filter to signal C{y}. 62 63 For some filters, it is possible to set an adjustable window: that is, 64 you can give an array to the C{window} parameter, but it needs to be the 65 same length as C{x_template}. 66 67 Example usage: 68 69 Generate some data: 70 71 >>> import time 72 >>> x = np.linspace(0,150,10000) 73 >>> y = np.sin(2*np.pi/20.*x)+np.random.normal(size=len(x)) 74 75 Apply some filter: 76 77 >>> c0 = time.time() 78 >>> x1,y1,pnts = filter_signal(x,y,"gauss",sigma=1) 79 >>> print time.time()-c0 80 >>> c0 = time.time() 81 >>> x1,y1,pnts = filter_signal(x,y,"gauss",sigma=1,threads=2) 82 >>> print time.time()-c0 83 84 Output parameter C{pnts} shows you how many points of the original signal 85 were used to compute the given point. E.g. for a box filter, it will give 86 the number of points over which was averaged. 87 88 @param ftype: one of 'gauss','pijpers','box','inl' 89 @type ftype: string 90 @rtype: tuple 91 @return: output from the used filter: typically (x, y, pnts) 92 """ 93 if x_template is None: 94 x_template = x + 0. 95 if f0 is None: f0 = 0 96 if fn is None: fn = len(x_template) 97 #-- set window and kernel function 98 ftype = ftype.lower() 99 window = globals()[ftype+'_window'] 100 kernel = globals()[ftype+'_kernel'] 101 logger.info("Applying filter: %s"%(ftype)) 102 103 #-- set window width 104 lower_window,higher_window = window(**kwargs) 105 106 #-- perhaps the window is a number, perhaps an array. Make sure it is an 107 # array 108 if not isinstance(lower_window,np.ndarray): 109 lower_window = lower_window*np.ones(len(x_template))[f0:fn:step] 110 if not isinstance(higher_window,np.ndarray): 111 higher_window = higher_window*np.ones(len(x_template))[f0:fn:step] 112 113 #-- start running through timeseries (make searchsorted local for speedup) 114 logger.debug("FILTER between index %d-%d with step %d"%(f0,fn,step)) 115 x_template = x_template[f0:fn:step] 116 searchsorted = x.searchsorted 117 out = [kernel(x,y,t,index=searchsorted([t-1e-30-lower_window[i],t+1e-30+higher_window[i]]), 118 **kwargs) for i,t in enumerate(x_template)] 119 #-- that's it! 120 return tuple([x_template] + list(np.array(out).T))

121

122 #} 123 124 #{ Basic functions for convolution-based filters 125 126 -def _fixed_window(window_width=0):

127 """ 128 Define general window 129 @rtype: (float,float) 130 @return: window limits 131 """ 132 lower_window = window_width/2. 133 higher_window = window_width/2. 134 #logger.debug("Window: width=%f"%(window_width)) 135 return lower_window,higher_window

136 #}

137 #{ Convolution-based filter kernels and Windows 138 139 -def gauss_window(sigma=1.,limit=4.):

140 """ 141 Define Gaussian_window 142 @rtype: (float,float) 143 @return: window limits 144 """ 145 if isinstance(sigma,np.ndarray): 146 sigma = sigma.max() 147 return _fixed_window(window_width=2*sigma*limit)

148

149 150 -def gauss_kernel(times, signal, t,sigma=1.,index=(0,-1),norm_weights=True):

151 """ 152 Define Gaussian kernel. 153 154 #Import necessary modules: 155 #>>> from pylab import plot,figure,title,subplot 156 #>>> from sigproc.timeseries import tfreq 157 158 #Generate some data: 159 #>>> x = linspace(0,150,10000) 160 #>>> y = random.normal(size=len(x)) 161 162 #Apply Gaussian filter twice, subtract and calculate periodogram: 163 #>>> y1,pnts = filter_signal(x,y,"gauss",sigma=1) 164 #>>> y2,pnts = filter_signal(x,y,"gauss",sigma=3) 165 166 #Plot the results: 167 #>>> p=figure(figsize=(8,18)) 168 #>>> p=subplot(311);p=title("GAUSSIAN FILTER") 169 #>>> p=plot(x,y,'ko') 170 #>>> p=plot(x,y1,'r-',linewidth=2) 171 #>>> p=plot(x,y2,'b-',linewidth=2) 172 #>>> p=subplot(312) 173 #>>> p=plot(x,y1-y2,'ko') 174 #>>> p=subplot(313) 175 #>>> p=plot(per1[0],per1[1],'k-') 176 #>>> p=plot(per2[0],per2[1],'r-') 177 178 @rtype: (ndarray,) 179 """ 180 index0,indexn = index 181 #-- maybe sigma is changing 182 if isinstance(sigma,np.ndarray): 183 sigma = sigma[(indexn+index0)/2] 184 #-- compute convolution kernel -- should we normalize? so that sum(weights)=1 185 times_ = times[index0:indexn] 186 signal_ = signal[index0:indexn] 187 weights = 1./(sqrt(2.*pi)*sigma) * exp( -(times_-t)**2./(2.*sigma**2.)) 188 if norm_weights and len(times)>1: 189 convolved = trapz(weights*signal_,x=times_-t) 190 norm_fact = trapz(weights,x=times_-t) 191 else: 192 convolved = average(signal_,weights=weights) 193 norm_fact = 1 194 convolved_signal = convolved/norm_fact 195 return convolved_signal,len(times_)

196

197 -def mad(a, c=0.6745):

198 """ 199 Median Absolute Deviation of a 1D array: 200 201 median(abs(a - median(a))) / c 202 203 @rtype: float 204 """ 205 return np.median(np.abs(a-np.median(a)))/c

206

207 -def inl_window(**kwargs):

208 """ 209 Defines the INL window 210 211 @rtype: (float,float) 212 @return: window limits 213 """ 214 return _fixed_window(**kwargs)

215

216 -def inl_kernel(times_in,signal_in,t,c=0.6745,sig_level=3.,tolerance=0.01,index=(0,-1),window_width=None):

217 """ 218 Iterative Nonlinear Filter (Aigrain). 219 220 Example usage: 221 222 #Import necessary modules: 223 #>>> from pylab import plot,figure,title,subplot 224 #>>> from sigproc.timeseries import tfreq 225 226 #Generate some data: 227 #>>> x = linspace(0,150,10000) 228 #>>> y = sin(2*pi*x/20.*x)+random.normal(size=len(x)) 229 230 #Apply INL filter twice, subtract and calculate periodogram: 231 #>>> y1,pnts1,mads = filter_signal(x,y,"inl",window_width=1) 232 #>>> y2,pnts2,mads = filter_signal(x,y,"inl",window_width=10) 233 #>>> freq1,per1,stat1 = tfreq.scargle(x,y,norm='amplitude',fn=4) 234 #>>> freq2,per2,stat2 = tfreq.scargle(x,y1-y2,norm='amplitude',fn=4) 235 236 #Plot the results: 237 #>>> p=figure(figsize=(8,18)) 238 #>>> p=subplot(311);p=title("INL FILTER") 239 #>>> p=plot(x,y,'ko') 240 #>>> p=plot(x,y1,'r-',linewidth=2) 241 #>>> p=plot(x,y2,'b-',linewidth=2) 242 #>>> p=subplot(312) 243 #>>> p=plot(x,y1-y2,'ko') 244 #>>> p=subplot(313) 245 #>>> p=plot(per1[0],per1[1],'k-') 246 #>>> p=plot(per2[0],per2[1],'r-') 247 248 @rtype: (ndarray,ndarray) 249 @return: continuum,sigma 250 """ 251 index0,indexn = index 252 times = times_in[index0:indexn] 253 signal = signal_in[index0:indexn] 254 sigma = mad(signal,c=c) 255 continuum = np.median(signal) 256 outliers = (np.abs(continuum-signal)>sig_level*sigma) 257 max_iter = 3 258 for iteration in range(max_iter): 259 signal_ = signal[outliers==0] 260 sigma = mad(signal_,c=c) 261 continuum_ = np.median(signal_) 262 if abs((continuum_-continuum)/continuum_)<tolerance: 263 break 264 else: 265 continuum = continuum_ 266 outliers = (np.abs(continuum - signal_)>sig_level*sigma) 267 268 return continuum_, sigma,len(times)

269

270 -def pijpers_window(delta=1.):

271 """ 272 Defines the window for the Pijpers filter. 273 274 @rtype: (float,float) 275 @return: window limits 276 """ 277 return _fixed_window(window_width=100*delta) 278 return _fixed_window(**kwargs)

279

280 281 -def pijpers_kernel(times_in, signal_in,t,limit=0.000001,delta=1.,sigma=1.,gamma=0, 282 r=0.0001,index0=0,indexn=-1,index=(0,-1),norm_weights=True):

283 """ 284 Defines the Pijpers (2006) filter kernel. 285 286 Equals box-multiplying in the frequency domain. 287 288 Example usage: 289 290 #Import necessary modules: 291 #>>> from pylab import plot,figure,title,subplot 292 #>>> from sigproc.timeseries import tfreq 293 294 #Generate some data: 295 #>>> x = linspace(0,150,10000) 296 #>>> y = random.normal(size=len(x)) 297 298 #Apply Gaussian filter twice, subtract and calculate periodogram: 299 #>>> y1,pnts1 = filter_signal(x,y,"pijpers",delta=1) 300 #>>> y2,pnts2 = filter_signal(x,y,"pijpers",delta=3) 301 #>>> freq1,per1,stat1 = tfreq.scargle(x,y,norm='amplitude',fn=4) 302 #>>> freq2,per2,stat2 = tfreq.scargle(x,y1-y2,norm='amplitude',fn=4) 303 304 #Plot the results: 305 #>>> p=figure(figsize=(8,18)) 306 #>>> p=subplot(311);p=title("PIJPERS FILTER") 307 #>>> p=plot(x,y,'ko') 308 #>>> p=plot(x,y1,'r-',linewidth=2) 309 #>>> p=plot(x,y2,'b-',linewidth=2) 310 #>>> p=subplot(312) 311 #>>> p=plot(x,y1-y2,'ko') 312 #>>> p=subplot(313) 313 #>>> p=plot(per1[0],per1[1],'k-') 314 #>>> p=plot(per2[0],per2[1],'r-') 315 316 @rtype: (ndarray,) 317 @return: convolved signal 318 """ 319 index0,indexn = index 320 times = times_in[index0:indexn] 321 signal = signal_in[index0:indexn] 322 323 delta = delta*2*pi 324 325 x = delta * (times-t) 326 #-- compute convolution kernel -- should we normalize? so that sum(weights)=1 327 weights = delta/pi * sinc(x/pi) * exp(-1/4. * r**2 * x**2) * cos(gamma*x) 328 #-- compensate for severe unequidistant sampling 329 convolved = np.trapz(weights*signal,x=times-t) 330 norm_fact = np.trapz(weights,x=times-t) 331 convolved_signal = convolved/norm_fact 332 return convolved_signal,len(times)

333

334 -def box_window(**kwargs):

335 """ 336 Defines the window for the box-filter. 337 338 @rtype: (float,float) 339 @return: window limits 340 """ 341 return _fixed_window(**kwargs)

342

343 -def box_kernel(times_in,signal_in,t,window_width=None,index=(0,-1),norm_weights=True):

344 """ 345 Defines the Box filter kernel (moving average). 346 347 Equals sinc-multiplication of frequency domain. 348 349 Example usage: 350 351 #Import necessary modules: 352 #>>> from pylab import plot,figure,title,subplot 353 #>>> from sigproc.timeseries import tfreq 354 355 #Generate some data: 356 #>>> x = linspace(0,150,10000) 357 #>>> y = random.normal(size=len(x)) 358 359 #Apply Gaussian filter twice, subtract and calculate periodogram: 360 #>>> y1,pnts = filter_signal(x,y,"box",window_width=1) 361 #>>> y2,pnts = filter_signal(x,y,"box",window_width=3) 362 363 @rtype: (ndarray,) 364 @return: convolved signal 365 """ 366 index0,indexn = index 367 times = times_in[index0:indexn] 368 signal = signal_in[index0:indexn] 369 weights = np.ones(len(times)) 370 if norm_weights: 371 weights /= np.sum(weights) 372 373 convolved_signal = np.sum(weights*signal) 374 return convolved_signal,len(times)

375

376 377 378 379 #def box_kernel(times_in,signal_in,t,window_width=None,index=(0,-1), 380 #norm_weights=True,is_error_array=False): 381 #""" 382 #Defines the Box filter kernel (moving average). 383 384 #Equals sinc-multiplication in frequency domain. 385 386 #When C{is_error_array} is C{True}, the integration will be done differently. 387 #Instead of the usual 388 389 #F_j = sum_i (F_i dx_{i,j}) / sum_i (dx_{i,j}) 390 391 #it will be 392 393 #s_j = sqrt(sum_i (s_i dx_{i,j})**2) / sum_i (dx_{i,j}) 394 395 #You can only set C{is_error_array} to C{True} when C{norm_weights=True}. 396 397 #Example usage: 398 399 ##Import necessary modules: 400 ##>>> from pylab import plot,figure,title,subplot 401 ##>>> from sigproc.timeseries import tfreq 402 403 ##Generate some data: 404 ##>>> x = linspace(0,150,10000) 405 ##>>> y = random.normal(size=len(x)) 406 407 ##Apply Gaussian filter twice, subtract and calculate periodogram: 408 ##>>> y1,pnts = filter_signal(x,y,"box",window_width=1) 409 ##>>> y2,pnts = filter_signal(x,y,"box",window_width=3) 410 411 #@rtype: (ndarray,) 412 #@return: convolved signal 413 #""" 414 #index0,indexn = index 415 #times = times_in[index0:indexn] 416 #signal = signal_in[index0:indexn] 417 #weights = np.ones(len(times)) 418 #if norm_weights: 419 #if is_error_array: 420 #convolved = sqrt(trapz((weights*signal)**2,x=times-t))/np.trapz(weights,x=times-t) 421 #else: 422 #convolved = trapz( weights*signal, x=times-t) /np.trapz(weights,x=times-t) 423 #norm_fact = 1 424 ##trapz(weights*signal,x=times-t) 425 ##norm_fact = trapz(weights,x=times-t) 426 #else: 427 #assert(is_error_array==False) 428 #convolved = average(signal,weights=weights) 429 #norm_fact = 1 430 431 #convolved_signal = convolved/norm_fact 432 ##weights /= np.sum(weights) 433 ##convolved_signal = np.sum(weights*signal) 434 #return convolved_signal,len(times) 435 436 437 438 439 440 441 442 443 444 #} 445 446 -def test():

447 """ 448 449 >>> from pylab import show 450 >>> p=show() 451 452 """ 453 import doctest 454 doctest.testmod()

455 456 457 if __name__ == "__main__": 458 test() 459

Source Code for Module ivs.sigproc.filtering