ivs.timeseries.freqanalyse

1 # -*- coding: utf-8 -*- 2 """ 3 Frequency Analysis Routines. 4 5 Here are some examples: 6 7 Section 1. Pulsation frequency analysis 8 ======================================= 9 10 Do a frequency analysis of the star HD129929, after Aerts 2004. We reproduce 11 their Fig. 8 here. 12 13 Import necessary modules: 14 15 >>> from ivs.catalogs import vizier 16 17 Read in the data and do an iterative prewhitening frequency analysis: 18 19 >>> data,units,comms = vizier.search('J/A+A/415/241/table1') 20 >>> params = iterative_prewhitening(data.HJD,data.Umag-data.Umag.mean(),f0=6.2,fn=7.2,maxiter=6,\ 21 stopcrit=(stopcrit_scargle_snr,4.,6,)) 22 >>> print pl.mlab.rec2txt(params,precision=6) 23 const ampl freq phase e_const e_ampl e_freq e_phase stopcrit 24 0.000310 0.014722 6.461700 0.443450 0.000401 0.001323 0.000006 0.089865 4.950496 25 0.000000 0.014866 6.978306 -0.050189 0.000000 0.001224 0.000006 0.082305 5.677022 26 0.000000 0.011687 6.449591 0.016647 0.000000 0.001090 0.000007 0.093280 5.747565 27 0.000000 0.011546 6.990431 -0.482814 0.000000 0.001001 0.000006 0.086678 6.454955 28 0.000000 0.009380 6.590938 -0.382048 0.000000 0.000938 0.000007 0.100016 6.510729 29 0.000000 0.007645 6.966174 0.323627 0.000000 0.000863 0.000008 0.112876 5.703801 30 31 Plot the results: 32 33 >>> p = pl.figure() 34 >>> p = pl.vlines(params['freq'],0,params['ampl'],color='k',linestyle='-') 35 >>> p = pl.xlabel('Frequency [c d$^{-1}$]') 36 >>> p = pl.ylabel('Amplitude [magnitude]') 37 >>> p = pl.annotate('$\ell=2$',(6.46,0.015),xycoords='data',va='bottom',ha='center',fontsize='large') 38 >>> p = pl.annotate('$\ell=0$',(6.59,0.010),xycoords='data',va='bottom',ha='center',fontsize='large') 39 >>> p = pl.annotate('$\ell=1$',(6.99,0.015),xycoords='data',va='bottom',ha='center',fontsize='large') 40 >>> p = pl.xlim(6.2,7.2) 41 >>> p = pl.ylim(0,0.018) 42 43 44 Section 2: Line profile variability 45 =================================== 46 47 We generate some line profiles as Gaussian profiles, with varying average. This 48 as an analogue for purely radial variability. 49 50 >>> wave = np.linspace(4950.,5000.,200.) 51 >>> times = np.linspace(0,1.,100) 52 >>> A,sigma,mu = 0.5, 5., 4975. 53 >>> frequency = 1/0.2 54 >>> wshifts = 10.*np.sin(2*np.pi*frequency*times) 55 >>> spectra = np.zeros((len(times),len(wave))) 56 >>> for i,wshift in enumerate(wshifts): 57 ... spectra[i] = 1.-evaluate.gauss(wave,[A,mu+wshift,sigma])+np.random.normal(size=len(wave),scale=0.01) 58 59 Once the line profiles are constructed, we can compute the Fourier transform 60 of these lines: 61 62 >>> output = spectrum_2D(times,wave,spectra,f0=frequency/2.,fn=3*frequency,df=0.001,threads=2,method='scargle',subs_av=True,full_output=True) 63 64 With this output, we can find retrieve the frequency. We plot the periodograms 65 for each wavelength bin on the left, and on the right the average over all 66 wavelength bins: 67 68 >>> p = pl.figure() 69 >>> p = pl.subplot(121) 70 >>> p = pl.imshow(output['pergram'][1][::-1],extent=[wave[0],wave[-1],frequency/2,3*frequency],aspect='auto') 71 >>> p = pl.xlabel('Wavelength (A)') 72 >>> p = pl.ylabel('Frequency (c/d)') 73 >>> cbar = pl.colorbar() 74 >>> cbar.set_label('Amplitude') 75 >>> p = pl.subplot(122) 76 >>> p = pl.plot(output['pergram'][1].mean(axis=1),output['pergram'][0],'k-') 77 >>> p = pl.ylim(frequency/2,3*frequency) 78 >>> p = pl.xlabel('Amplitude') 79 >>> p = pl.ylabel('Frequency (c/d)') 80 81 ]]include figure]]timeseries_freqanalyse_02.png] 82 83 We can then fix the frequency and compute all the parameters. However, we now 84 choose not to subtract the average profile, but instead use the GLS periodogram 85 to include the constant 86 87 >>> output = spectrum_2D(times,wave,spectra,f0=frequency-0.05,fn=frequency+0.05,df=0.001,threads=2,method='gls',subs_av=False,full_output=True) 88 89 >>> p = pl.figure() 90 >>> p = pl.subplot(221) 91 >>> p = pl.errorbar(wave,output['pars']['const'],yerr=output['pars']['e_const'],fmt='ko-') 92 >>> p,q = pl.xlabel('Wavelength (A)'),pl.ylabel('Constant') 93 >>> p = pl.subplot(222) 94 >>> p = pl.errorbar(wave,output['pars']['ampl'],yerr=output['pars']['e_ampl'],fmt='ro-') 95 >>> p,q = pl.xlabel('Wavelength (A)'),pl.ylabel('Amplitude') 96 >>> p = pl.subplot(223) 97 >>> p = pl.errorbar(wave,output['pars']['freq'],yerr=output['pars']['e_freq'],fmt='ko-') 98 >>> p,q = pl.xlabel('Wavelength (A)'),pl.ylabel('Frequency (c/d)') 99 >>> p = pl.subplot(224) 100 >>> p = pl.errorbar(wave,output['pars']['phase'],yerr=output['pars']['e_phase'],fmt='ro-') 101 >>> p,q = pl.xlabel('Wavelength (A)'),pl.ylabel('Phase') 102 103 ]]include figure]]timeseries_freqanalyse_03.png] 104 105 Section 3. Time-frequency analysis 106 ================================== 107 108 Example usage: we generate some data where the frequency jumps to double 109 its value in the middle of the time series. The first half is thus given by 110 111 >>> params = np.rec.fromarrays([[10.],[1.],[0.],[0.]],names=['freq','ampl','const','phase']) 112 >>> times = np.linspace(0,15,1000) 113 >>> signal = evaluate.sine(times,params) 114 115 And the second half by 116 117 >>> params['freq'] = 20. 118 >>> signal[:len(times)/2] = evaluate.sine(times[:len(times)/2],params) 119 120 The time-frequency analysis is calculate via the command 121 122 >>> output = time_frequency(times,signal,fn=30.) 123 124 And the outcome can be plotted via 125 126 >>> p = pl.figure() 127 >>> p = pl.imshow(output['pergram'][1].T[::-1],aspect='auto',extent=[times[0],times[-1],output['pergram'][0][0],output['pergram'][0][-1]]) 128 >>> p,q = pl.xlabel('Time (d)'),pl.ylabel('Frequency (c/d)') 129 130 ]]include figure]]timeseries_freqanalyse_04.png] 131 132 or 133 134 >>> p = pl.figure() 135 >>> p = pl.subplot(221) 136 >>> p = pl.errorbar(output['times'],output['pars']['const'],yerr=output['pars']['e_const'],fmt='ko-') 137 >>> p,q = pl.xlabel('Time (d)'),pl.ylabel('Constant') 138 >>> p = pl.subplot(222) 139 >>> p = pl.errorbar(output['times'],output['pars']['ampl'],yerr=output['pars']['e_ampl'],fmt='ro-') 140 >>> p,q = pl.xlabel('Time (d)'),pl.ylabel('Amplitude') 141 >>> p = pl.subplot(223) 142 >>> p = pl.errorbar(output['times'],output['pars']['freq'],yerr=output['pars']['e_freq'],fmt='ko-') 143 >>> p,q = pl.xlabel('Time (d)'),pl.ylabel('Frequency (c/d)') 144 >>> p = pl.subplot(224) 145 >>> p = pl.errorbar(output['times'],output['pars']['phase'],yerr=output['pars']['e_phase'],fmt='ro-') 146 >>> p,q = pl.xlabel('Time (d)'),pl.ylabel('Phase') 147 148 ]]include figure]]timeseries_freqanalyse_05.png] 149 150 Author: Pieter Degroote 151 """ 152 import logging 153 import numpy as np 154 import pylab as pl 155 from ivs.sigproc import fit 156 from ivs.sigproc import evaluate 157 from ivs.timeseries import pergrams 158 from ivs.timeseries.decorators import defaults_pergram 159 from ivs.aux import numpy_ext 160 import os 161 162 logger = logging.getLogger("TS.FREQANAL")

