Home       Trees       Indices       Help   
Package ivs :: Package timeseries :: Module pergrams
[hide private]
[frames] | no frames]

Module pergrams

source code

Contains many different periodogram calculations

Section 1. Basic usage

Given a time series of the form

>>> times = np.linspace(0,1,1000)
>>> signal = np.sin(times)

The basic interface is

>>> freq,ampl = scargle(times,signal)

If you give no extra information, default values for the start, end and step frequency will be chosen. See the 'defaults_pergram' decorator for a list of keywords common to all periodogram calculations.

Many periodograms can be computed in parallel, by supplying an extra keyword 'threads':

>>> freq,ampl = scargle(times,signal,threads=2)

Warning: the timeseries must be sorted in time and cannot contain the same timepoint twice. Otherwise, a 'ValueError, concatenation problem' can occur.

If something goes wrong in the periodogram computation, be sure to run check_input on your input data. This will print out some basic diagnostics to see if your data are valid.

Section 2. Nyquist frequency

The periodogram functions are written such that they never exceed the value of the Nyquist frequency. This behaviour can be changed.

By default, the Nyquist frequency is defined as half of the inverse of the smallest time step in the data. That is, the nyq_stat variable is set to the function np.min. If you prefer the Nyquist to be defined as the median, set the nyq_stat variable for any periodogram to np.median. If you want a more complex or self-defined function, that is also acceptable. If you give a number to nyq_stat, nothing will be computed but that value will be considered the nyquist frequency.

Section 3. Speed comparison

>>> import time

The timeseries is generated for N=500,1000,2000,5000,10000,20000,50000

>>> techniques = ['scargle','deeming','gls','fasper',
...               'schwarzenberg_czerny','pdm','box']
>>> Ns = [1000,2000,5000,10000,20000,30000]
>>> for tech in techniques[:0]:
...     freqs = []
...     clock = []
...     for N in Ns:
...         times = np.linspace(0,1,N)
...         signal = np.sin(times) + np.random.normal(size=N)
...         c0 = time.time()
...         freq,ampl = locals()[tech](times,signal)
...         clock.append(time.time()-c0)
...         freqs.append(len(freq))
...     p = pl.plot(freqs,clock,'o-',label=tech)
...     print tech,np.polyfit(np.log10(np.array(freqs)),np.log10(np.array(clock)),1)
>>> #p = pl.legend(loc='best',fancybox=True)
>>> #q = p.get_frame().set_alpha(0.5)
>>> #p,q = pl.xlabel('Number of frequencies'),pl.ylabel('Seconds')

]]include figure]]ivs_timeseries_pergram_speeds.png]

Section 2. Periodogram comparison

We generate a sinusoidal signal, add some noise and compute the periodogram with the different methods.

>>> times = np.sort(np.random.uniform(size=500,low=0,high=100))
>>> signal = np.sin(2*np.pi/10.*times)
>>> signal += np.random.normal(size=500)

Fourier based periodgrams suited for unequidistant data: scargle, deeming, fasper and schwarzenberg_czerny.

>>> f1,a1 = scargle(times,signal,fn=0.35)
>>> f2,a2 = fasper(times,signal,fn=0.35)
>>> f3,a3 = deeming(times,signal,fn=0.35)
>>> f4,a4 = scargle(times,signal,norm='power',fn=0.35)
>>> f5,a5 = schwarzenberg_czerny(times,signal,fn=0.35,nh=1,mode=2)
>>> p = pl.figure()
>>> p = pl.subplot(121)
>>> p = pl.plot(f1,a1,'k-',lw=4,label='scargle')
>>> p = pl.plot(f2,a2,'r--',lw=4,label='fasper')
>>> p = pl.plot(f3,a3,'b-',lw=2,label='deeming')
>>> p = pl.xlim(f1[0],f1[-1])
>>> leg = pl.legend(fancybox=True,loc='best')
>>> leg.get_frame().set_alpha(0.5)
>>> p = pl.subplot(122)
>>> p = pl.plot(f4,a4,'k-',lw=2,label='scargle')
>>> p = pl.plot(f5,a5/2.,'r--',lw=2,label='schwarzenberg-czerny')
>>> p = pl.xlim(f1[0],f1[-1])
>>> leg = pl.legend(fancybox=True,loc='best')
>>> leg.get_frame().set_alpha(0.5)

]]include figure]]ivs_timeseries_pergrams_fourier.png]

