Source Code for Module ivs.timeseries.eacf

  1  # -*- coding: utf-8 -*- 
  2  """ 
  3  Computes the envelope auto-correlation function. This function is estimated by computing the power spectrum of  
  4  a smoothed power spectrum of a time series. It is ideally suited to discover periodicities in the power spectrum, 
  5  as is often the case for solar-like oscillators (cf the asymptotic relation of Tassoul), and can thus be used to  
  6  estimate the mean large separation (mean frequency spacing between two consecutive radial modes). 
  7   
  8  As an example we take the Kepler red giant KIC3744043, which shows a beautiful spectrum of solar-like oscillations: 
  9   
 10  ]include figure]]powerspec_KIC3744043.png] 
 11   
 12  To compute the EACF, we first import: 
 13   
 14  >>> import numpy as np 
 15  >>> from ivs.timeseries.eacf import eacf 
 16   
 17  Red giants have spacings in the power spectrum between, say, 1.0 and 15.0, so this is the interval in which we will 
 18  compute the EACF: 
 19   
 20  >>> spacings = np.arange(1.0, 15.0, 0.005) 
 21   
 22  The figure above shows that the oscillations are only present between (roughly) 80 microHz and 135 microHz, so we will limit ourselves to this frequency interval: 
 23   
 24  >>> minFreq, maxFreq = 80.0, 135.0 
 25   
 26  The EACF is now computed as: 
 27   
 28  >>> autoCorrelation, croppedFreqs, smoothedSpectrum = eacf(freq, spec, spacings, 2.0, minFreq, maxFreq) 
 29   
 30  The value '2.0' is related to the FWHM of the (Hann) kernel that is used to smooth the power spectrum of the signal, on which we will come back later. 
 31   
 32  The EACF looks as follows: 
 33   
 34  ]include figure]]eacf_smooth15.png] 
 35   
 36  The highest peak is caused by the main periodicity in the power spectrum. This does not, however, correspond directly  
 37  to the mean large separation, but to half the large separation. The reason is that almost in the middle between each  
 38  two radial modes there is a dipole mode. To compute the large separation we therefore have to enter: 
 39   
 40  >>> meanLargeSeparation = 2.0 * spacings[np.argmax(autoCorrelation)] 
 41  >>> meanLargeSeparation 
 42  9.8599999999998325 
 43   
 44  How did we choose the value 2.0 (microHz) for the FWHM of the smoothing kernel? The smoothing is done by convolving 
 45  the power spectrum with a Hann kernel with a certain width. If you take the width too large, the spectrum will be too 
 46  heavily smoothed so that no periodicities can be detected. If you take the width too small, the EACF will contain too  
 47  many features, which makes it difficult to recognize the right peak. As a rule-of-thumb you can take the kernel width 
 48  to be 1/8 of your initial estimate of the mean large separation. If you don't have such an initial estimate, you can 
 49  estimate nu_max, which is the frequency of maximum power (= the location of the gaussian envelope of the  
 50  power excess), and derive from that value a first estimate for the large separation: 
 51   
 52  >>> nuMax = 120.0                                      # for KIC3744043 
 53  >>> initialEstimateLargeSep = 0.24 * nuMax**0.78       # general approximation for red giants 
 54  >>> kernelWidth = initialEstimateLargeSep / 8.0        # rule of thumb 
 55  >>> autoCorrelation, croppedFreqs, smoothedSpectrum = eacf(freq, spec, spacings, kernelWidth, minFreq, maxFreq) 
 56   
 57  """ 
 58   
 59   
 60  __all__ = ['eacf'] 
 61   
 62   
 63  import numpy as np 
 64  from ivs.timeseries.pergrams import DFTpower2 
 65   
 66 -def eacf(freqs, spectrum, spacings, kernelWidth, minFreq=None, maxFreq=None, doSanityCheck=False): 
 67   
 68      """ 
 69      Compute the Envelope Auto-Correlation Function (EACF) of a signal. 
 70   
 71      The EACF is derived by computing the power spectrum of a smoothed version of the  
 72      power spectrum of the time series. It allows to identify equispaced patterns in the  
 73      power spectrum 
 74   
 75      Source:  
 76          - Mosser & Appourchaux, 2009, A&A 508, p. 877 
 77          - Mosser, 2010, Astron. Nachr. 331, p. 944 
 78       
 79      @param freqs: frequencies in which the power specturm is given. Assumed to be equidistant. 
 80      @type freqs: ndarray 
 81      @param spectrum: power spectrum corresponding to the given frequency points. 
 82      @type spectrum: ndarray 
 83      @param spacings: array of frequency spacings for which the EACF needs to be computed 
 84      @type spacings: ndarray 
 85      @param kernelWidth: the width (in freq units) of the convolution kernel used to first smooth the power spectrum 
 86      @type kernelWidth: float 
 87      @param minFreq: only the part [minFreq, maxFreq] of the spectrum is taken into account. 
 88                      If it is None, the 1st point of the 'freqs' array is taken. 
 89      @type minFreq: float 
 90      @param maxFreq: only the part [minFreq, maxFreq] of the spectrum is taken into account 
 91                      If it is None, the last point of the 'freqs' array is taken. 
 92      @type maxFreq: float 
 93      @param doSanityCheck: if True: make a few sanity checks of the arguments (e.g. equidistancy of freqs). 
 94                            If False, do nothing. 
 95      @type doSanityCheck: boolean 
 96      @return: autoCorrelation, croppedFreqs, smoothedSpectrum 
 97                  - autoCorrelation:  the EACF evaluated in the values of 'spacings' 
 98                  - croppedFreqs:     the frequencies of the selected part [minFreq, maxFreq] 
 99                  - smoothedSpectrum: the smoothed power spectrum used to compute the EACF. Same length as croppedFreqs.     
100      @rtype: (ndarray, ndarray, ndarray) 
101      """ 
102       
103      # If requested, perform sanity checks 
104      # The 1000. in the check on equidistancy was needed to let pass the arrays created by  
105      # linspace() and arange(). 
106   
107      if doSanityCheck: 
108          if len(freqs) != len(spectrum): 
109              raise ValueError("freqs and spectrum don't have the same length") 
110          if np.fabs(np.diff(np.diff(freqs)).max()) > 1000. * np.finfo(freqs.dtype).eps: 
111              raise ValueError("freqs array is not equidistant") 
112          if np.sometrue(spacings <= 0.0): 
113              raise ValueError("spacings are not all strictly positive") 
114          if kernelWidth <= 0.0: 
115              raise ValueError("kernel width is not > 0.0") 
116          if minFreq > freqs[-1]: 
117              raise ValueError("minimum frequency 'minFreq' is larger than largest frequency in freqs array") 
118          if maxFreq < freqs[0]: 
119              raise ValueError("maximum frequency 'maxFreq' is smaller than smallest frequency in freqs array") 
120          if minFreq >= maxFreq: 
121              raise ValueError("minFreq >= maxFreq") 
122   
123   
124      # Set the default values 
125       
126      if minFreq == None: minFreq = freqs[0] 
127      if maxFreq == None: maxFreq = freqs[-1] 
128      freqStep = freqs[1]-freqs[0] 
129      
130      # Crop the spectrum to the specified range 
131       
132      croppedSpectrum = spectrum[(freqs >= minFreq) & (freqs <= maxFreq)] 
133      croppedFreqs = freqs[(freqs >= minFreq) & (freqs <= maxFreq)] 
134       
135      # Set up a normalized Hann smoothing kernel. 
136      # Note: a Hann window of size N has FWHM = (N-1)/2. The FWHM is given by  
137      # the user in frequency units, so the corresponding N can be derived. 
138        
139      kernelSize = 1 + 2 * int(round(kernelWidth/freqStep)) 
140      kernel = np.hanning(kernelSize) 
141      kernel = kernel/kernel.sum() 
142       
143      # Smooth the spectrum using a convolution 
144       
145      smoothedSpectrum = np.convolve(kernel, croppedSpectrum, mode='same') 
146      smoothedSpectrum -= smoothedSpectrum.mean() 
147       
148      # Compute the power spectrum of the power spectrum 
149       
150      autoCorrelation = DFTpower2(croppedFreqs, smoothedSpectrum, 1.0/spacings) 
151       
152      # That's it. 
153       
154      return autoCorrelation, croppedFreqs, smoothedSpectrum 
155