Find one frequency, automatically going to maximum precision and return parameters & error estimates.

This routine makes the frequency grid finer until it is finer than the estimated error on the frequency. After that, it will compute (harmonic) parameters and estimate errors.

There is a possibility to escape this optimization by setting scale_df=0 or freqregscale=0.

You can include a nonlinear least square update of the parameters, by setting the keyword optimize=1 (optimization outside loop) or optimize=2 (optimization after each iteration).

The method with which to find the frequency can be set with the keyword method, the model used to fit and optimize should be set with model. Extra keywords for the model functions should go in model_kwargs. If method is a tuple, the first method will be used for the first frequency search only. This could be useful to take advantage of such methods as fasper which do not allow for iterative zoom-ins. By default, the function looks for the highest (or deepest in the case of the pdm method) peak, but instead it is possible to go for the peak having the highest SNR before prewhitening by setting prewhiteningorder_snr to True. In this case, the noise spectrum is calculated using a convolution with a prewhiteningorder_snr_window wide box.

Possible extra keywords: see definition of the used periodogram function.

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

Example keywords:

• 'correlation_correction', default=True
• 'freqregscale', default=0.5: factor for zooming in on frequency
• 'dfscale', default = 0.25: factor for optimizing frequency resolution

Example usage: We generate a sine signal

>>> times = np.linspace(0,100,1000)
>>> signal = np.sin(2*np.pi*2.5*times) + np.random.normal(size=len(times))

Compute the frequency

>>> parameters, pgram, model = find_frequency(times,signal,full_output=True)

Make a plot:

>>> p = pl.figure()
>>> p = pl.axes([0.1,0.3,0.85,0.65])
>>> p = pl.plot(pgram[0],pgram[1],'k-')
>>> p = pl.xlim(2.2,2.8)
>>> p = pl.ylabel('Amplitude')
>>> p = pl.axes([0.1,0.1,0.85,0.2])
>>> p = pl.plot(pgram[0][:-1],np.diff(pgram[0]),'k-')
>>> p = pl.xlim(2.2,2.8)
>>> p,q = pl.xlabel('Frequency (c/d)'),pl.ylabel('Frequency resolution $\Delta F$')

]]include figure]]timeseries_freqanalyse_06.png]

Returns: record array(, 2x1Darray, 1Darray)

parameters and errors(, periodogram, model function)

Fit one or more functions to a timeseries via iterative prewhitening.

This function will use find_frequency to fit the function parameters to the original signal (including any previously found parameters), optimize the parameters if needed, and remove it from the data to search for a new frequency in the residuals.