Least-square fitting: gls and kepler.

>>> f1,a1 = gls(times,signal,fn=0.35)
>>> f2,a2 = kepler(times,signal,fn=0.35)
>>> p = pl.figure()
>>> p = pl.plot(f1,a1,'k-',lw=2,label='gls')
>>> p = pl.plot(f2,a2,'r--',lw=2,label='kepler')
>>> p = pl.xlim(f1[0],f1[-1])
>>> leg = pl.legend(fancybox=True,loc='best')
>>> leg.get_frame().set_alpha(0.5)

]]include figure]]ivs_timeseries_pergrams_lsq.png]

Phase folding techniques: box and pdm

>>> f1,a1 = box(times,signal,fn=0.35)
>>> f2,a2 = pdm(times,signal,fn=0.35)
>>> p = pl.figure()
>>> p = pl.plot(f1,a1,'k-',lw=2,label='box')
>>> p = pl.plot(f2,a2,'r--',lw=2,label='pdm')
>>> p = pl.xlim(f1[0],f1[-1])
>>> leg = pl.legend(fancybox=True,loc='best')
>>> leg.get_frame().set_alpha(0.5)

]]include figure]]ivs_timeseries_pergrams_phase.png]

Functions [hide private]
    Periodograms
array,array
scargle(times, signal, f0=None, fn=None, df=None, norm='amplitude', weights=None, single=False)
Scargle periodogram of Scargle (1982).		 source code
array,array
fasper(times, signal, f0=None, fn=None, df=None, single=True, norm='amplitude')
Fasper periodogram from Numerical Recipes.		 source code
array,array
deeming(times, signal, f0=None, fn=None, df=None, norm='amplitude')
Deeming periodogram of Deeming et al (1975).		 source code
array,array
gls(times, signal, f0=None, fn=None, df=None, errors=None, wexp=2)
Generalised Least Squares periodogram of Zucher et al (2010).		 source code
array,array
clean(times, signal, f0=None, fn=None, df=None, freqbins=None, niter=10., gain=1.0)
Cleaned Fourier periodogram of Roberts et al (1987)		 source code
array,array
schwarzenberg_czerny(times, signal, f0=None, fn=None, df=None, nh=2, mode=1)
Multi harmonic periodogram of Schwarzenberg-Czerny (1996).		 source code
array
DFTpower(time, signal, f0=None, fn=None, df=None, full_output=False)
Computes the modulus square of the fourier transform.		 source code
ndarray
DFTpower2(time, signal, freqs)
Computes the power spectrum of a signal using a discrete Fourier transform.		 source code
 
DFTscargle(times, signal, f0, fn, df)
Compute Discrete Fourier Transform for unevenly spaced data ( Scargle, 1989).		 source code
array,array
FFTpower(signal, timestep)
Computes power spectrum of an equidistant time series 'signal' using the FFT algorithm.		 source code
array,array
FFTpowerdensity(signal, timestep)
Computes the power density of an equidistant time series 'signal', using the FFT algorithm.		 source code
array
weightedpower(time, signal, weight, freq)
Compute the weighted power spectrum of a time signal.		 source code
array,array
pdm(times, signal, f0=None, fn=None, df=None, Nbin=5, Ncover=2, D=0, forbit=None, asini=None, e=None, omega=None, nmax=10)
Phase Dispersion Minimization of Jurkevich-Stellingwerf (1978)		 source code
array,array
box(times, signal, f0=None, fn=None, df=None, Nbin=10, qmi=0.005, qma=0.75)
Box-Least-Squares spectrum of Kovacs et al (2002).		 source code
array,array
kepler(times, signal, f0=None, fn=None, df=None, e0=0., en=0.91, de=0.1, errors=None, wexp=2, x00=0., x0n=359.9)
Keplerian periodogram of Zucker et al (2010).		 source code
array
Zwavelet(time, signal, freq, position, sigma=10.0)
Weighted Wavelet Z-transform of Foster (1996)		 source code
    Pure Python versions