163 164 #{ Convenience functions 165 166 -def find_frequency(times,signal,method='scargle',model='sine',full_output=False, 167 optimize=0,max_loops=20, scale_region=0.1, scale_df=0.20, model_kwargs=None, 168 correlation_correction=True,prewhiteningorder_snr=False, 169 prewhiteningorder_snr_window=1.,**kwargs):

170 """ 171 Find one frequency, automatically going to maximum precision and return 172 parameters & error estimates. 173 174 This routine makes the frequency grid finer until it is finer than the 175 estimated error on the frequency. After that, it will compute (harmonic) 176 parameters and estimate errors. 177 178 There is a possibility to escape this optimization by setting scale_df=0 or 179 freqregscale=0. 180 181 You can include a nonlinear least square update of the parameters, by 182 setting the keyword C{optimize=1} (optimization outside loop) or 183 C{optimize=2} (optimization after each iteration). 184 185 The method with which to find the frequency can be set with the keyword 186 C{method}, the model used to fit and optimize should be set with C{model}. 187 Extra keywords for the model functions should go in C{model_kwargs}. If 188 C{method} is a tuple, the first method will be used for the first frequency 189 search only. This could be useful to take advantage of such methods as 190 fasper which do not allow for iterative zoom-ins. By default, the function looks 191 for the highest (or deepest in the case of the pdm method) peak, but instead it is 192 possible to go for the peak having the highest SNR before prewhitening by setting 193 C{prewhiteningorder_snr} to True. In this case, the noise spectrum is calculated 194 using a convolution with a C{prewhiteningorder_snr_window} wide box. 195 196 Possible extra keywords: see definition of the used periodogram function. 197 198 B{Warning}: the timeseries must be B{sorted in time} and B{cannot contain 199 the same timepoint twice}. Otherwise, a 'ValueError, concatenation problem' 200 can occur. 201 202 Example keywords: 203 - 'correlation_correction', default=True 204 - 'freqregscale', default=0.5: factor for zooming in on frequency 205 - 'dfscale', default = 0.25: factor for optimizing frequency resolution 206 207 Example usage: We generate a sine signal 208 209 >>> times = np.linspace(0,100,1000) 210 >>> signal = np.sin(2*np.pi*2.5*times) + np.random.normal(size=len(times)) 211 212 Compute the frequency 213 214 >>> parameters, pgram, model = find_frequency(times,signal,full_output=True) 215 216 Make a plot: 217 218 >>> p = pl.figure() 219 >>> p = pl.axes([0.1,0.3,0.85,0.65]) 220 >>> p = pl.plot(pgram[0],pgram[1],'k-') 221 >>> p = pl.xlim(2.2,2.8) 222 >>> p = pl.ylabel('Amplitude') 223 >>> p = pl.axes([0.1,0.1,0.85,0.2]) 224 >>> p = pl.plot(pgram[0][:-1],np.diff(pgram[0]),'k-') 225 >>> p = pl.xlim(2.2,2.8) 226 >>> p,q = pl.xlabel('Frequency (c/d)'),pl.ylabel('Frequency resolution $\Delta F$') 227 228 ]]include figure]]timeseries_freqanalyse_06.png] 229 230 @rtype: record array(, 2x1Darray, 1Darray) 231 @return: parameters and errors(, periodogram, model function) 232 """ 233 if model_kwargs is None: 234 model_kwargs = dict() 235 #-- initial values 236 e_f = 0 237 freq_diff = np.inf 238 prev_freq = -np.inf 239 counter = 0 240 241 f_max = np.inf 242 f_min = 0.#-np.inf 243 244 #-- calculate periodogram until frequency precision is 245 # under 1/10th of correlation corrected version of frequency error 246 method_kwargs = kwargs.copy() # don't modify the dictionary the user gave 247 248 while freq_diff>e_f/10.: 249 #-- possibly, we might want to use different periodograms for the first 250 # calculation than for the zoom in, since some periodograms are faster 251 # than others but do not have the ability to 'zoom in' (e.g. the FFT) 252 if freq_diff==np.inf and not isinstance(method,str): 253 method_ = method[1] 254 method = method[0] # override method to be a string the next time 255 #-- calculate periodogram 256 freqs,ampls = getattr(pergrams,method)(times,signal,**method_kwargs) 257 f0,fn,df = freqs[0],freqs[-1],freqs[1]-freqs[0] 258 #-- now use the second method for the zoom-ins from now on 259 if freq_diff==np.inf and not isinstance(method,str): 260 method = method_ 261 #-- extract the frequency: this part should be generalized, but for now, 262 # it will do: 263 if method in ['pdm']: 264 frequency = freqs[np.argmin(ampls)] 265 #-- instead of going for the highest peak, let's get the most significant one 266 if prewhiteningorder_snr: 267 if counter == 0: #we calculate a noise spectrum with a convolution in a 1 d-1 window 268 windowlength = float(prewhiteningorder_snr_window)/(freqs[1]-freqs[0]) 269 window = np.ones(int(windowlength))/float(windowlength) 270 ampls_ = np.concatenate((ampls[::-1],ampls,ampls[::-1])) #we mirror the amplitude spectrum on both ends so the convolution will be better near the edges 271 noises_ = np.convolve(ampls_, window, 'same') 272 noises = np.split(noises_,3)[1] #and we recut the resulted convolution to match the original frequency range 273 freqs_old = np.copy(freqs) 274 noises_old = np.copy(noises) 275 else: 276 noises = np.interp(freqs,freqs_old,noises_old) #we use the original noise spectrum in this narrower windows too, which should save some time, and avoid the problem of having a wider window for the SNR calculation than the width of the zoom-in window 277 frequency = freqs[np.argmax(ampls/noises)] 278 else: 279 frequency = freqs[np.argmax(ampls)] 280 if full_output and counter==0: 281 freqs_,ampls_ = freqs,ampls 282 #-- estimate parameters and calculate a fit, errors and residuals 283 params = getattr(fit,model)(times,signal,frequency,**model_kwargs) 284 if hasattr(fit,'e_'+model): 285 errors = getattr(fit,'e_'+model)(times,signal,params,correlation_correction=correlation_correction) 286 e_f = errors['e_freq'][-1] 287 #-- possibly there are not errors defined for this fitting functions 288 else: 289 errors = None 290 #-- optimize inside loop if necessary and if we gained prediction 291 # value: 292 if optimize==2: 293 params_,errors_,gain = fit.optimize(times,signal,params,model) 294 if gain>0: 295 params = params_ 296 logger.info('Accepted optimization (gained %g%%)'%gain) 297 298 #-- improve precision 299 freq_diff = abs(frequency-prev_freq) 300 prev_freq = frequency 301 freq_region = fn-f0 302 f0 = max(f_min,frequency-freq_region*scale_region/2.) 303 fn = min(f_max,frequency+freq_region*scale_region/2.) 304 df *= scale_df 305 method_kwargs['f0'] = f0 306 method_kwargs['fn'] = fn 307 method_kwargs['df'] = df 308 #-- possibilities to escape iterative zoom in 309 #print '---> {counter}/{max_loops}: freq={frequency} ({f0}-->{fn}/{df}), e_f={e_f}, freq_diff={freq_diff}'.format(**locals()),max(ampls) 310 if scale_region==0 or scale_df==0: 311 break 312 if counter >= max_loops: 313 logger.error("Frequency precision not reached in %d steps, breaking loop"%(max_loops)) 314 break 315 if (fn-f0)/df<5: 316 logger.error("Frequency precision not reached with stepsize %e , breaking loop"%(df/scale_df)) 317 break 318 counter += 1 319 #-- optimize parameters outside of loop if necessary: 320 if optimize==1: 321 params_,errors_,gain = fit.optimize(times,signal,params,model) 322 if gain>0: 323 params = params_ 324 logger.info('Accepted optimization (gained %g%%)'%gain) 325 #-- add the errors to the parameter array if possible 326 if errors is not None: 327 params = numpy_ext.recarr_join(params,errors) 328 logger.info("%s model parameters via %s periodogram:\n"%(model,method)+pl.mlab.rec2txt(params,precision=8)) 329 #-- when full output is required, return parameters, periodogram and fitting 330 # function 331 if full_output: 332 mymodel = getattr(evaluate,model)(times,params) 333 return params,(freqs_,ampls_),mymodel 334 else: 335 return params

