用LombScargle方法进行峰值检测

""" from numpy import * from numpy.fft import * def __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. The weights used are coefficients of the Lagrange interpolating polynomial Arguments: y : yy : n : x : m : Returns: """ nfac=[0,1,1,2,6,24,120,720,5040,40320,362880] if m > 10. : print 'factorial table too small in spread' return ix=long(x) if x == float(ix): yy[ix]=yy[ix]+y else: ilo = long(x-0.5*float(m)+1.0) ilo = min( max( ilo , 1 ), n-m+1 ) ihi = ilo+m-1 nden = nfac[m] fac=x-ilo for j in range(ilo+1,ihi+1): fac = fac*(x-j) yy[ihi] = yy[ihi] + y*fac/(nden*(x-ihi)) for j in range(ihi-1,ilo-1,-1): nden=(nden/(j+1-ilo))*(j-ihi) yy[j] = yy[j] + y*fac/(nden*(x-j)) def fasper(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). 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) """ #Check dimensions of input arrays n = long(len(x)) if n != len(y): print 'Incompatible arrays.' return nout = 0.5*ofac*hifac*n nfreqt = long(ofac*hifac*n*MACC) #Size the FFT as next power nfreq = 64L # of 2 above nfreqt. while nfreq < nfreqt: nfreq = 2*nfreq ndim = long(2*nfreq) #Compute the mean, variance ave = y.mean() ##sample variance because the divisor is N-1 var = ((y-y.mean())**2).sum()/(len(y)-1) # and range of the data. xmin = x.min() xmax = x.max() xdif = xmax-xmin #extirpolate the data into the workspaces wk1 = zeros(ndim, dtype='complex') wk2 = zeros(ndim, dtype='complex') fac = ndim/(xdif*ofac) fndim = ndim ck = ((x-xmin)*fac) % fndim ckk = (2.0*ck) % fndim for j in range(0L, n): __spread__(y[j]-ave,wk1,ndim,ck[j],MACC) __spread__(1.0,wk2,ndim,ckk[j],MACC) #Take the Fast Fourier Transforms wk1 = ifft( wk1 )*len(wk1) wk2 = ifft( wk2 )*len(wk1) wk1 = wk1[1:nout+1] wk2 = wk2[1:nout+1] rwk1 = wk1.real iwk1 = wk1.imag rwk2 = wk2.real iwk2 = wk2.imag df = 1.0/(xdif*ofac) #Compute the Lomb value for each frequency hypo2 = 2.0 * abs( wk2 ) hc2wt = rwk2/hypo2 hs2wt = iwk2/hypo2 cwt = sqrt(0.5+hc2wt) swt = sign(hs2wt)*(sqrt(0.5-hc2wt)) den = 0.5*n+hc2wt*rwk2+hs2wt*iwk2 cterm = (cwt*rwk1+swt*iwk1)**2./den sterm = (cwt*iwk1-swt*rwk1)**2./(n-den) wk1 = df*(arange(nout, dtype='float')+1.) wk2 = (cterm+sterm)/(2.0*var) pmax = wk2.max() jmax = wk2.argmax() #Significance estimation #expy = exp(-wk2) #effm = 2.0*(nout)/ofac #sig = effm*expy #ind = (sig > 0.01).nonzero() #sig[ind] = 1.0-(1.0-expy[ind])**effm #Estimate significance of largest peak value expy = exp(-pmax) effm = 2.0*(nout)/ofac prob = effm*expy if prob > 0.01: prob = 1.0-(1.0-expy)**effm return wk1,wk2,nout,jmax,prob def getSignificance(wk1, wk2, nout, ofac): """ returns the peak false alarm probabilities Hence the lower is the probability and the more significant is the peak """ expy = exp(-wk2) effm = 2.0*(nout)/ofac sig = effm*expy ind = (sig > 0.01).nonzero() sig[ind] = 1.0-(1.0-expy[ind])**effm return sig

13.5945121951 13.5945121951 12.6615853659 12.6615853659 12.6615853659 4.10975609756 4.10975609756 4.10975609756 7.99695121951 7.99695121951 16.237804878 16.237804878 16.237804878 16.0823170732 16.237804878 16.237804878 8.92987804878 8.92987804878 10.6402439024 10.6402439024 28.0548780488 28.0548780488 28.0548780488 27.8993902439 27.8993902439 41.5823170732 41.5823170732 41.5823170732 41.5823170732 41.5823170732 41.5823170732 18.7256097561 15.9268292683 15.9268292683 15.9268292683 15.9268292683 15.9268292683 15.9268292683 14.0609756098 14.0609756098 14.0609756098 14.0609756098 14.0609756098 23.8567073171 23.8567073171 23.8567073171 23.8567073171 25.4115853659 25.4115853659 28.0548780488 40.0274390244 40.0274390244 40.0274390244 40.0274390244 40.0274390244 40.0274390244 20.5914634146 20.5914634146 20.4359756098 19.6585365854 18.2591463415 19.3475609756 18.2591463415 10.3292682927 27.743902439 27.743902439 27.743902439 27.743902439 27.743902439 27.743902439 22.3018292683 22.3018292683 21.368902439 21.368902439 21.368902439 21.5243902439 20.4359756098 20.4359756098 20.4359756098 20.4359756098 20.4359756098 20.4359756098 20.4359756098 11.8841463415 11.8841463415 1.0 11.1067073171 10.1737804878 14.5274390244 14.5274390244 14.5274390244 14.5274390244 14.5274390244 14.5274390244 11.7286585366 11.7286585366 12.6615853659 11.7286585366 8.15243902439 1.0 7.84146341463 6.90853658537 12.6615853659 12.6615853659 12.6615853659 12.6615853659 12.6615853659 12.6615853659 12.6615853659 12.6615853659 12.6615853659 13.1280487805 12.9725609756 12.9725609756 12.9725609756 10.3292682927 10.3292682927 10.3292682927 10.3292682927 9.55182926829 10.4847560976 29.9207317073 29.9207317073 29.9207317073 29.9207317073 30.0762195122 30.0762195122 26.1890243902 7.99695121951 25.256097561 7.99695121951 7.99695121951 7.99695121951 6.59756097561 6.59756097561 6.59756097561 6.59756097561 7.53048780488 7.53048780488 7.53048780488 7.53048780488 7.53048780488 7.53048780488 7.53048780488 7.53048780488 10.0182926829 10.0182926829 10.0182926829 10.0182926829 10.0182926829 10.0182926829 10.4847560976 15.9268292683 15.9268292683 15.9268292683 15.9268292683 15.9268292683 16.8597560976 15.9268292683 15.9268292683 16.8597560976 16.7042682927 16.7042682927 16.7042682927 9.08536585366 8.46341463415 8.46341463415 8.46341463415 8.46341463415 6.90853658537 7.84146341463 6.90853658537 4.26524390244 12.3506097561 12.3506097561 12.3506097561 12.3506097561 12.3506097561 12.3506097561 12.3506097561 12.3506097561 12.3506097561 12.3506097561 12.3506097561 14.2164634146 14.2164634146 14.2164634146 14.0609756098 14.0609756098 14.0609756098 14.0609756098 16.8597560976 16.8597560976 16.7042682927 16.7042682927 16.7042682927 16.7042682927 17.9481707317 17.9481707317 19.6585365854 19.6585365854 19.6585365854 19.6585365854 10.7957317073 10.7957317073 10.7957317073 10.7957317073 10.7957317073 12.1951219512 12.1951219512 22.9237804878 22.9237804878 22.9237804878 22.9237804878 22.9237804878 22.9237804878 22.9237804878 7.84146341463 7.84146341463 7.84146341463 7.84146341463 8.7743902439 8.7743902439 7.84146341463 8.61890243902 8.61890243902 8.61890243902 8.61890243902 18.2591463415 18.2591463415 18.2591463415 18.2591463415 18.2591463415 18.2591463415 18.2591463415 18.2591463415 18.2591463415 9.39634146341 9.39634146341 9.24085365854 9.24085365854 9.24085365854 9.24085365854 9.08536585366 9.08536585366 9.08536585366 9.08536585366 9.55182926829 9.55182926829 9.55182926829 9.55182926829 9.55182926829 16.5487804878 16.5487804878 16.5487804878 16.5487804878 16.5487804878 16.5487804878 16.5487804878 16.5487804878 16.5487804878 16.5487804878 16.5487804878 16.5487804878 16.5487804878 16.5487804878 1.0 16.0823170732 16.0823170732 16.0823170732 16.0823170732 16.0823170732 16.0823170732 16.0823170732 16.0823170732 16.0823170732 17.1707317073 17.0152439024 21.9908536585 21.9908536585 21.9908536585 21.9908536585 21.9908536585 21.9908536585 21.9908536585 7.84146341463 8.7743902439 7.84146341463 6.75304878049 5.9756097561 5.9756097561 5.9756097561 5.9756097561 5.9756097561 5.9756097561 3.95426829268 7.06402439024 7.06402439024 7.06402439024 11.262195122 11.262195122 11.262195122 11.262195122 11.262195122 11.262195122 9.08536585366 9.86280487805 7.99695121951 7.99695121951 14.2164634146 14.0609756098 14.0609756098 14.0609756098 14.0609756098 14.0609756098 2.24390243902 2.08841463415 3.02134146341 3.02134146341 2.08841463415 4.73170731707 4.73170731707 4.73170731707 4.73170731707 6.44207317073 6.44207317073 6.44207317073 6.44207317073 6.44207317073 6.44207317073 6.44207317073 6.44207317073 6.44207317073 6.44207317073 6.59756097561 6.59756097561 6.59756097561 6.75304878049 1.0 6.28658536585 6.28658536585 7.21951219512 6.28658536585 10.6402439024 10.6402439024 10.6402439024 10.6402439024 10.6402439024 10.6402439024 10.6402439024 14.3719512195 14.3719512195 15.6158536585 15.6158536585 15.6158536585 35.6737804878 35.6737804878 35.6737804878 35.6737804878 35.6737804878 35.6737804878 35.6737804878 35.6737804878 35.6737804878 35.6737804878 35.6737804878 28.6768292683 28.6768292683 28.6768292683 28.6768292683 28.6768292683 51.8445121951 51.8445121951 51.8445121951 51.8445121951 51.8445121951 52.0 52.0 4.42073170732 4.42073170732 5.9756097561 5.9756097561 5.9756097561 5.9756097561 5.9756097561 5.9756097561 4.10975609756 3.95426829268 3.64329268293 3.64329268293 4.73170731707 4.73170731707 6.28658536585 6.28658536585 6.28658536585 6.28658536585 6.28658536585 6.28658536585 6.28658536585 5.9756097561 5.82012195122 5.82012195122 5.82012195122 5.82012195122 5.82012195122 12.1951219512 12.1951219512 12.1951219512 12.1951219512 12.1951219512 12.1951219512 12.1951219512 12.1951219512 1.0 11.7286585366 11.7286585366 11.7286585366 11.7286585366 11.7286585366 11.7286585366 11.1067073171 11.1067073171 11.1067073171 11.1067073171 11.1067073171 11.1067073171 11.1067073171 11.1067073171 10.0182926829 10.0182926829 16.7042682927 16.7042682927 16.7042682927 16.7042682927 16.7042682927 16.7042682927 29.1432926829 29.1432926829 29.1432926829 29.1432926829 29.1432926829 29.1432926829 29.1432926829 29.1432926829 29.1432926829 1.15548780488 2.71036585366 2.71036585366 2.71036585366 2.71036585366 2.71036585366 2.71036585366 2.71036585366 3.17682926829 4.10975609756 4.10975609756 5.9756097561 5.9756097561 5.9756097561 6.90853658537 5.9756097561 10.1737804878 10.1737804878 10.1737804878 8.61890243902 8.46341463415 8.46341463415 9.39634146341 8.46341463415 8.46341463415 5.35365853659 5.35365853659 5.35365853659 5.35365853659 5.35365853659 5.35365853659 3.33231707317 4.42073170732 3.33231707317 6.59756097561 6.44207317073 5.82012195122 6.75304878049 5.82012195122 5.82012195122 5.82012195122 4.73170731707 5.66463414634 5.66463414634 4.73170731707 4.73170731707 5.66463414634 5.66463414634 5.50914634146 2.71036585366 5.50914634146 2.71036585366 2.71036585366 5.50914634146 5.50914634146 5.50914634146 6.28658536585 6.28658536585 5.9756097561 5.9756097561 7.06402439024 5.9756097561 7.53048780488 8.46341463415 8.46341463415 13.2835365854 13.2835365854 13.2835365854 13.2835365854 2.55487804878 2.55487804878 2.55487804878 2.55487804878 4.10975609756 3.17682926829 3.17682926829 4.26524390244 3.64329268293 3.64329268293 3.64329268293 3.33231707317 3.33231707317 3.33231707317 2.24390243902 3.33231707317 2.24390243902 2.24390243902 3.64329268293 3.64329268293 3.64329268293 3.64329268293 3.64329268293 3.64329268293 7.53048780488 7.53048780488 7.53048780488 7.53048780488 7.53048780488 7.53048780488 7.53048780488 7.53048780488 7.53048780488 6.28658536585 6.28658536585 7.21951219512 6.28658536585 6.28658536585 6.28658536585 6.28658536585 6.28658536585 6.28658536585 3.7987804878 4.73170731707 3.7987804878 3.7987804878 3.7987804878 3.7987804878 3.7987804878 3.7987804878 4.26524390244 4.26524390244 5.19817073171 5.19817073171 6.28658536585 6.28658536585 6.28658536585 6.28658536585 6.28658536585 6.28658536585 6.28658536585 6.28658536585 7.53048780488 7.53048780488 7.53048780488 7.53048780488 7.53048780488 7.53048780488 3.7987804878 3.7987804878 3.95426829268 3.02134146341 3.02134146341 3.02134146341 1.0 1.93292682927 2.55487804878 2.55487804878 5.9756097561 5.9756097561 5.9756097561 5.9756097561 5.9756097561 5.9756097561 5.9756097561 5.9756097561 5.9756097561 5.9756097561 5.9756097561 6.28658536585 6.28658536585 6.28658536585 6.28658536585 6.28658536585 6.28658536585 16.0823170732 16.0823170732 31.3201219512 31.3201219512 31.3201219512 31.3201219512 31.3201219512 31.3201219512 31.3201219512 31.3201219512 3.64329268293 3.64329268293 4.26524390244 4.26524390244 3.7987804878 4.73170731707 3.7987804878 3.7987804878 2.55487804878 3.48780487805 2.55487804878 2.55487804878 3.17682926829 3.17682926829 3.17682926829 3.17682926829 3.17682926829 3.17682926829 3.17682926829 3.17682926829 3.17682926829 3.17682926829 3.17682926829 3.17682926829 3.17682926829 3.17682926829 3.17682926829 3.17682926829 3.17682926829 3.17682926829 3.17682926829 3.17682926829 3.33231707317 12.3506097561 12.3506097561 12.3506097561 12.3506097561 12.3506097561 12.3506097561 4.73170731707 4.73170731707 4.73170731707 4.73170731707 4.73170731707 4.73170731707 4.73170731707 4.73170731707 2.86585365854 2.86585365854 1.46646341463 1.46646341463 1.46646341463 1.46646341463 1.46646341463 1.46646341463 1.62195121951 1.62195121951 1.62195121951 1.77743902439 1.77743902439 4.42073170732 4.42073170732 4.42073170732 4.42073170732 4.42073170732 4.42073170732 4.42073170732 3.95426829268 3.95426829268 2.71036585366 2.71036585366 2.71036585366 2.71036585366 2.71036585366 1.77743902439 2.86585365854 3.02134146341 2.86585365854 2.86585365854 3.17682926829 3.17682926829

1条回答

网友

1楼 · 发布于 2024-05-15 10:30:20

经过一番挖掘，看来AstroML方法是最好的。在

import numpy as np
from matplotlib import pyplot as plt
from astroML.time_series import lomb_scargle, search_frequencies
import pandas as pd

df = pd.read_csv('extinct.csv',header=None)
Y = df[0]

dy = 0.5 + 0.5 * np.random.random(len(df))
omega = np.linspace(10, 100, 1000)

sig = np.array([0.1, 0.01, 0.001])
PS, z = lomb_scargle(df.index, Y, dy, omega, generalized=True, significance=sig)

plt.plot(omega,PS)
plt.hold(True)

xlim = (omega[0], omega[-1])
for zi, pi in zip(z, sig):
    plt.plot(xlim, (zi, zi), ':k', lw=1)
    plt.text(xlim[-1] - 0.001, zi - 0.02, "$%.1g$" % pi, ha='right', va='top')
    plt.hold(True)
plt.show()

这给了

显著性水平也显示在图表上。我以前用广义最小二乘法，不使用平滑法。在

相关问题更多 >

编程相关推荐

热门问题

热门文章