array
pdm_py(time, signal, f0=None, fn=None, df=None, Nbin=10, Ncover=5, D=0.)
Computes the theta-statistics to do a Phase Dispersion Minimisation.		 source code
 
fasper_py(x, y, ofac, hifac, MACC=4)
function fasper Given abscissas x (which need not be equally spaced) and ordinates y, and given a desired oversampling factor ofac (a typical value being 4 or larger).		 source code
    Helper functions
array
windowfunction(time, freq)
Computes the modulus square of the window function of a set of time points at the given frequencies.		 source code
 
check_input(times, signal, **kwargs)
Check the input arguments for periodogram calculations for mistakes.		 source code
 
__spread__(y, yy, n, x, m)
Given an array yy(0:n-1), extirpolate (spread) a value y into m actual array elements that best approximate the "fictional" (i.e., possible noninteger) array element number x.		 source code
 
__ane__(n, e) source code
 
__bne__(n, e) source code
 
getSignificance(wk1, wk2, nout, ofac)
Returns the peak false alarm probabilities		 source code
 
scargle_probability(peak_value, times, freqs, correct_for_frange=False, **kwargs)
Compute the probability to observe a peak in the Scargle periodogram.		 source code
Variables [hide private]
  logger = logging.getLogger("TS.PERGRAMS")
Function Details [hide private]

scargle(times, signal, f0=None, fn=None, df=None, norm='amplitude', weights=None, single=False)

source code 

Scargle periodogram of Scargle (1982).

Several options are available (possibly combined):

  1. weighted Scargle
  2. Amplitude spectrum
  3. Distribution power spectrum
  4. Traditional power spectrum
  5. Power density spectrum (see Kjeldsen, 2005 or Carrier, 2010)

This definition makes use of a Fortran-routine written by Jan Cuypers, Conny Aerts and Peter De Cat. A slightly adapted version is used for the weighted version (adapted by Pieter Degroote).

Through the option "norm", it's possible to norm the periodogram as to get a periodogram that has a known statistical distribution. Usually, this norm is the variance of the data (NOT of the noise or residuals, see Schwarzenberg- Czerny 1998!).

Also, it is possible to retrieve the power density spectrum in units of [ampl**2/frequency]. In this routine, the normalisation constant is taken to be the total time span T. Kjeldsen (2005) chooses to multiply the power by the 'effective length of the observing run', which is calculated as the reciprocal of the area under spectral window (in power, and take 2*Nyquist as upper frequency value).

REMARK: this routine does not automatically remove the average. It is the user's responsibility to do this adequately: e.g. subtract a weighted average if one computes the weighted periodogram!!

Parameters:
  • times (numpy array) - time points
  • signal (numpy array) - observations
  • weights (numpy array) - weights of the datapoints
  • norm (str) - type of normalisation
  • f0 (float) - start frequency
  • fn (float) - stop frequency
  • df (float) - step frequency
Returns: array,array
frequencies, amplitude spectrum
Decorators:
  • @defaults_pergram
  • @parallel_pergram
  • @make_parallel

fasper(times, signal, f0=None, fn=None, df=None, single=True, norm='amplitude')

source code 

Fasper periodogram from Numerical Recipes.

Normalisation here is not correct!!

Parameters:
  • times (numpy array) - time points
  • signal (numpy array) - observations
  • f0 (float) - start frequency
  • fn (float) - stop frequency
  • df (float) - step frequency
Returns: array,array
frequencies, amplitude spectrum
Decorators:
  • @defaults_pergram
  • @parallel_pergram
  • @make_parallel

deeming(times, signal, f0=None, fn=None, df=None, norm='amplitude')

source code 

Deeming periodogram of Deeming et al (1975).

Thanks to Jan Cuypers

Parameters:
  • times (numpy array) - time points
  • signal (numpy array) - observations
  • norm (str) - type of normalisation
  • f0 (float) - start frequency
  • fn (float) - stop frequency
  • df (float) - step frequency
Returns: array,array
frequencies, amplitude spectrum
Decorators:
  • @defaults_pergram
  • @parallel_pergram
  • @make_parallel

gls(times, signal, f0=None, fn=None, df=None, errors=None, wexp=2)

source code 

Generalised Least Squares periodogram of Zucher et al (2010).