336

337 338 339 -def iterative_prewhitening(times,signal,maxiter=1000,optimize=0,method='scargle', 340 model='sine',full_output=False,stopcrit=None,correlation_correction=True, 341 prewhiteningorder_snr=False,prewhiteningorder_snr_window=1.,**kwargs):

342 """ 343 Fit one or more functions to a timeseries via iterative prewhitening. 344 345 This function will use C{find_frequency} to fit the function parameters to 346 the original signal (including any previously found parameters), 347 optimize the parameters if needed, and remove it from the data to search 348 for a new frequency in the residuals. 349 350 It is always the original signal that is used to fit all parameters again; 351 B{only the (optimized) frequency is remembered from step to step} (Vanicek's 352 method). 353 354 You best set C{maxiter} to some sensable value, to hard-limit the number of 355 frequencies that will be searched for. You can additionally use a C{stopcrit} 356 and stop looking for frequencies once it is reached. C{stopcrit} should be 357 a tuple; the first argument is the function to call, the other arguments 358 are passed to the function, after the mandatory arguments C{times,signal, 359 modelfunc,allparams,pergram}. See L{stopcrit_scargle_snr} for an example of 360 such a function. 361 362 By default, the function looks for the highest (or deepest in the case of the pdm 363 method) peak, but instead it is possible to go for the peak having the highest 364 SNR before prewhitening by setting C{prewhiteningorder_snr} to True. In this case, 365 the noise spectrum is calculated using a convolution with a 366 C{prewhiteningorder_snr_window} wide box. Usage of this is strongly encouraged, 367 especially combined with L{stopcrit_scargle_snr} as C{stopcrit}. 368 369 @return: parameters, model(, model function) 370 @rtype: rec array(, ndarray) 371 """ 372 residuals = signal.copy() 373 frequencies = [] 374 stop_criteria = [] 375 while maxiter: 376 #-- compute the next frequency from the residuals 377 params,pergram,this_fit = find_frequency(times,residuals,method=method, 378 full_output=True,correlation_correction=correlation_correction, 379 prewhiteningorder_snr=prewhiteningorder_snr, 380 prewhiteningorder_snr_window=prewhiteningorder_snr_window,**kwargs) 381 382 #-- do the fit including all frequencies 383 frequencies.append(params['freq'][-1]) 384 allparams = getattr(fit,model)(times,signal,frequencies) 385 386 #-- if there's a need to optimize, optimize the last n parameters 387 if optimize>0: 388 residuals_for_optimization = residuals 389 if optimize<=len(params): 390 model_fixed_params = getattr(evaluate,model)(times,allparams[:-optimize]) 391 residuals_for_optimization -= model_fixed_params 392 uparams,e_uparams, gain = fit.optimize(times,residuals_for_optimization,allparams[-optimize:],model) 393 #-- only accept the optimization if we gained prediction power 394 if gain>0: 395 allparams[-optimize:] = uparams 396 logger.info('Accepted optimization (gained %g%%)'%gain) 397 398 #-- compute the residuals to use in the next prewhitening step 399 modelfunc = getattr(evaluate,model)(times,allparams) 400 residuals = signal - modelfunc 401 402 #-- exhaust the counter 403 maxiter -= 1 404 405 #-- check stop criterion 406 if stopcrit is not None: 407 func = stopcrit[0] 408 args = stopcrit[1:] 409 condition,value = func(times,signal,modelfunc,allparams,pergram,*args) 410 logger.info('Stop criterion (%s): %.3g'%(func.__name__,value)) 411 stop_criteria.append(value) 412 if condition: 413 logger.info('Stop criterion reached') 414 break 415 416 #-- calculate the errors 417 e_allparams = getattr(fit,'e_'+model)(times,signal,allparams,correlation_correction=correlation_correction) 418 419 allparams = numpy_ext.recarr_join(allparams,e_allparams) 420 if stopcrit is not None: 421 allparams = numpy_ext.recarr_join(allparams,np.rec.fromarrays([stop_criteria],names=['stopcrit'])) 422 423 if full_output: 424 return allparams,modelfunc 425 else: 426 return allparams

