1
2 """
3 Contains many different periodogram calculations
4
5 Section 1. Basic usage
6 ======================
7
8 Given a time series of the form
9
10 >>> times = np.linspace(0,1,1000)
11 >>> signal = np.sin(times)
12
13 The basic interface is
14
15 >>> freq,ampl = scargle(times,signal)
16
17 If you give no extra information, default values for the start, end and step
18 frequency will be chosen. See the 'defaults_pergram' decorator for a list of
19 keywords common to all periodogram calculations.
20
21 Many periodograms can be computed in parallel, by supplying an extra keyword
22 'threads':
23
24 >>> freq,ampl = scargle(times,signal,threads=2)
25
26 B{Warning}: the timeseries must be B{sorted in time} and B{cannot contain the
27 same timepoint twice}. Otherwise, a 'ValueError, concatenation problem' can
28 occur.
29
30 If something goes wrong in the periodogram computation, be sure to run
31 L{check_input} on your input data. This will print out some basic diagnostics
32 to see if your data are valid.
33
34 Section 2. Nyquist frequency
35 ============================
36
37 The periodogram functions are written such that they never exceed the value
38 of the Nyquist frequency. This behaviour can be changed.
39
40 By default, the Nyquist frequency is defined as half of the inverse of the smallest
41 time step in the data. That is, the C{nyq_stat} variable is set to the function
42 C{np.min}. If you prefer the Nyquist to be defined as the median, set the
43 C{nyq_stat} variable for any periodogram to C{np.median}. If you want a more
44 complex or self-defined function, that is also acceptable. If you give a number
45 to C{nyq_stat}, nothing will be computed but that value will be considered the
46 nyquist frequency.
47
48 Section 3. Speed comparison
49 ===========================
50
51 >>> import time
52
53 The timeseries is generated for N=500,1000,2000,5000,10000,20000,50000
54
55 >>> techniques = ['scargle','deeming','gls','fasper',
56 ... 'schwarzenberg_czerny','pdm','box']
57 >>> Ns = [1000,2000,5000,10000,20000,30000]
58 >>> for tech in techniques[:0]:
59 ... freqs = []
60 ... clock = []
61 ... for N in Ns:
62 ... times = np.linspace(0,1,N)
63 ... signal = np.sin(times) + np.random.normal(size=N)
64 ... c0 = time.time()
65 ... freq,ampl = locals()[tech](times,signal)
66 ... clock.append(time.time()-c0)
67 ... freqs.append(len(freq))
68 ... p = pl.plot(freqs,clock,'o-',label=tech)
69 ... print tech,np.polyfit(np.log10(np.array(freqs)),np.log10(np.array(clock)),1)
70 >>> #p = pl.legend(loc='best',fancybox=True)
71 >>> #q = p.get_frame().set_alpha(0.5)
72 >>> #p,q = pl.xlabel('Number of frequencies'),pl.ylabel('Seconds')
73
74 ]]include figure]]ivs_timeseries_pergram_speeds.png]
75
76 Section 2. Periodogram comparison
77 =================================
78
79 We generate a sinusoidal signal, add some noise and compute the periodogram with
80 the different methods.
81
82 >>> times = np.sort(np.random.uniform(size=500,low=0,high=100))
83 >>> signal = np.sin(2*np.pi/10.*times)
84 >>> signal += np.random.normal(size=500)
85
86 Fourier based periodgrams suited for unequidistant data: L{scargle}, L{deeming},
87 L{fasper} and L{schwarzenberg_czerny}.
88
89 >>> f1,a1 = scargle(times,signal,fn=0.35)
90 >>> f2,a2 = fasper(times,signal,fn=0.35)
91 >>> f3,a3 = deeming(times,signal,fn=0.35)
92 >>> f4,a4 = scargle(times,signal,norm='power',fn=0.35)
93 >>> f5,a5 = schwarzenberg_czerny(times,signal,fn=0.35,nh=1,mode=2)
94
95 >>> p = pl.figure()
96 >>> p = pl.subplot(121)
97 >>> p = pl.plot(f1,a1,'k-',lw=4,label='scargle')
98 >>> p = pl.plot(f2,a2,'r--',lw=4,label='fasper')
99 >>> p = pl.plot(f3,a3,'b-',lw=2,label='deeming')
100 >>> p = pl.xlim(f1[0],f1[-1])
101 >>> leg = pl.legend(fancybox=True,loc='best')
102 >>> leg.get_frame().set_alpha(0.5)
103 >>> p = pl.subplot(122)
104 >>> p = pl.plot(f4,a4,'k-',lw=2,label='scargle')
105 >>> p = pl.plot(f5,a5/2.,'r--',lw=2,label='schwarzenberg-czerny')
106 >>> p = pl.xlim(f1[0],f1[-1])
107 >>> leg = pl.legend(fancybox=True,loc='best')
108 >>> leg.get_frame().set_alpha(0.5)
109
110 ]]include figure]]ivs_timeseries_pergrams_fourier.png]
111
112 Least-square fitting: L{gls} and L{kepler}.
113
114 >>> f1,a1 = gls(times,signal,fn=0.35)
115 >>> f2,a2 = kepler(times,signal,fn=0.35)
116
117 >>> p = pl.figure()
118 >>> p = pl.plot(f1,a1,'k-',lw=2,label='gls')
119 >>> p = pl.plot(f2,a2,'r--',lw=2,label='kepler')
120 >>> p = pl.xlim(f1[0],f1[-1])
121 >>> leg = pl.legend(fancybox=True,loc='best')
122 >>> leg.get_frame().set_alpha(0.5)
123
124 ]]include figure]]ivs_timeseries_pergrams_lsq.png]
125
126 Phase folding techniques: L{box} and L{pdm}
127
128 >>> f1,a1 = box(times,signal,fn=0.35)
129 >>> f2,a2 = pdm(times,signal,fn=0.35)
130
131 >>> p = pl.figure()
132 >>> p = pl.plot(f1,a1,'k-',lw=2,label='box')
133 >>> p = pl.plot(f2,a2,'r--',lw=2,label='pdm')
134 >>> p = pl.xlim(f1[0],f1[-1])
135 >>> leg = pl.legend(fancybox=True,loc='best')
136 >>> leg.get_frame().set_alpha(0.5)
137
138 ]]include figure]]ivs_timeseries_pergrams_phase.png]
139
140 """
141 import logging
142 import numpy as np
143 from numpy import cos,sin,pi
144 from scipy.special import jn
145 from ivs.aux.decorators import make_parallel
146 from ivs.aux import loggers
147 from ivs.aux import termtools
148 from ivs.timeseries.decorators import parallel_pergram,defaults_pergram,getNyquist
149
150 import pyscargle
151 import pyscargle_single
152 import pyfasper
153 import pyfasper_single
154 import pyclean
155 import pyGLS
156 import pyKEP
157 import pydft
158 import multih
159 import deeming as fdeeming
160 import eebls
161
162 logger = logging.getLogger("TS.PERGRAMS")
163
164
165
166
167
168 @defaults_pergram
169 @parallel_pergram
170 @make_parallel
171 -def scargle(times, signal, f0=None, fn=None, df=None, norm='amplitude',
172 weights=None, single=False):
173 """
174 Scargle periodogram of Scargle (1982).
175
176 Several options are available (possibly combined):
177 1. weighted Scargle
178 2. Amplitude spectrum
179 3. Distribution power spectrum
180 4. Traditional power spectrum
181 5. Power density spectrum (see Kjeldsen, 2005 or Carrier, 2010)
182
183 This definition makes use of a Fortran-routine written by Jan Cuypers, Conny
184 Aerts and Peter De Cat. A slightly adapted version is used for the weighted
185 version (adapted by Pieter Degroote).
186
187 Through the option "norm", it's possible to norm the periodogram as to get a
188 periodogram that has a known statistical distribution. Usually, this norm is
189 the variance of the data (NOT of the noise or residuals, see Schwarzenberg-
190 Czerny 1998!).
191
192 Also, it is possible to retrieve the power density spectrum in units of
193 [ampl**2/frequency]. In this routine, the normalisation constant is taken
194 to be the total time span T. Kjeldsen (2005) chooses to multiply the power
195 by the 'effective length of the observing run', which is calculated as the
196 reciprocal of the area under spectral window (in power, and take 2*Nyquist
197 as upper frequency value).
198
199 REMARK: this routine does B{not} automatically remove the average. It is the
200 user's responsibility to do this adequately: e.g. subtract a B{weighted}
201 average if one computes the weighted periodogram!!
202
203 @param times: time points
204 @type times: numpy array
205 @param signal: observations
206 @type signal: numpy array
207 @param weights: weights of the datapoints
208 @type weights: numpy array
209 @param norm: type of normalisation
210 @type norm: str
211 @param f0: start frequency
212 @type f0: float
213 @param fn: stop frequency
214 @type fn: float
215 @param df: step frequency
216 @type df: float
217 @return: frequencies, amplitude spectrum
218 @rtype: array,array
219 """
220 if single: pyscargle_ = pyscargle_single
221 else:
222 pyscargle_ = pyscargle
223
224 sigma=0.;xgem=0.;xvar=0.;n=len(times)
225 T = times.ptp()
226 nf=int((fn-f0)/df+0.001)+1
227 f1=np.zeros(nf,'d');s1=np.zeros(nf,'d')
228 ss=np.zeros(nf,'d');sc=np.zeros(nf,'d');ss2=np.zeros(nf,'d');sc2=np.zeros(nf,'d')
229
230
231 if weights is None:
232 f1,s1=pyscargle_.scar2(signal,times,f0,df,f1,s1,ss,sc,ss2,sc2)
233 else:
234 w=np.array(weights,'float')
235 logger.debug('Weighed scargle')
236 f1,s1=pyscargle_.scar3(signal,times,f0,df,f1,s1,ss,sc,ss2,sc2,w)
237
238
239 if not s1[0]: s1[0]=0.
240 fact = np.sqrt(4./n)
241 if norm =='distribution':
242 s1 /= np.var(signal)
243 elif norm == "amplitude":
244 s1 = fact * np.sqrt(s1)
245 elif norm == "density":
246 s1 = fact**2 * s1 * T
247 return f1, s1
248
249
250
251 @defaults_pergram
252 @parallel_pergram
253 @make_parallel
254 -def fasper(times,signal, f0=None, fn=None, df=None, single=True, norm='amplitude'):
255 """
256 Fasper periodogram from Numerical Recipes.
257
258 Normalisation here is not correct!!
259
260 @param times: time points
261 @type times: numpy array
262 @param signal: observations
263 @type signal: numpy array
264 @param f0: start frequency
265 @type f0: float
266 @param fn: stop frequency
267 @type fn: float
268 @param df: step frequency
269 @type df: float
270 @return: frequencies, amplitude spectrum
271 @rtype: array,array
272 """
273
274 nyq = 1./(2*np.diff(times).mean())
275 mynyq = 1./(2*np.diff(times).min())
276 T = times.ptp()
277 ofac = 1./(df*T)
278 hifac = fn/mynyq*mynyq/nyq
279
280 n = len(times)
281 nout = int(4*ofac*hifac*n*4.)
282 wk1 = np.zeros(nout)
283 wk2 = np.zeros(nout)
284 jmax,prob = 0,0.
285
286 if not single:
287 wk1,wk2,nwk,nout,jmax,prob = pyfasper.fasper(times,signal,ofac,hifac,wk1,wk2,nout,jmax,prob)
288 else:
289 wk1,wk2,nwk,nout,jmax,prob = pyfasper_single.fasper(times,signal,ofac,hifac,wk1,wk2,nout,jmax,prob)
290
291 wk1,wk2 = wk1[:nout],wk2[:nout]*1.5
292 fact = np.sqrt(4./n)
293 if norm =='distribution':
294 wk2 /= np.var(signal)
295 elif norm == "amplitude":
296 wk2 = fact * np.sqrt(wk2)
297 elif norm == "density":
298 wk2 = fact**2 * wk2 * T
299 if f0 is not None:
300 keep = f0<wk1
301 wk1,wk2 = wk1[keep],wk2[keep]
302 return wk1,wk2
303
304
305
306
307
308
309
310
311
312 @defaults_pergram
313 @parallel_pergram
314 @make_parallel
315 -def deeming(times,signal, f0=None, fn=None, df=None, norm='amplitude'):
316 """
317 Deeming periodogram of Deeming et al (1975).
318
319 Thanks to Jan Cuypers
320
321 @param times: time points
322 @type times: numpy array
323 @param signal: observations
324 @type signal: numpy array
325 @param norm: type of normalisation
326 @type norm: str
327 @param f0: start frequency
328 @type f0: float
329 @param fn: stop frequency
330 @type fn: float
331 @param df: step frequency
332 @type df: float
333 @return: frequencies, amplitude spectrum
334 @rtype: array,array
335 """
336
337 nf=int((fn-f0)/df+0.001)+1
338 n = len(times)
339 T = times.ptp()
340 f1,s1 = fdeeming.deeming1(times,signal,f0,df,nf)
341 s1 /= n
342 fact = np.sqrt(4./n)
343 fact = np.sqrt(4./n)
344 if norm =='distribution':
345 s1 /= np.var(signal)
346 elif norm == "amplitude":
347 s1 = fact * np.sqrt(s1)
348 elif norm == "density":
349 s1 = fact**2 * s1 * T
350
351 return f1,s1
352
353
354 @defaults_pergram
355 @parallel_pergram
356 @make_parallel
357 -def gls(times,signal, f0=None, fn=None, df=None, errors=None, wexp=2):
358 """
359 Generalised Least Squares periodogram of Zucher et al (2010).
360
361 @param times: time points
362 @type times: numpy array
363 @param signal: observations
364 @type signal: numpy array
365 @param f0: start frequency
366 @type f0: float
367 @param fn: stop frequency
368 @type fn: float
369 @param df: step frequency
370 @type df: float
371 @return: frequencies, amplitude spectrum
372 @rtype: array,array
373 """
374 T = times.ptp()
375 n = len(times)
376 if errors is None:
377 errors = np.ones(n)
378 maxstep = int((fn-f0)/df+1)
379
380
381 f1 = np.zeros(maxstep)
382 s1 = np.zeros(maxstep)
383 p1 = np.zeros(maxstep)
384 l1 = np.zeros(maxstep)
385
386
387 pyGLS.gls(times+0.,signal+0.,errors,f0,fn,df,wexp,f1,s1,p1,l1)
388 return f1,s1
389
390
391
392
393
394 @defaults_pergram
395 @parallel_pergram
396 @make_parallel
397 -def clean(times,signal, f0=None, fn=None, df=None, freqbins=None, niter=10.,
398 gain=1.0):
399 """
400 Cleaned Fourier periodogram of Roberts et al (1987)
401
402 Parallization probably isn't such a good idea here because of the frequency
403 bins.
404
405 Fortran module probably from John Telting.
406
407 Should always start from zero, so f0 is not an option
408
409 Generate some signal with heavy window side lobes:
410
411 >>> times_ = np.linspace(0,150,1000)
412 >>> times = np.array([times_[i] for i in xrange(len(times_)) if (i%10)>7])
413 >>> signal = np.sin(2*pi/10*times) + np.random.normal(size=len(times))
414
415 Compute the scargle periodogram as a reference, and compare with the
416 the CLEAN periodogram with different gains.
417
418 >>> niter,freqbins = 10,[0,1.2]
419 >>> p1 = scargle(times,signal,fn=1.2,norm='amplitude',threads=2)
420 >>> p2 = clean(times,signal,fn=1.2,gain=1.0,niter=niter,freqbins=freqbins)
421 >>> p3 = clean(times,signal,fn=1.2,gain=0.1,niter=niter,freqbins=freqbins)
422
423 Make a figure of the result:
424
425 >>> p=pl.figure()
426 >>> p=pl.plot(p1[0],p1[1],'k-',label="Scargle")
427 >>> p=pl.plot(p2[0],p2[1],'r-',label="Clean (g=1.0)")
428 >>> p=pl.plot(p3[0],p3[1],'b-',label="Clean (g=0.1)")
429 >>> p=pl.legend()
430
431 ]]include figure]]ivs_timeseries_pergrams_clean.png]
432
433 @keyword freqbins: frequency bins for clean computation
434 @type freqbins: list or array
435 @keyword niter: number of iterations
436 @type niter: integer
437 @keyword gain: gain for clean computation
438 @type gain: float between 0 (no cleaning) and 1 (full cleaning)
439 @return: frequencies, amplitude spectrum
440 @rtype: array,array
441 """
442 T = times.ptp()
443 n = len(times)
444 if freqbins is None:
445 freqbins = [f0,fn]
446
447 startfreqs = np.array(freqbins[0::2])
448 endfreqs = np.array(freqbins[1::2])
449 nbins = len(freqbins)-1
450
451 nf = int(fn/df)
452
453
454 f,wpow,wpha = pyclean.main_clean(times,signal,fn,nf,gain,niter,nbins,\
455 startfreqs,endfreqs)
456
457 return f,wpow
458
459
460
461
462
463
464
465
466
467 @defaults_pergram
468 @parallel_pergram
469 @make_parallel
470 -def schwarzenberg_czerny(times, signal, f0=None, fn=None, df=None, nh=2, mode=1):
471 """
472 Multi harmonic periodogram of Schwarzenberg-Czerny (1996).
473
474 This periodogram follows an F-distribution, so it is possible to perform
475 hypothesis testing.
476
477 If the number of the number of harmonics is 1, then this peridogram reduces
478 to the Lomb-Scargle periodogram except for its better statistic behaviour.
479 This script uses a Fortran procedure written by Schwarzenberg-Czerny.
480
481 Modes:
482 - mode=1: AoV Fisher-Snedecor F(df1,df2) statistic
483 - mode=2: total power fitted in all harmonics for mode=2
484
485 @param times: list of observations times
486 @type times: numpy 1d array
487 @param signal: list of observations
488 @type signal: numpy 1d array
489 @keyword f0: start frequency (cycles per day) (default: 0.)
490 @type f0: float
491 @keyword fn: stop frequency (cycles per day) (default: 10.)
492 @type fn: float
493 @keyword df: step frequency (cycles per day) (default: 0.001)
494 @type df: float
495 @keyword nh: number of harmonics to take into account
496 @type nh: integer
497 @return: frequencies, f-statistic
498 @rtype: array,array
499 """
500 T = times.ptp()
501 n = len(times)
502 frequencies = np.arange(f0, fn+df, df)
503 ll = len(frequencies)
504 th = np.zeros(len(frequencies))
505
506 th = multih.sfou(n,times,signal,ll,f0,df,nh,mode,th)
507
508
509
510 return frequencies,th
511
512 -def DFTpower(time, signal, f0=None, fn=None, df=None, full_output=False):
513
514 """
515 Computes the modulus square of the fourier transform.
516
517 Unit: square of the unit of signal. Time points need not be equidistant.
518 The normalisation is such that a signal A*sin(2*pi*nu_0*t)
519 gives power A^2 at nu=nu_0
520
521 @param time: time points [0..Ntime-1]
522 @type time: ndarray
523 @param signal: signal [0..Ntime-1]
524 @type signal: ndarray
525 @param f0: the power is computed for the frequencies freq = arange(f0,fn,df)
526 @type f0: float
527 @param fn: see f0
528 @type fn: float
529 @param df: see f0
530 @type df: float
531 @return: power spectrum of the signal
532 @rtype: array
533 """
534
535 freqs = np.arange(f0,fn,df)
536 Ntime = len(time)
537 Nfreq = int(np.ceil((fn-f0)/df))
538
539 A = np.exp(1j*2.*pi*f0*time) * signal
540 B = np.exp(1j*2.*pi*df*time)
541 ft = np.zeros(Nfreq, complex)
542 ft[0] = A.sum()
543 for k in range(1,Nfreq):
544 A *= B
545 ft[k] = np.sum(A)
546
547 if full_output:
548 return freqs,ft**2*4.0/Ntime**2
549 else:
550 return freqs,(ft.real**2 + ft.imag**2) * 4.0 / Ntime**2
551
554
555 """
556 Computes the power spectrum of a signal using a discrete Fourier transform.
557
558 The main difference between DFTpower and DFTpower2, is that the latter allows for non-equidistant
559 frequencies for which the power spectrum will be computed.
560
561 @param time: time points, not necessarily equidistant
562 @type time: ndarray
563 @param signal: signal corresponding to the given time points
564 @type signal: ndarray
565 @param freqs: frequencies for which the power spectrum will be computed. Unit: inverse of 'time'.
566 @type freqs: ndarray
567 @return: power spectrum. Unit: square of unit of 'signal'
568 @rtype: ndarray
569 """
570
571 powerSpectrum = np.zeros(len(freqs))
572
573 for i, freq in enumerate(freqs):
574 arg = 2.0 * np.pi * freq * time
575 powerSpectrum[i] = np.sum(signal * np.cos(arg))**2 + np.sum(signal * np.sin(arg))**2
576
577 powerSpectrum = powerSpectrum * 4.0 / len(time)**2
578 return(powerSpectrum)
579
583
584 """
585 Compute Discrete Fourier Transform for unevenly spaced data ( Scargle, 1989).
586
587 Doesn't work yet!
588
589 It is recommended to start f0 at 0.
590 It is recommended to stop fn at the nyquist frequency
591
592 This makes use of a FORTRAN algorithm written by Scargle (1989).
593
594 @param times: observations times
595 @type times: numpy array
596 @param signal: observations
597 @type signal: numpy array
598 @param f0: start frequency
599 @type f0: float
600 @param fn: end frequency
601 @type fn: float
602 @param df: frequency step
603 @type df: float
604 @return: frequencies, dft, Re(dft), Im(dft)
605 """
606 f0 = 0.
607
608
609
610 nfreq = int((fn-f0)/(df))
611 print "Nfreq=",nfreq
612 tzero = times[0]
613 si = 1.
614 lfreq = 2*nfreq+1
615 mm = 2*nfreq
616 ftrx = np.zeros(mm)
617 ftix = np.zeros(mm)
618 om = np.zeros(mm)
619 w = np.zeros(mm)
620 wz = df
621 nn = len(times)
622
623
624 ftrx,ftix,om,w = pydft.ft(signal,times,wz,nfreq,si,lfreq,tzero,df,ftrx,ftix,om,w,nn,mm)
625
626 if f0==0:
627 ftrx[1:] *= np.sqrt(2)
628 ftix[1:] *= np.sqrt(2)
629 w[1:] *= np.sqrt(2)
630
631 cut_off = len(w)
632 for i in range(0,len(w))[::-1]:
633 if w[i] != 0: cut_off = i+1;break
634
635 om = om[:cut_off]
636 ftrx = ftrx[:cut_off]
637 ftix = ftix[:cut_off]
638 w = w[:cut_off]
639
640
641 T = times[-1]-times[0]
642 N = len(times)
643
644 w *= T/(2.*N)
645 ftrx *= T/(2.*N)
646 ftix *= T/(2.*N)
647
648 return om*2*np.pi,w,ftrx,ftix
649
652
653 """
654 Computes power spectrum of an equidistant time series 'signal'
655 using the FFT algorithm. The length of the time series need not
656 be a power of 2 (zero padding is done automatically).
657 Normalisation is such that a signal A*sin(2*pi*nu_0*t)
658 gives power A^2 at nu=nu_0 (IF nu_0 is in the 'freq' array)
659
660 @param signal: the time series [0..Ntime-1]
661 @type signal: ndarray
662 @param timestep: time step fo the equidistant time series
663 @type timestep: float
664 @return: frequencies and the power spectrum
665 @rtype: array,array
666
667 """
668
669
670
671
672 fourier = np.fft.rfft(signal)
673 Ntime = len(signal)
674 Nfreq = len(fourier)
675
676
677
678 power = np.abs(fourier)**2 * 4.0 / Ntime**2
679
680
681
682
683
684
685
686 freq = np.arange(float(Nfreq)) / (Nfreq-1) * 0.5 / timestep
687
688
689
690 return (freq, power)
691
700
701 """
702 Computes the power density of an equidistant time series 'signal',
703 using the FFT algorithm. The length of the time series need not
704 be a power of 2 (zero padding is done automatically).
705
706 @param signal: the time series [0..Ntime-1]
707 @type signal: ndarray
708 @param timestep: time step fo the equidistant time series
709 @type timestep: float
710 @return: frequencies and the power density spectrum
711 @rtype: array,array
712
713 """
714
715
716
717
718 fourier = np.fft.rfft(signal)
719 Ntime = len(signal)
720 Nfreq = len(fourier)
721
722
723
724 powerdensity = np.abs(fourier)**2 / Ntime * timestep
725
726
727
728
729
730
731
732 freq = np.arange(float(Nfreq)) / (Nfreq-1) * 0.5 / timestep
733
734
735
736 return (freq, powerdensity)
737
738
739
740
741
742
743
744
745
746 -def weightedpower(time, signal, weight, freq):
747
748 """
749 Compute the weighted power spectrum of a time signal.
750 For each given frequency a weighted sine fit is done using
751 chi-square minimization.
752
753 @param time: time points [0..Ntime-1]
754 @type time: ndarray
755 @param signal: observations [0..Ntime-1]
756 @type signal: ndarray
757 @param weight: 1/sigma_i^2 of observation
758 @type weight: ndarray
759 @param freq: frequencies [0..Nfreq-1] for which the power
760 needs to be computed
761 @type freq: ndarray
762 @return: weighted power [0..Nfreq-1]
763 @rtype: array
764
765 """
766
767 result = np.zeros(len(freq))
768
769 for i in range(len(freq)):
770 if (freq[i] != 0.0):
771 sine = np.sin(2.0*pi*freq[i]*time)
772 cosine = np.cos(2.0*pi*freq[i]*time)
773 a11= np.sum(weight*sine*sine)
774 a12 = np.sum(weight*sine*cosine)
775 a21 = a12
776 a22 = np.sum(weight*cosine*cosine)
777 b1 = np.sum(weight*signal*sine)
778 b2 = np.sum(weight*signal*cosine)
779 denominator = a11*a22-a12*a21
780 A = (b1*a22-b2*a12)/denominator
781 B = (b2*a11-b1*a21)/denominator
782 result[i] = A*A+B*B
783 else:
784 result[i] = np.sum(signal)/len(signal)
785
786 return(result)
787
788
789
790
791
792
793
794
795
796
797 @defaults_pergram
798 @parallel_pergram
799 @make_parallel
800 -def pdm(times, signal,f0=None,fn=None,df=None,Nbin=5,Ncover=2,
801 D=0,forbit=None,asini=None,e=None,omega=None,nmax=10):
802 """
803 Phase Dispersion Minimization of Jurkevich-Stellingwerf (1978)
804
805 This definition makes use of a Fortran routine written by Jan Cuypers and
806 Conny Aerts.
807
808 Inclusion of linear frequency shift by Pieter Degroote (see Cuypers 1986)
809
810 Inclusion of binary orbital motion by Pieter Degroote (see Shibahashi &
811 Kurtz 2012). When orbits are added, times must be in days, then asini is
812 in AU.
813
814 For circular orbits, give only forbit and asini.
815
816 @param times: time points
817 @type times: numpy array
818 @param signal: observations
819 @type signal: numpy array
820 @param f0: start frequency
821 @type f0: float
822 @param fn: stop frequency
823 @type fn: float
824 @param df: step frequency
825 @type df: float
826 @param Nbin: number of bins (default: 5)
827 @type Nbin: int
828 @param Ncover: number of covers (default: 1)
829 @type Ncover: int
830 @param D: linear frequency shift parameter
831 @type D: float
832 @return: frequencies, theta statistic
833 @rtype: array,array
834 """
835 T = times.ptp()
836 n = len(times)
837
838
839 xvar = signal.std()**2.
840 xx = (n-1) * xvar
841 nf = int((fn-f0) / df + 0.001) + 1
842 f1 = np.zeros(nf,'d')
843 s1 = np.zeros(nf,'d')
844
845
846
847 if D is None and asini is None:
848 f1, s1 = pyscargle.justel(signal,times,f0,df,Nbin,Ncover,xvar,xx,f1,s1,n,nf)
849
850 elif asini is None:
851 f1, s1 = pyscargle.justel2(signal,times,f0,df,Nbin,Ncover,xvar,xx,D,f1,s1,n,nf)
852
853 elif asini is not None and (e is None or e==0):
854 f1, s1 = pyscargle.justel3(signal,times,f0,df,Nbin,Ncover,xvar,xx,asini,
855 forbit,f1,s1,n,nf)
856
857 elif e>0:
858 forbit = 2*pi*forbit
859 ans,bns = np.array([[__ane__(n,e),__bne__(n,e)] for n in range(1,nmax+1)]).T
860 ksins = np.sqrt(ans**2*np.cos(omega)**2+bns**2*np.sin(omega)**2)
861 thns = np.arctan(bns/ans*np.tan(omega))
862 tau = -np.sum(bns*np.sin(omega))
863 f1, s1 = pyscargle.justel4(signal,times,f0,df,Nbin,Ncover,xvar,xx,asini,
864 forbit,e,omega,ksins,thns,tau,f1,s1,n,nf,nmax)
865
866
867
868 if not s1[0]: s1[0] = 1.
869
870 return f1, s1
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885 @defaults_pergram
886 @parallel_pergram
887 @make_parallel
888 -def box(times, signal, f0=None, fn=None, df=None, Nbin=10, qmi=0.005, qma=0.75 ):
889 """
890 Box-Least-Squares spectrum of Kovacs et al (2002).
891
892 [ see Kovacs, Zucker & Mazeh 2002, A&A, Vol. 391, 369 ]
893
894 This is the slightly modified version of the original BLS routine
895 by considering Edge Effect (EE) as suggested by
896 Peter R. McCullough [ pmcc@stsci.edu ].
897
898 This modification was motivated by considering the cases when
899 the low state (the transit event) happened to be devided between
900 the first and last bins. In these rare cases the original BLS
901 yields lower detection efficiency because of the lower number of
902 data points in the bin(s) covering the low state.
903
904 For further comments/tests see www.konkoly.hu/staff/kovacs.html
905
906 Transit fraction and precision are given by nb,qmi and qma
907
908 Remark: output parameter parameter contains:
909 [frequency,depth,transit fraction width,fractional start, fraction end]
910
911 @param times: observation times
912 @type times: numpy 1D array
913 @param signal: observations
914 @type signal: numpy 1D array
915 @param f0: start frequency
916 @type f0: float
917 @param fn: end frequency
918 @type fn: float
919 @param df: frequency step
920 @type df: float
921 @param Nbin: number of bins in the folded time series at any test period
922 @type Nbin: integer
923 @param qmi: minimum fractional transit length to be tested
924 @type qmi: 0<float<qma<1
925 @param qma: maximum fractional transit length to be tested
926 @type qma: 0<qmi<float<1
927 @return: frequencies, amplitude spectrum
928 @rtype: array,array
929 """
930
931 n = len(times)
932 T = times.ptp()
933 u = np.zeros(n)
934 v = np.zeros(n)
935
936
937 nf = (fn-f0)/df
938 if f0<2./T: f0=2./T
939
940
941 power,depth,qtran,in1,in2 = eebls.eebls(times,signal,u,v,nf,f0,df,Nbin,qmi,qma,n)
942 frequencies = np.linspace(f0,fn,nf)
943
944
945
946 return frequencies,power
947
948
949
950
951
952 @defaults_pergram
953 @parallel_pergram
954 @make_parallel
955 -def kepler(times,signal, f0=None, fn=None, df=None, e0=0., en=0.91, de=0.1,
956 errors=None, wexp=2, x00=0.,x0n=359.9):
957 """
958 Keplerian periodogram of Zucker et al (2010).
959
960 @param times: observation times
961 @type times: numpy 1D array
962 @param signal: observations
963 @type signal: numpy 1D array
964 @param f0: start frequency
965 @type f0: float
966 @param fn: end frequency
967 @type fn: float
968 @param df: frequency step
969 @type df: float
970 @param e0: start eccentricity
971 @type e0: float
972 @param en: end eccentricity
973 @type en: float
974 @param de: eccentricity step
975 @type de: float
976 @param x00: start x0
977 @type x00: float
978 @param x0n: end x0
979 @type x0n: float
980 @return: frequencies, amplitude spectrum
981 @rtype: array,array
982 """
983 T = times.ptp()
984 n = len(times)
985 if errors is None:
986 errors = np.ones(n)
987 maxstep = int((fn-f0)/df+1)
988
989
990 f1 = np.zeros(maxstep)
991 s1 = np.zeros(maxstep)
992 p1 = np.zeros(maxstep)
993 l1 = np.zeros(maxstep)
994 s2 = np.zeros(maxstep)
995 k2 = np.zeros(6)
996
997
998 pyKEP.kepler(times+0,signal+0,errors,f0,fn,df,wexp,e0,en,de,\
999 x00,x0n,f1,s1,p1,l1,s2,k2)
1000 return f1,s2
1001
1002
1003
1004
1005
1006 -def Zwavelet(time, signal, freq, position, sigma=10.0):
1007
1008 """
1009 Weighted Wavelet Z-transform of Foster (1996)
1010
1011 Computes "Weighted Wavelet Z-Transform"
1012 which is a type of time-frequency diagram more suitable
1013 for non-equidistant time series.
1014 See: G. Foster, 1996, Astronomical Journal, 112, 1709
1015
1016 The result can be plotted as a colorimage (imshow in pylab).
1017
1018 Mind the max, min order for position. It's often useful to try log(Z)
1019 and/or different sigmas.
1020
1021 @param time: time points [0..Ntime-1]
1022 @type time: ndarray
1023 @param signal: observed data points [0..Ntime-1]
1024 @type signal: ndarray
1025 @param freq: frequencies (omega/2pi in Foster's paper) [0..Nfreq-1]
1026 the array should not contain 0.0
1027 @type freq: ndarray
1028 @param position: time parameter: tau in Foster's paper [0..Npos-1]
1029 @type position: ndarray
1030 @param sigma: smoothing parameter in time domain: sigma in Foster's paper
1031 @type sigma: float
1032 @return: Z[0..Npos-1, 0..Nfreq-1]: the Z-transform: time-freq diagram
1033 @rtype: array
1034
1035 """
1036
1037 y = np.zeros(3)
1038 S = np.zeros([3,3])
1039 Z = np.zeros([len(position),len(freq)])
1040
1041 for i in range(len(position)):
1042
1043 tau = position[i]
1044 arg = 2.0*pi*(time-tau)
1045
1046 for j in range(len(freq)):
1047
1048 nu = freq[j]
1049
1050
1051
1052 weight = np.exp(-(time-tau)**2 * (nu / 2./sigma)**2)
1053 W = np.sum(weight)
1054
1055
1056
1057 cosine = np.cos(arg*nu)
1058 sine = np.sin(arg*nu)
1059
1060 print arg, nu
1061 print sine, cosine
1062
1063
1064
1065
1066 S[0,0] = 1.0
1067 S[0,1] = S[1,0] = np.sum(weight * 1.0 * cosine) / W
1068 S[0,2] = S[2,0] = np.sum(weight * 1.0 * sine) / W
1069 S[1,1] = np.sum(weight * cosine * cosine) / W
1070 S[1,2] = S[2,1] = np.sum(weight * cosine * sine) / W
1071 S[2,2] = np.sum(weight * sine * sine) / W
1072 print S
1073 invS = np.linalg.inv(S)
1074
1075
1076
1077 for k in range(3):
1078 y[k] = np.sum(weight * 1.0 * signal) / W * invS[k,0] \
1079 + np.sum(weight * cosine * signal) / W * invS[k,1] \
1080 + np.sum(weight * sine * signal) / W * invS[k,2]
1081
1082
1083
1084 model = y[0] + y[1] * cosine + y[2] * sine
1085
1086
1087
1088 Vsignal = np.sum(weight * signal**2) / W - (np.sum(weight * signal) / W)**2
1089 Vmodel = np.sum(weight * model**2) / W - (np.sum(weight * model) / W)**2
1090
1091
1092
1093 Neff = W**2 / np.sum(weight**2)
1094 Z[i,j] = (Neff - 3) * Vmodel / 2. / (Vsignal - Vmodel)
1095
1096
1097
1098 return Z
1099
1100
1101
1102
1103
1104 @defaults_pergram
1105 @parallel_pergram
1106 @make_parallel
1107 -def pdm_py(time, signal, f0=None, fn=None, df=None, Nbin=10, Ncover=5, D=0.):
1108
1109 """
1110 Computes the theta-statistics to do a Phase Dispersion Minimisation.
1111 See Stellingwerf R.F., 1978, ApJ, 224, 953)
1112
1113 Joris De Ridder
1114
1115 Inclusion of linear frequency shift by Pieter Degroote (see Cuypers 1986)
1116
1117 @param time: time points [0..Ntime-1]
1118 @type time: ndarray
1119 @param signal: observed data points [0..Ntime-1]
1120 @type signal: ndarray
1121 @param f0: start frequency
1122 @type f0: float
1123 @param fn: stop frequency
1124 @type fn: float
1125 @param df: step frequency
1126 @type df: float
1127 @param Nbin: the number of phase bins (with length 1/Nbin)
1128 @type Nbin: integer
1129 @param Ncover: the number of covers (i.e. bin shifts)
1130 @type Ncover: integer
1131 @param D: linear frequency shift parameter
1132 @type D: float
1133 @return: theta-statistic for each given frequency [0..Nfreq-1]
1134 @rtype: array
1135 """
1136 freq = np.arange(f0,fn+df,df)
1137
1138 Ntime = len(time)
1139 Nfreq = len(freq)
1140
1141 binsize = 1.0 / Nbin
1142 covershift = 1.0 / (Nbin * Ncover)
1143
1144 theta = np.zeros(Nfreq)
1145
1146 for i in range(Nfreq):
1147
1148
1149 phase = np.fmod((time - time[0]) * freq[i] + D/2.*time**2, 1.0)
1150
1151
1152
1153 Nempty = 0
1154
1155
1156
1157 for k in range(Nbin):
1158 for n in range(Ncover):
1159
1160
1161
1162
1163
1164 left = np.fmod(k * binsize + n * covershift, 1.0)
1165 right = np.fmod((k+1) * binsize + n * covershift, 1.0)
1166
1167
1168
1169 if (left < right):
1170 bindata = np.compress((left <= phase) & (phase < right), signal)
1171 else:
1172 bindata = np.compress(~((right <= phase) & (phase < left)), signal)
1173
1174
1175
1176 if (len(bindata) != 0):
1177 theta[i] += (len(bindata) - 1) * bindata.var()
1178 else:
1179 Nempty += 1
1180
1181
1182
1183 theta[i] /= Ncover * Ntime - (Ncover * Nbin - Nempty)
1184
1185
1186
1187 theta /= signal.var()
1188
1189
1190
1191 return freq,theta
1192
1193
1194
1195
1196 -def fasper_py(x,y,ofac,hifac, MACC=4):
1197 """ function fasper
1198 Given abscissas x (which need not be equally spaced) and ordinates
1199 y, and given a desired oversampling factor ofac (a typical value
1200 being 4 or larger). this routine creates an array wk1 with a
1201 sequence of nout increasing frequencies (not angular frequencies)
1202 up to hifac times the "average" Nyquist frequency, and creates
1203 an array wk2 with the values of the Lomb normalized periodogram at
1204 those frequencies. The arrays x and y are not altered. This
1205 routine also returns jmax such that wk2(jmax) is the maximum
1206 element in wk2, and prob, an estimate of the significance of that
1207 maximum against the hypothesis of random noise. A small value of prob
1208 indicates that a significant periodic signal is present.
1209
1210 Reference:
1211 Press, W. H. & Rybicki, G. B. 1989
1212 ApJ vol. 338, p. 277-280.
1213 Fast algorithm for spectral analysis of unevenly sampled data
1214 (1989ApJ...338..277P)
1215
1216 Arguments:
1217 X : Abscissas array, (e.g. an array of times).
1218 Y : Ordinates array, (e.g. corresponding counts).
1219 Ofac : Oversampling factor.
1220 Hifac : Hifac * "average" Nyquist frequency = highest frequency
1221 for which values of the Lomb normalized periodogram will
1222 be calculated.
1223
1224 Returns:
1225 Wk1 : An array of Lomb periodogram frequencies.
1226 Wk2 : An array of corresponding values of the Lomb periodogram.
1227 Nout : Wk1 & Wk2 dimensions (number of calculated frequencies)
1228 Jmax : The array index corresponding to the MAX( Wk2 ).
1229 Prob : False Alarm Probability of the largest Periodogram value
1230 MACC : Number of interpolation points per 1/4 cycle
1231 of highest frequency
1232
1233 History:
1234 02/23/2009, v1.0, MF
1235 Translation of IDL code (orig. Numerical recipies)
1236 """
1237
1238 n = long(len(x))
1239 if n != len(y):
1240 print 'Incompatible arrays.'
1241 return
1242
1243 nout = 0.5*ofac*hifac*n
1244 nfreqt = long(ofac*hifac*n*MACC)
1245 nfreq = 64L
1246
1247 while nfreq < nfreqt:
1248 nfreq = 2*nfreq
1249
1250 ndim = long(2*nfreq)
1251
1252
1253 ave = y.mean()
1254
1255 var = ((y-y.mean())**2).sum()/(len(y)-1)
1256
1257 xmin = x.min()
1258 xmax = x.max()
1259 xdif = xmax-xmin
1260
1261
1262 wk1 = np.zeros(ndim, dtype='complex')
1263 wk2 = np.zeros(ndim, dtype='complex')
1264
1265 fac = ndim/(xdif*ofac)
1266 fndim = ndim
1267 ck = ((x-xmin)*fac) % fndim
1268 ckk = (2.0*ck) % fndim
1269
1270 for j in range(0L, n):
1271 __spread__(y[j]-ave,wk1,ndim,ck[j],MACC)
1272 __spread__(1.0,wk2,ndim,ckk[j],MACC)
1273
1274
1275 wk1 = np.fft.ifft( wk1 )*len(wk1)
1276 wk2 = np.fft.ifft( wk2 )*len(wk1)
1277
1278 wk1 = wk1[1:nout+1]
1279 wk2 = wk2[1:nout+1]
1280 rwk1 = wk1.real
1281 iwk1 = wk1.imag
1282 rwk2 = wk2.real
1283 iwk2 = wk2.imag
1284
1285 df = 1.0/(xdif*ofac)
1286
1287
1288 hypo2 = 2.0 * abs( wk2 )
1289 hc2wt = rwk2/hypo2
1290 hs2wt = iwk2/hypo2
1291
1292 cwt = np.sqrt(0.5+hc2wt)
1293 swt = np.sign(hs2wt)*(np.sqrt(0.5-hc2wt))
1294 den = 0.5*n+hc2wt*rwk2+hs2wt*iwk2
1295 cterm = (cwt*rwk1+swt*iwk1)**2./den
1296 sterm = (cwt*iwk1-swt*rwk1)**2./(n-den)
1297
1298 wk1 = df*(np.arange(nout, dtype='float')+1.)
1299 wk2 = (cterm+sterm)/(2.0*var)
1300 pmax = wk2.max()
1301 jmax = wk2.argmax()
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312 expy = np.exp(-pmax)
1313 effm = 2.0*(nout)/ofac
1314 prob = effm*expy
1315
1316 if prob > 0.01:
1317 prob = 1.0-(1.0-expy)**effm
1318
1319 return wk1,wk2,nout,jmax,prob
1320
1325
1326 """
1327 Computes the modulus square of the window function of a set of
1328 time points at the given frequencies. The time point need not be
1329 equidistant. The normalisation is such that 1.0 is returned at
1330 frequency 0.
1331
1332 @param time: time points [0..Ntime-1]
1333 @type time: ndarray
1334 @param freq: frequency points. Units: inverse unit of 'time' [0..Nfreq-1]
1335 @type freq: ndarray
1336 @return: |W(freq)|^2 [0..Nfreq-1]
1337 @rtype: array
1338
1339 """
1340
1341 Ntime = len(time)
1342 Nfreq = len(freq)
1343 winkernel = np.empty_like(freq)
1344
1345 for i in range(Nfreq):
1346 winkernel[i] = np.sum(np.cos(2.0*pi*freq[i]*time))**2 \
1347 + np.sum(np.sin(2.0*pi*freq[i]*time))**2
1348
1349
1350
1351 return winkernel/Ntime**2
1352
1413
1416 """
1417 Given an array yy(0:n-1), extirpolate (spread) a value y into
1418 m actual array elements that best approximate the "fictional"
1419 (i.e., possible noninteger) array element number x. The weights
1420 used are coefficients of the Lagrange interpolating polynomial
1421 Arguments:
1422 y :
1423 yy :
1424 n :
1425 x :
1426 m :
1427 Returns:
1428
1429 """
1430 nfac=[0,1,1,2,6,24,120,720,5040,40320,362880]
1431 if m > 10. :
1432 print 'factorial table too small in spread'
1433 return
1434
1435 ix=long(x)
1436 if x == float(ix):
1437 yy[ix]=yy[ix]+y
1438 else:
1439 ilo = long(x-0.5*float(m)+1.0)
1440 ilo = min( max( ilo , 1 ), n-m+1 )
1441 ihi = ilo+m-1
1442 nden = nfac[m]
1443 fac=x-ilo
1444 for j in range(ilo+1,ihi+1): fac = fac*(x-j)
1445 yy[ihi] = yy[ihi] + y*fac/(nden*(x-ihi))
1446 for j in range(ihi-1,ilo-1,-1):
1447 nden=(nden/(j+1-ilo))*(j-ihi)
1448 yy[j] = yy[j] + y*fac/(nden*(x-j))
1449
1451 return 2.*np.sqrt(1-e**2)/e/n*jn(n,n*e)
1452
1454 return 1./n*(jn(n-1,n*e)-jn(n+1,n*e))
1455
1458 """
1459 Returns the peak false alarm probabilities
1460
1461 Hence the lower is the probability and the more significant is the peak
1462 """
1463 expy = np.exp(-wk2)
1464 effm = 2.0*(nout)/ofac
1465 sig = effm*expy
1466 ind = (np.sig > 0.01).nonzero()
1467 sig[ind] = 1.0-(1.0-expy[ind])**effm
1468 return sig
1469
1472 """
1473 Compute the probability to observe a peak in the Scargle periodogram.
1474
1475 If C{correct_for_frange=True}, the Bonferroni correction will be applied
1476 to a smaller number of frequencies (i.e. the independent number of
1477 frequencies in C{freqs}). To be conservative, set C{correct_for_frange=False}.
1478
1479 Example simulation:
1480
1481 >>> times = np.linspace(0,1,5000)
1482 >>> N = 500
1483 >>> probs = np.zeros(N)
1484 >>> peaks = np.zeros(N)
1485 >>> for i in range(N):
1486 ... signal = np.random.normal(size=len(times))
1487 ... f,s = scargle(times,signal,threads='max',norm='distribution')
1488 ... peaks[i] = s.max()
1489 ... probs[i] = scargle_probability(s.max(),times,f)
1490
1491 Now make a plot:
1492
1493 >>> p = pl.figure()
1494 >>> p = pl.subplot(131)
1495 >>> p = pl.plot(probs,'ko')
1496 >>> p = pl.plot([0,N],[0.01,0.01],'r-',lw=2)
1497 >>> p = pl.subplot(132)
1498 >>> p = pl.plot(peaks[np.argsort(peaks)],probs[np.argsort(peaks)],'ro')
1499 >>> p = pl.plot(peaks[np.argsort(peaks)],1-(1-np.exp(-np.sort(peaks)))**(10000.),'g-')
1500 >>> #p = pl.plot(peaks[np.argsort(peaks)],1-(1-(1-np.sort(peaks)/2500.)**2500.)**(10000.),'b--')
1501 >>> p = pl.subplot(133)
1502 >>> for i in np.logspace(-3,0,100):
1503 ... p = pl.plot([i*100],[np.sum(probs<i)/float(N)*100],'ko')
1504 >>> p = pl.plot([1e-6,100],[1e-6,100],'r-',lw=2)
1505 >>> p = pl.xlabel('Should observe this many points below threshold')
1506 >>> p = pl.ylabel('Observed this many points below threshold')
1507
1508 ]]include figure]]ivs_timeseries_pergrams_prob.png]
1509
1510 """
1511
1512 nr_obs = len(times)
1513 ni = 2*nr_obs
1514
1515
1516 if correct_for_frange:
1517 nyqstat = kwargs.pop('nyqstat',np.min)
1518 nyquist = getNyquist(times,nyqstat=nyqstat)
1519 ni = int(freqs.ptp()/nyquist*ni)
1520
1521 p_value = 1. - (1.- np.exp(-peak_value))**ni
1522 return p_value
1523
1524
1525 if __name__=="__main__":
1526 import sys
1527 import pylab as pl
1528 from ivs.aux import loggers
1529 logger = loggers.get_basic_logger()
1530
1531
1532 if '--test' in sys.argv[1] or '-t' in sys.argv[1]:
1533 import doctest
1534 doctest.testmod()
1535 pl.show()
1536
1537 else:
1538
1539 method,args,kwargs = argkwargparser.parse()
1540 print "Running method %s with arguments %s and keyword arguments %s"%(method,args,kwargs)
1541 if '--help' in args or 'help' in args or 'help' in kwargs:
1542 sys.exit()
1543 times,signal = ascii.read2array(kwargs.pop('infile')).T[:2]
1544 freqs,ampls = globals()[method](times,signal, **kwargs)
1545 pl.plot(freqs,ampls,'k-')
1546 pl.show()
1547