Parameters:
  • times (numpy array) - time points
  • signal (numpy array) - observations
  • f0 (float) - start frequency
  • fn (float) - stop frequency
  • df (float) - step frequency
Returns: array,array
frequencies, amplitude spectrum
Decorators:
  • @defaults_pergram
  • @parallel_pergram
  • @make_parallel

clean(times, signal, f0=None, fn=None, df=None, freqbins=None, niter=10., gain=1.0)

source code 

Cleaned Fourier periodogram of Roberts et al (1987)

Parallization probably isn't such a good idea here because of the frequency bins.

Fortran module probably from John Telting.

Should always start from zero, so f0 is not an option

Generate some signal with heavy window side lobes:

>>> times_ = np.linspace(0,150,1000)
>>> times = np.array([times_[i] for i in xrange(len(times_)) if (i%10)>7])
>>> signal = np.sin(2*pi/10*times) + np.random.normal(size=len(times))

Compute the scargle periodogram as a reference, and compare with the the CLEAN periodogram with different gains.

>>> niter,freqbins = 10,[0,1.2]
>>> p1 = scargle(times,signal,fn=1.2,norm='amplitude',threads=2)
>>> p2 = clean(times,signal,fn=1.2,gain=1.0,niter=niter,freqbins=freqbins)
>>> p3 = clean(times,signal,fn=1.2,gain=0.1,niter=niter,freqbins=freqbins)

Make a figure of the result:

>>> p=pl.figure()
>>> p=pl.plot(p1[0],p1[1],'k-',label="Scargle")
>>> p=pl.plot(p2[0],p2[1],'r-',label="Clean (g=1.0)")
>>> p=pl.plot(p3[0],p3[1],'b-',label="Clean (g=0.1)")
>>> p=pl.legend()

]]include figure]]ivs_timeseries_pergrams_clean.png]

Parameters:
  • freqbins (list or array) - frequency bins for clean computation
  • niter (integer) - number of iterations
  • gain (float between 0 (no cleaning) and 1 (full cleaning)) - gain for clean computation
Returns: array,array
frequencies, amplitude spectrum
Decorators:
  • @defaults_pergram
  • @parallel_pergram
  • @make_parallel

schwarzenberg_czerny(times, signal, f0=None, fn=None, df=None, nh=2, mode=1)

source code 

Multi harmonic periodogram of Schwarzenberg-Czerny (1996).

This periodogram follows an F-distribution, so it is possible to perform hypothesis testing.

If the number of the number of harmonics is 1, then this peridogram reduces to the Lomb-Scargle periodogram except for its better statistic behaviour. This script uses a Fortran procedure written by Schwarzenberg-Czerny.

Modes:

  • mode=1: AoV Fisher-Snedecor F(df1,df2) statistic
  • mode=2: total power fitted in all harmonics for mode=2
Parameters:
  • times (numpy 1d array) - list of observations times
  • signal (numpy 1d array) - list of observations
  • f0 (float) - start frequency (cycles per day) (default: 0.)
  • fn (float) - stop frequency (cycles per day) (default: 10.)
  • df (float) - step frequency (cycles per day) (default: 0.001)
  • nh (integer) - number of harmonics to take into account
Returns: array,array
frequencies, f-statistic
Decorators:
  • @defaults_pergram
  • @parallel_pergram
  • @make_parallel

DFTpower(time, signal, f0=None, fn=None, df=None, full_output=False)

source code 

Computes the modulus square of the fourier transform.

Unit: square of the unit of signal. Time points need not be equidistant. The normalisation is such that a signal A*sin(2*pi*nu_0*t) gives power A^2 at nu=nu_0

Parameters:
  • time (ndarray) - time points [0..Ntime-1]
  • signal (ndarray) - signal [0..Ntime-1]
  • f0 (float) - the power is computed for the frequencies freq = arange(f0,fn,df)
  • fn (float) - see f0
  • df (float) - see f0
Returns: array
power spectrum of the signal

DFTpower2(time, signal, freqs)

source code 

Computes the power spectrum of a signal using a discrete Fourier transform.

The main difference between DFTpower and DFTpower2, is that the latter allows for non-equidistant frequencies for which the power spectrum will be computed.