It is always the original signal that is used to fit all parameters again; only the (optimized) frequency is remembered from step to step (Vanicek's method).

You best set maxiter to some sensable value, to hard-limit the number of frequencies that will be searched for. You can additionally use a stopcrit and stop looking for frequencies once it is reached. stopcrit should be a tuple; the first argument is the function to call, the other arguments are passed to the function, after the mandatory arguments times,signal, modelfunc,allparams,pergram. See stopcrit_scargle_snr for an example of such a function.

By default, the function looks for the highest (or deepest in the case of the pdm method) peak, but instead it is possible to go for the peak having the highest SNR before prewhitening by setting prewhiteningorder_snr to True. In this case, the noise spectrum is calculated using a convolution with a prewhiteningorder_snr_window wide box. Usage of this is strongly encouraged, especially combined with stopcrit_scargle_snr as stopcrit.

Returns: rec array(, ndarray)

parameters, model(, model function)

Compute a 2D periodogram.

Rows (first axis) should be spectra chronologically

x are time points (length N) y are second axis units (e.g. wavelengths) (length NxM)

If the periodogram/wavelength combination has a large range, this script can produce a ValueError or MemoryError. To solve this, you could iterate this function over a subset of wavelength bins yourself, and write the results to a file.

This function also outputs a model, which you can use to subtract from the data (probably together with the average profile, which is also returned) and look further for any frequencies:

   output = spectrum_2D(x,y,matrix,subs_av=True)
   new_matrix = matrix - output['avprof'] - output['model'
   output2 = spectrum_2D(x,y,new_matrix,subs_av=False)

If you want to prewhiten a model, you'd probably want to use the same frequency across the whole line profile. You can hack this by setting f0=frequency and fn=frequency+df with df the size of the frequency step.

Example usage: first we generate some variable line profiles. In this case, this is just a radial velocity variation

>>> times = np.linspace(0,150,100)
>>> wavel = np.r_[4500:4600:1.0]
>>> matrix = np.zeros((len(times),len(wavel)))
>>> for i in xrange(len(times)):
...   central_wave = 5*np.sin(2*np.pi/10*times[i])
...   matrix[i,:] = 1 - 0.5*np.exp( -(wavel-4550-central_wave)**2/10**2)

Once the line profiles are constructed, we can compute the Fourier transform of these lines:

>>> output = spectrum_2D(times,wavel,matrix,method='scargle',model='sine',f0=0.05,fn=0.3,subs_av=True,full_output=True)

With this output, we can find retrieve the frequency. We plot the periodograms for each wavelength bin on the left, and on the right the average over all wavelength bins:

>>> p = pl.figure()
>>> p = pl.subplot(121)
>>> p = pl.imshow(output['pergram'][1][::-1],extent=[wavel[0],wavel[-1],0.05,0.3],aspect='auto')
>>> p = pl.subplot(122)
>>> p = pl.plot(output['pergram'][1].mean(axis=1),output['pergram'][0],'k-')
>>> p = pl.ylim(0.05,0.3)

We can then fix the frequency and compute all the parameters. However, we now choose not to subtract the average profile, but instead use the GLS periodogram to include the constant

>>> output = spectrum_2D(times,wavel,matrix,method='gls',model='sine',f0=0.095,fn=0.105,full_output=True)

>>> p = pl.figure()
>>> p = pl.subplot(221)
>>> p = pl.errorbar(wavel,output['pars']['const'],yerr=output['pars']['e_const'],fmt='ko-')
>>> p = pl.subplot(222)
>>> p = pl.errorbar(wavel,output['pars']['ampl'],yerr=output['pars']['e_ampl'],fmt='ro-')
>>> p = pl.subplot(223)
>>> p = pl.errorbar(wavel,output['pars']['freq'],yerr=output['pars']['e_freq'],fmt='ko-')
>>> p = pl.subplot(224)
>>> p = pl.errorbar(wavel,output['pars']['phase'],yerr=output['pars']['e_phase'],fmt='ro-')

Returns: dict

dict with keys avprof (2D array), pars (rec array), model (1D array), pergram (freqs,2Darray)

Short Time (Fourier) Transform.

Slide a window through the timeseries, multiply the timeseries with the window, and perform a Fourier Transform.

It is best to fix explicitly f0, fn and df, to limit the computation time!

Extra kwargs go to find_frequency

Parameters:

• n_windows (integer) - number of slices
• window_width (float) - width of each slice (defaults to T/20)
• detrend (callable) - detrending function, accepting times and signal as args
• window (string) - window function to apply

Returns: dict

spectrogram, times used, parameters, errors, points used per slice

Decorators:

• @defaults_pergram

Compute the autocorrelation.

The formulate to estimate the autocorrelation in the periodogram is

R(k) = 1 / (n * s^2) * \sum_{t=1}^n (X_t - mu) * (X_{t+k} - mu)

where n is the number of points in the interval, mu an estimator for the sample mean and s is an estimator for the standard deviation of the sample.

When computed over a large interval, we have to take into account that we cannot sample more points than we have in the spectrum, so we have to average out over fewer and fewer points, the more we shift the spectrum.

Above formula becomes

R(k) = 1 / ((n-k) * s^2) * \sum_{t=1}^{n-k} (X_t - mu) * (X_{t+k} - mu)

Highs and lows are cut of. The lower cut off value is the sample mean, the higher cut off is a user-defined multiple of the sample mean.

This function can also be used to compute period spacings: just invert the frequency spacing, reverse the order and reinterpolate to be equidistant: Example:

>>> periods = linspace(1./p[0][-1],1./p[0][0],2*len(p[0]))
>>> ampl = interpol.linear_interpolation(1./p[0][::-1],p[1][::-1],periods)
>>> ampl = where(isnan(ampl),hstack([ampl[1:],ampl[:1]]),ampl)

Parameters:

• frequencies (numpy 1d array) - frequency array
• power (numpy 1d array) - power spectrum
• max_step (float) - maximum time shift
• interval (tuple of floats) - tuple of frequencies (start, end)
• threshold (float) - high cut off value (in units of sample mean)

Returns: (ndarray,ndarray)

domain of autocorrelation and autocorrelation