427

428 429 430 431 -def spectrum_2D(x,y,matrix,weights_2d=None,show_progress=False, 432 subs_av=True,full_output=False,**kwargs):

433 """ 434 Compute a 2D periodogram. 435 436 Rows (first axis) should be spectra chronologically 437 438 x are time points (length N) 439 y are second axis units (e.g. wavelengths) (length NxM) 440 441 If the periodogram/wavelength combination has a large range, this script 442 can produce a B{ValueError} or B{MemoryError}. To solve this, you could 443 iterate this function over a subset of wavelength bins yourself, and write 444 the results to a file. 445 446 This function also outputs a model, which you can use to subtract from the 447 data (probably together with the average profile, which is also returned) 448 and look further for any frequencies:: 449 450 output = spectrum_2D(x,y,matrix,subs_av=True) 451 new_matrix = matrix - output['avprof'] - output['model' 452 output2 = spectrum_2D(x,y,new_matrix,subs_av=False) 453 454 If you want to prewhiten a model, you'd probably want to use the same 455 frequency across the whole line profile. You can hack this by setting 456 C{f0=frequency} and C{fn=frequency+df} with C{df} the size of the frequency 457 step. 458 459 B{Example usage}: first we generate some variable line profiles. In this case, 460 this is just a radial velocity variation 461 462 >>> times = np.linspace(0,150,100) 463 >>> wavel = np.r_[4500:4600:1.0] 464 >>> matrix = np.zeros((len(times),len(wavel))) 465 >>> for i in xrange(len(times)): 466 ... central_wave = 5*np.sin(2*np.pi/10*times[i]) 467 ... matrix[i,:] = 1 - 0.5*np.exp( -(wavel-4550-central_wave)**2/10**2) 468 469 Once the line profiles are constructed, we can compute the Fourier transform 470 of these lines: 471 472 >>> output = spectrum_2D(times,wavel,matrix,method='scargle',model='sine',f0=0.05,fn=0.3,subs_av=True,full_output=True) 473 474 With this output, we can find retrieve the frequency. We plot the periodograms 475 for each wavelength bin on the left, and on the right the average over all 476 wavelength bins: 477 478 >>> p = pl.figure() 479 >>> p = pl.subplot(121) 480 >>> p = pl.imshow(output['pergram'][1][::-1],extent=[wavel[0],wavel[-1],0.05,0.3],aspect='auto') 481 >>> p = pl.subplot(122) 482 >>> p = pl.plot(output['pergram'][1].mean(axis=1),output['pergram'][0],'k-') 483 >>> p = pl.ylim(0.05,0.3) 484 485 We can then fix the frequency and compute all the parameters. However, we now 486 choose not to subtract the average profile, but instead use the GLS periodogram 487 to include the constant 488 489 >>> output = spectrum_2D(times,wavel,matrix,method='gls',model='sine',f0=0.095,fn=0.105,full_output=True) 490 491 >>> p = pl.figure() 492 >>> p = pl.subplot(221) 493 >>> p = pl.errorbar(wavel,output['pars']['const'],yerr=output['pars']['e_const'],fmt='ko-') 494 >>> p = pl.subplot(222) 495 >>> p = pl.errorbar(wavel,output['pars']['ampl'],yerr=output['pars']['e_ampl'],fmt='ro-') 496 >>> p = pl.subplot(223) 497 >>> p = pl.errorbar(wavel,output['pars']['freq'],yerr=output['pars']['e_freq'],fmt='ko-') 498 >>> p = pl.subplot(224) 499 >>> p = pl.errorbar(wavel,output['pars']['phase'],yerr=output['pars']['e_phase'],fmt='ro-') 500 501 @return: dict with keys C{avprof} (2D array), C{pars} (rec array), C{model} (1D array), C{pergram} (freqs,2Darray) 502 @rtype: dict 503 """ 504 #-- compute average profile 505 if subs_av: 506 matrix_av = np.outer(np.ones(len(matrix)),matrix.mean(axis=0)) 507 matrix = matrix - matrix_av 508 else: 509 matrix_av = 0. 510 511 #-- prepare output of sine-parameters 512 params = [] 513 freq_spectrum = [] 514 mymodel = [] 515 #-- do frequency analysis 516 for iwave in xrange(len(matrix[0])): 517 signal = matrix[:,iwave] 518 if weights_2d is not None: 519 weights = weights_2d[:,iwave] 520 kwargs['weights'] = weights 521 522 #-- make sure output is always a tuple, in case full output was asked 523 # we don't want iterative zoom in so set scale_df=0 524 out = find_frequency(x,signal,full_output=full_output,scale_df=0,**kwargs) 525 526 #-- add the parameters of this wavelength bin to the list 527 if full_output: 528 params.append(out[0]) 529 freq_spectrum.append(out[1][1]) 530 mymodel.append(out[2]) 531 else: 532 params.append(out) 533 534 #-- prepare output 535 output = {} 536 output['avprof'] = matrix_av 537 output['pars'] = np.hstack(params) 538 if full_output: 539 output['pergram'] = out[1][0],np.vstack(freq_spectrum).T 540 output['model'] = np.vstack(mymodel).T 541 542 return output