Parameters:
  • time (ndarray) - time points, not necessarily equidistant
  • signal (ndarray) - signal corresponding to the given time points
  • freqs (ndarray) - frequencies for which the power spectrum will be computed. Unit: inverse of 'time'.
Returns: ndarray
power spectrum. Unit: square of unit of 'signal'

DFTscargle(times, signal, f0, fn, df)

source code 

Compute Discrete Fourier Transform for unevenly spaced data ( Scargle, 1989).

Doesn't work yet!

It is recommended to start f0 at 0. It is recommended to stop fn at the nyquist frequency

This makes use of a FORTRAN algorithm written by Scargle (1989).

Parameters:
  • times (numpy array) - observations times
  • signal (numpy array) - observations
  • f0 (float) - start frequency
  • fn (float) - end frequency
  • df (float) - frequency step
Returns:
frequencies, dft, Re(dft), Im(dft)

FFTpower(signal, timestep)

source code 

Computes power spectrum of an equidistant time series 'signal' using the FFT algorithm. The length of the time series need not be a power of 2 (zero padding is done automatically). Normalisation is such that a signal A*sin(2*pi*nu_0*t) gives power A^2 at nu=nu_0 (IF nu_0 is in the 'freq' array)

Parameters:
  • signal (ndarray) - the time series [0..Ntime-1]
  • timestep (float) - time step fo the equidistant time series
Returns: array,array
frequencies and the power spectrum

FFTpowerdensity(signal, timestep)

source code 

Computes the power density of an equidistant time series 'signal', using the FFT algorithm. The length of the time series need not be a power of 2 (zero padding is done automatically).

Parameters:
  • signal (ndarray) - the time series [0..Ntime-1]
  • timestep (float) - time step fo the equidistant time series
Returns: array,array
frequencies and the power density spectrum

weightedpower(time, signal, weight, freq)

source code 

Compute the weighted power spectrum of a time signal. For each given frequency a weighted sine fit is done using chi-square minimization.

Parameters:
  • time (ndarray) - time points [0..Ntime-1]
  • signal (ndarray) - observations [0..Ntime-1]
  • weight (ndarray) - 1/sigma_i^2 of observation
  • freq (ndarray) - frequencies [0..Nfreq-1] for which the power needs to be computed
Returns: array
weighted power [0..Nfreq-1]

pdm(times, signal, f0=None, fn=None, df=None, Nbin=5, Ncover=2, D=0, forbit=None, asini=None, e=None, omega=None, nmax=10)

source code 

Phase Dispersion Minimization of Jurkevich-Stellingwerf (1978)

This definition makes use of a Fortran routine written by Jan Cuypers and Conny Aerts.

Inclusion of linear frequency shift by Pieter Degroote (see Cuypers 1986)

Inclusion of binary orbital motion by Pieter Degroote (see Shibahashi & Kurtz 2012). When orbits are added, times must be in days, then asini is in AU.

For circular orbits, give only forbit and asini.

Parameters:
  • times (numpy array) - time points
  • signal (numpy array) - observations
  • f0 (float) - start frequency
  • fn (float) - stop frequency
  • df (float) - step frequency
  • Nbin (int) - number of bins (default: 5)
  • Ncover (int) - number of covers (default: 1)
  • D (float) - linear frequency shift parameter
Returns: array,array
frequencies, theta statistic
Decorators:
  • @defaults_pergram
  • @parallel_pergram
  • @make_parallel

box(times, signal, f0=None, fn=None, df=None, Nbin=10, qmi=0.005, qma=0.75)

source code 

Box-Least-Squares spectrum of Kovacs et al (2002).

[ see Kovacs, Zucker & Mazeh 2002, A&A, Vol. 391, 369 ]

This is the slightly modified version of the original BLS routine by considering Edge Effect (EE) as suggested by Peter R. McCullough [ pmcc@stsci.edu ].

This modification was motivated by considering the cases when the low state (the transit event) happened to be devided between the first and last bins. In these rare cases the original BLS yields lower detection efficiency because of the lower number of data points in the bin(s) covering the low state.

For further comments/tests see www.konkoly.hu/staff/kovacs.html

Transit fraction and precision are given by nb,qmi and qma

Remark: output parameter parameter contains: [frequency,depth,transit fraction width,fractional start, fraction end]

Parameters:
  • times (numpy 1D array) - observation times
  • signal (numpy 1D array) - observations
  • f0 (float) - start frequency
  • fn (float) - end frequency
  • df (float) - frequency step
  • Nbin (integer) - number of bins in the folded time series at any test period
  • qmi (0<float<qma<1) - minimum fractional transit length to be tested
  • qma (0<qmi<float<1) - maximum fractional transit length to be tested
Returns: array,array
frequencies, amplitude spectrum
Decorators:
  • @defaults_pergram
  • @parallel_pergram
  • @make_parallel

kepler(times, signal, f0=None, fn=None, df=None, e0=0., en=0.91, de=0.1, errors=None, wexp=2, x00=0., x0n=359.9)

source code 

Keplerian periodogram of Zucker et al (2010).

Parameters:
  • times (numpy 1D array) - observation times
  • signal (numpy 1D array) - observations
  • f0 (float) - start frequency
  • fn (float) - end frequency
  • df (float) - frequency step
  • e0 (float) - start eccentricity
  • en (float) - end eccentricity
  • de (float) - eccentricity step
  • x00 (float) - start x0
  • x0n (float) - end x0
Returns: array,array
frequencies, amplitude spectrum
Decorators:
  • @defaults_pergram
  • @parallel_pergram
  • @make_parallel

Zwavelet(time, signal, freq, position, sigma=10.0)

source code 

Weighted Wavelet Z-transform of Foster (1996)

Computes "Weighted Wavelet Z-Transform" which is a type of time-frequency diagram more suitable for non-equidistant time series. See: G. Foster, 1996, Astronomical Journal, 112, 1709

The result can be plotted as a colorimage (imshow in pylab).

Mind the max, min order for position. It's often useful to try log(Z) and/or different sigmas.

Parameters:
  • time (ndarray) - time points [0..Ntime-1]
  • signal (ndarray) - observed data points [0..Ntime-1]
  • freq (ndarray) - frequencies (omega/2pi in Foster's paper) [0..Nfreq-1] the array should not contain 0.0
  • position (ndarray) - time parameter: tau in Foster's paper [0..Npos-1]
  • sigma (float) - smoothing parameter in time domain: sigma in Foster's paper
Returns: array
Z[0..Npos-1, 0..Nfreq-1]: the Z-transform: time-freq diagram

pdm_py(time, signal, f0=None, fn=None, df=None, Nbin=10, Ncover=5, D=0.)

source code 

Computes the theta-statistics to do a Phase Dispersion Minimisation. See Stellingwerf R.F., 1978, ApJ, 224, 953)

Joris De Ridder

Inclusion of linear frequency shift by Pieter Degroote (see Cuypers 1986)

Parameters:
  • time (ndarray) - time points [0..Ntime-1]
  • signal (ndarray) - observed data points [0..Ntime-1]
  • f0 (float) - start frequency
  • fn (float) - stop frequency
  • df (float) - step frequency
  • Nbin (integer) - the number of phase bins (with length 1/Nbin)
  • Ncover (integer) - the number of covers (i.e. bin shifts)
  • D (float) - linear frequency shift parameter
Returns: array
theta-statistic for each given frequency [0..Nfreq-1]
Decorators:
  • @defaults_pergram
  • @parallel_pergram
  • @make_parallel

fasper_py(x, y, ofac, hifac, MACC=4)

source code  
function fasper
    Given abscissas x (which need not be equally spaced) and ordinates
    y, and given a desired oversampling factor ofac (a typical value
    being 4 or larger). this routine creates an array wk1 with a
    sequence of nout increasing frequencies (not angular frequencies)
    up to hifac times the "average" Nyquist frequency, and creates
    an array wk2 with the values of the Lomb normalized periodogram at
    those frequencies. The arrays x and y are not altered. This
    routine also returns jmax such that wk2(jmax) is the maximum
    element in wk2, and prob, an estimate of the significance of that
    maximum against the hypothesis of random noise. A small value of prob
    indicates that a significant periodic signal is present.

Reference: 
    Press, W. H. & Rybicki, G. B. 1989
    ApJ vol. 338, p. 277-280.
    Fast algorithm for spectral analysis of unevenly sampled data
    (1989ApJ...338..277P)