543

544 @defaults_pergram 545 -def time_frequency(times,signal,window_width=None,n_windows=100, 546 window='rectangular',detrend=None,**kwargs):

547 """ 548 Short Time (Fourier) Transform. 549 550 Slide a window through the timeseries, multiply the timeseries with the 551 window, and perform a Fourier Transform. 552 553 It is best to fix explicitly C{f0}, C{fn} and C{df}, to limit the computation 554 time! 555 556 Extra kwargs go to L{find_frequency} 557 558 @param n_windows: number of slices 559 @type n_windows: integer 560 @param window_width: width of each slice (defaults to T/20) 561 @type window_width: float 562 @param detrend: detrending function, accepting times and signal as args 563 @type detrend: callable 564 @param window: window function to apply 565 @type window: string 566 @return: spectrogram, times used, parameters, errors, points used per slice 567 @rtype: dict 568 """ 569 if window_width is None: 570 window_width = kwargs.get('window_width',times.ptp()/20.) 571 #-- cut light curve until window fits exactly 572 stft_times = np.linspace(times[0]+window_width/2.,times[-1]-window_width/2.,n_windows) 573 574 #-- get f0,fn and df to fix in all computations. If they are not given, they 575 # are set to default values by the decorator 576 f0 = kwargs.pop('f0') 577 fn = kwargs.pop('fn') 578 df = kwargs.pop('df') 579 nyq_stat = kwargs.pop('nyq_stat',fn) 580 581 #-- prepare arrays for parameters, points and spectrum 582 pars = [] 583 pnts = np.zeros(n_windows) 584 spec = None 585 #-- compute the periodogram for each slice 586 for i,t in enumerate(stft_times): 587 region = (abs(times-t) <= (window_width/2.)) 588 times_ = times[region] 589 signal_ = signal[region] 590 if detrend: 591 times_,signal_ = detrend(times_,signal_) 592 pnts[i] = len(times_) 593 if len(times_)>1: 594 output = find_frequency(times_,signal_,full_output=True,f0=f0,fn=fn,df=df,nyq_stat=nyq_stat,scale_df=0,**kwargs) 595 if spec is None: 596 spec = np.ones((n_windows,len(output[1][1]))) 597 pars.append(output[0]) 598 spec[i,:len(output[1][1])] = output[1][1] 599 else: 600 nanpars = np.rec.array(np.nan*np.ones(len(pars[-1].dtype.names)),dtype=pars[-1].dtype) 601 pars.append(nanpars) 602 #try: 603 # pars.append(np.nan*pars[-1]) 604 #except: 605 # print pars 606 # print pars[-1] 607 # print pars[-1].dtype 608 # print 0*pars[-1] 609 # print np.nan*pars[-1] 610 # raise 611 spec[i] = np.nan 612 out = {} 613 out['times'] = stft_times 614 out['pars'] = np.hstack(pars) 615 out['pergram'] = (output[1][0],spec) 616 return out