Arguments:
    X   : Abscissas array, (e.g. an array of times).
    Y   : Ordinates array, (e.g. corresponding counts).
    Ofac : Oversampling factor.
    Hifac : Hifac * "average" Nyquist frequency = highest frequency
        for which values of the Lomb normalized periodogram will
        be calculated.
    
Returns:
    Wk1 : An array of Lomb periodogram frequencies.
    Wk2 : An array of corresponding values of the Lomb periodogram.
    Nout : Wk1 & Wk2 dimensions (number of calculated frequencies)
    Jmax : The array index corresponding to the MAX( Wk2 ).
    Prob : False Alarm Probability of the largest Periodogram value
    MACC : Number of interpolation points per 1/4 cycle
            of highest frequency

History:
    02/23/2009, v1.0, MF
    Translation of IDL code (orig. Numerical recipies)

windowfunction(time, freq)

source code 

Computes the modulus square of the window function of a set of time points at the given frequencies. The time point need not be equidistant. The normalisation is such that 1.0 is returned at frequency 0.

Parameters:
  • time (ndarray) - time points [0..Ntime-1]
  • freq (ndarray) - frequency points. Units: inverse unit of 'time' [0..Nfreq-1]
Returns: array
|W(freq)|^2 [0..Nfreq-1]

check_input(times, signal, **kwargs)

source code 

Check the input arguments for periodogram calculations for mistakes.

If you get an error when trying to compute a periodogram, and you don't understand it, just feed the input you gave to this function, and it will perform some basic checks.

__spread__(y, yy, n, x, m)

source code  

Given an array yy(0:n-1), extirpolate (spread) a value y into
m actual array elements that best approximate the "fictional"
(i.e., possible noninteger) array element number x. The weights
used are coefficients of the Lagrange interpolating polynomial
Arguments:
    y : 
    yy : 
    n : 
    x : 
    m : 
Returns:

getSignificance(wk1, wk2, nout, ofac)

source code 

Returns the peak false alarm probabilities

Hence the lower is the probability and the more significant is the peak

scargle_probability(peak_value, times, freqs, correct_for_frange=False, **kwargs)

source code 

Compute the probability to observe a peak in the Scargle periodogram.

If correct_for_frange=True, the Bonferroni correction will be applied to a smaller number of frequencies (i.e. the independent number of frequencies in freqs). To be conservative, set correct_for_frange=False.

Example simulation:

>>> times = np.linspace(0,1,5000)
>>> N = 500
>>> probs = np.zeros(N)
>>> peaks = np.zeros(N)
>>> for i in range(N):
...     signal = np.random.normal(size=len(times))
...     f,s = scargle(times,signal,threads='max',norm='distribution')
...     peaks[i] = s.max()
...     probs[i] = scargle_probability(s.max(),times,f)

Now make a plot:

>>> p = pl.figure()
>>> p = pl.subplot(131)
>>> p = pl.plot(probs,'ko')
>>> p = pl.plot([0,N],[0.01,0.01],'r-',lw=2)
>>> p = pl.subplot(132)
>>> p = pl.plot(peaks[np.argsort(peaks)],probs[np.argsort(peaks)],'ro')
>>> p = pl.plot(peaks[np.argsort(peaks)],1-(1-np.exp(-np.sort(peaks)))**(10000.),'g-')
>>> #p = pl.plot(peaks[np.argsort(peaks)],1-(1-(1-np.sort(peaks)/2500.)**2500.)**(10000.),'b--')
>>> p = pl.subplot(133)
>>> for i in np.logspace(-3,0,100):
...     p = pl.plot([i*100],[np.sum(probs<i)/float(N)*100],'ko')
>>> p = pl.plot([1e-6,100],[1e-6,100],'r-',lw=2)
>>> p = pl.xlabel('Should observe this many points below threshold')
>>> p = pl.ylabel('Observed this many points below threshold')

]]include figure]]ivs_timeseries_pergrams_prob.png]


   Home       Trees       Indices       Help   
Generated by Epydoc 3.0.1 on Fri Mar 30 10:45:19 2018 http://epydoc.sourceforge.net