617

618 #{ Convenience stop-criteria 619 620 -def stopcrit_scargle_prob(times,signal,modelfunc,allparams,pergram,crit_value):

621 """ 622 Stop criterium based on probability. 623 """ 624 value = pergrams.scargle_probability(pergram[1].max(),times,pergram[0]) 625 print value 626 return value>crit_value,value

627

628 -def stopcrit_scargle_snr(times,signal,modelfunc,allparams,pergram,crit_value,width=6.):

629 """ 630 Stop criterium based on signal-to-noise ratio. 631 """ 632 width = width/2. 633 argmax = np.argmax(pergram[1]) 634 ampls = pergram[1] 635 start = max(0,pergram[0][argmax]-width) 636 stop = min(pergram[0][argmax]+width,pergram[0][-1]) 637 if start==0: 638 stop += width-pergram[0][argmax] 639 if stop==pergram[0][-1]: 640 start = pergram[0][-1]-pergram[0][argmax]+width 641 ampls = ampls[(start<=pergram[0]) & (pergram[0]<=stop)] 642 value = pergram[1][argmax]/ampls.mean() 643 return value<crit_value,value

644 #}

645 646 -def autocorrelation(frequencies,power,max_step=1.5,interval=(),threshold=None,method=1):

647 """ 648 Compute the autocorrelation. 649 650 The formulate to estimate the autocorrelation in the periodogram is 651 652 R(k) = 1 / (n * s^2) * \sum_{t=1}^n (X_t - mu) * (X_{t+k} - mu) 653 654 where n is the number of points in the interval, mu an estimator for the 655 sample mean and s is an estimator for the standard deviation of the sample. 656 657 When computed over a large interval, we have to take into account that 658 we cannot sample more points than we have in the spectrum, so we have to 659 average out over fewer and fewer points, the more we shift the spectrum. 660 661 Above formula becomes 662 663 R(k) = 1 / ((n-k) * s^2) * \sum_{t=1}^{n-k} (X_t - mu) * (X_{t+k} - mu) 664 665 Highs and lows are cut of. The lower cut off value is the sample mean, the 666 higher cut off is a user-defined multiple of the sample mean. 667 668 This function can also be used to compute period spacings: just invert the 669 frequency spacing, reverse the order and reinterpolate to be equidistant: 670 Example: 671 >>> periods = linspace(1./p[0][-1],1./p[0][0],2*len(p[0])) 672 >>> ampl = interpol.linear_interpolation(1./p[0][::-1],p[1][::-1],periods) 673 >>> ampl = where(isnan(ampl),hstack([ampl[1:],ampl[:1]]),ampl) 674 675 @param frequencies: frequency array 676 @type frequencies: numpy 1d array 677 @param power: power spectrum 678 @type power: numpy 1d array 679 @keyword max_step: maximum time shift 680 @type max_step: float 681 @keyword interval: tuple of frequencies (start, end) 682 @type interval: tuple of floats 683 @keyword threshold: high cut off value (in units of sample mean) 684 @type threshold: float 685 @return: domain of autocorrelation and autocorrelation 686 @rtype: (ndarray,ndarray) 687 """ 688 #-- cut out the interesting part of the spectrum 689 if interval is not (): 690 cut_out = frequencies[(interval[0]<=frequencies) & (frequencies<=interval[1])] 691 start_freq = cut_out[0] 692 stop_freq = cut_out[-1] 693 start = np.argmin(abs(frequencies-interval[0])) 694 stop = np.argmin(abs(frequencies-interval[1])) 695 else: 696 start = 1 697 stop = len(frequencies)-1 698 699 #-- compute the frequency step 700 Dfreq = (frequencies[start+1] - frequencies[start+0]) 701 max_step = int(max_step/Dfreq) 702 autocorr = [] 703 variance = [] 704 mean = np.average(power) 705 706 #-- cut of high peaks at a signal to noise level of 6 707 # cut of low values at a signal to noise level of 1 708 if threshold is not None: 709 power[power>=(threshold*mean)] = threshold*mean 710 #power = where(less(power, mean), mean, power) 711 712 #-- normalize power as to arrive at the AUTO-correlation function. 713 mean = np.average(power) 714 power = power-mean 715 716 #-- compute autocorrelation. If nessecary, take border effects into 717 # account. 718 for i in xrange(2,max_step): 719 end_s = min(len(power), stop+i) 720 end_o = start + end_s - (start+i) 721 original = power[start :end_o] 722 shifted = power[start+i:end_s] 723 if len(original) < 10: 724 logger.error("AUTOCORR: too few points left in interval, breaking up.") 725 break 726 if method==1: 727 autocorr.append(np.average(original*shifted)) 728 else: 729 autocorr.append(np.correlate(original,shifted)) 730 variance.append(np.average(original*original)) 731 732 domain = np.arange(2,max_step) * Dfreq 733 domain = domain[0:len(autocorr)] 734 735 #-- normalize 736 autocorr = np.array(autocorr)/np.array(variance) 737 logger.info("Computed autocorrelation in interval %s with maxstep %s"%(interval,max_step)) 738 return domain, autocorr

739 740 if __name__=="__main__": 741 import doctest 742 import pylab as pl 743 import sys 744 from ivs.aux import argkwargparser 745 from ivs.inout import ascii 746 747 #-- if no arguments are given, we just do a test run 748 if not sys.argv[1:]: 749 doctest.testmod() 750 pl.show() 751 sys.exit() 752 753 #-- if arguments are given, we assume the user wants to run one of the 754 # functions with arguments given in the command line 755 # EXAMPLES: 756 # $:> python freqanalyse.py find_frequency infile=test.dat full_output=True 757 # $:> python freqanalyse.py time_frequency infile=test.dat full_output=True 758 else: 759 method,args,kwargs = argkwargparser.parse() 760 print "Running method %s with arguments %s and keyword arguments %s"%(method,args,kwargs) 761 if '--help' in args or 'help' in args or 'help' in kwargs: 762 sys.exit() 763 full_output = kwargs.get('full_output',False) 764 times,signal = ascii.read2array(kwargs.pop('infile')).T[:2] 765 out = globals()[method](times,signal, **kwargs) 766 767 #-- when find_frequency is called 768 if method=='find_frequency' and full_output: 769 print pl.mlab.rec2txt(out[0],precision=8) 770 pl.figure() 771 pl.subplot(211) 772 pl.plot(out[1][0],out[1][1],'k-') 773 pl.subplot(212) 774 pl.plot(times,signal,'ko',ms=2) 775 pl.plot(times,out[2],'r-',lw=2) 776 pl.show() 777 elif method=='find_frequency': 778 print pl.mlab.rec2txt(out) 779 780 #-- when time_frequency is called 781 elif method=='time_frequency': 782 print pl.mlab.rec2txt(out['pars'],precision=8) 783 pl.figure() 784 pl.imshow(out['pergram'][1].T[::-1],aspect='auto',extent=[out['times'][0],out['times'][-1],out['pergram'][0][0],out['pergram'][0][-1]]) 785 786 787 pl.show() 788

Source Code for Module ivs.timeseries.freqanalyse