脚本参考：连续功率谱估计。

qxsw2016

按教课书分析连续谱估计：
#打开1月上旬温度数据文件
table = readtable('f:/temp.csv', delimiter=',', format='%s%f')
xt = table['temp']
#print xt
n = len(xt)
m = 9 #n/2-n/3
rtao = zeros(m+1)
sl = zeros(m+1)
#step1
for tao in arange(0,m+1):
    x1 = xt[0:n-tao-1]
    x2 = xt[tao:n-1]
    slope, intercept, r, p, std,number = linregress(x1, x2)
    rtao[tao] = r
#step3     
for l in arange(0,m+1):
    sl[l] = rtao[0]
    for tao in arange(1,m):
           sl[l] = sl[l]+rtao[tao]*(1+cos(pi*tao/m))*cos(pi*l*tao/m)
    if(l == 0 or l == m):        
        sl[l] = 0.5*sl[l]/m
    else:
        sl[l] = sl[l]/m
print sl
x = arange(0,m+1)
xtop = zeros(9)
xtop = 18.0/x  
ax1=axes()
plot(x,sl,'black',linewidth=2)   
xticks(x)   
xlabel('t',fontsize=16)
ylabel(r'S(l)',fontsize=16)
xlim(0,9)
ax2 = twiny(ax1)
xaxis(location='top',tickin=False)
xlabel('T(year)',bold=False)
xlim(0,9)
xticks(x, ['%.1f' % xx for xx in xtop])
title('Continuous power spectrum estimation')

yinwangyi01

谢谢分享，很有用，感谢！！

qxsw2016

如果序列的滞相关系数rtao(1)为较大正值是，表明序列具有持续性，用红噪声标准标谱检验，参考代码如下：
#红噪声标准谱alpha = 0.05
sm =  (sl[0]+sl[m])/(2*m)+sum(sl[1:m-1])/m #m+1个谱估计值得均值
nju = (2*n-m/2.)/m
for k in arange(0,m+1):
    s0k[k] = sm*((1-rtao[1]**2)/(1+rtao[1]**2-2*rtao[1]*cos(pi*k/m)))
kai2 = 11.5  #查表估计值
s0k= s0k*kai2/nju

qxsw2016

最大熵谱脚本参考：
以奇异谱分析的重建成分为基础，基于Burg 算法的最大熵谱估计（数组索引开始位置为1，有点晕），过程中用到预测误差准则FPE确定最佳回归阶数。rc3+rc4的最大熵谱分析如下：
fn= r'F:\dataset\rc34.txt'
ncol = numasciicol(fn)
nrow =  numasciirow(fn)
fndata = asciiread(fn,shape=(nrow,ncol))
x = fndata[:,0]
n = 49
m = 24 #最大波数
delta2 = zeros(n+1)
a = zeros((n+1,n+1)) #处理方便，递推公式索引从1开始
r = zeros(n+1)
fpe = zeros(n-1)
shl = zeros(m+1)
#基于Burg 算法的最大熵谱估计
# k = 0
qsum = sum(x*x)
delta2[0] = qsum/float(n)
r[0]= delta2[0]
# k = 1
s12 = sum(x[0:n-1]*x[1:n])
s11 = sum(x[0:n-1]*x[0:n-1])
s22 = sum(x[1:n]*x[1:n])

a[1,1] = 2*s12/(s11+s22)        
delta2[1]=(1-a[1,1]*a[1,1])*delta2[0]
r[1] = a[1,1]*delta2[0]
# k = 2
x1 = x[2:n]-a[1,1]*x[1:n-1]
x2 = x[0:n-2]-a[1,1]*x[1:n-1]
s12 = sum(x1*x2)
s11 = sum(x1*x1)
s22 = sum(x2*x2)

a[2,2] = 2*s12/(s11+s22)
a[2,1] = a[1,1]-a[2,2]*a[1,1]
delta2[2] = (1-a[2,2]*a[2,2])*delta2[1]
r[2] = a[1,1]*r[1]+a[2,2]*delta2[1]
# k >2
for k in arange(2,n-1):      
    s12= 0   
    s11 = 0
    s22 = 0         
    for i in arange(k,n):
        v1 =0
        v2 =0
        for j in arange(1,k+1):
            v1 = v1 + a[k,j]*x[i-j]
            v2 = v2 + a[k,j]*x[i-k+j-1]            
        v1 =x - v1  
        v2 =x[i-k-1] - v2              
        s12 = s12 + 2*v1*v2   
        s11 = s11 + v1*v1
        s22 = s22 + v2*v2
    a[k+1,k+1] = s12/(s11+s22)

    for j in arange(1,k+1):
        a[k+1,j] = a[k,j]-a[k+1,k+1]*a[k,k+1-j]
    delta2[k+1] = (1-a[k+1,k+1]*a[k+1,k+1])*delta2[k]

    v1 = sum(a[k,1:k+1]*r[k+1:1:-1])
    r[k+1] = v1+a[k+1,k+1]*delta2[k]
fpe[0] = 999999
#fpe 准则确定最佳回归阶数k0
for k in arange(1,n):
    fpe[k-1] = (n+k)*delta2[k]/(n-k)
k0 = fpe.index(min(fpe))
print k0
print delta2[11]
for l in arange(0,m+1):
    v1 = 0
    v2 = 0
    for k in arange(1,k0+1):
        v1 = v1 + a[k0,k]*cos(pi*l*k/m)
        v2 = v2 + a[k0,k]*sin(pi*l*k/m)
    v1 = (1-v1)*(1-v1)
    v2 = v2*v2
    shl[l] = delta2[k0]/(v1+v2)
xtick = arange(0,m+1)
xtop = zeros(m+1)
xtop[1:m+1] = 2*m*1./xtick[1:m+1]  

ax1=axes()
plot(xtick,shl,'blacks--',linewidth=2)   
xticks(xtick)   
xlabel('L',fontsize=16)
xlim(0,m+1)
ax2 = twiny(ax1)
xaxis(location='top',tickin=False)
xlabel('T(year)',bold=False)
xlim(0,m+1)
xticks(xtick, ['%.1f' % xx for xx in xtop])
title('RC3+RC4 maximum entropy spectra')

qxsw2016

连续谱分析，编了一个函数，原来的落后自相关系数计算与书籍上有差异，但仅在计算数值上有微小的系统差异，最终谱估计的结果也基本一致，另外增加了功率谱检验。
def compute_acf(x, max_lag=None):
   
    if max_lag is None:
        max_lag = len(x) // 3
    n = len(x)
   
    x_mean = np.mean(x)
    variance = np.sum((x - x_mean)**2) / n   
   
    r = np.zeros(max_lag + 1)      
    r[0] = 1.0   
   
    for tau in range(1, max_lag + 1):        
        numerator = np.sum((x[:n-tau] - x_mean) * (x[tau:] - x_mean)) / n        
        denominator = variance      
        r[tau] = numerator / denominator
        
    return r

def spectrumAna(x,m,alpha1,alpha2,isDetrend=False):

    n = x.shape[0]
    x_detrend = np.signal.detrend(x, type='linear')

    x_used = x
    if(isDetrend):
        x_used = x_detrend  

    rtau = compute_acf(x_used, max_lag=m)
   
    sl = np.zeros(m+1)
    l_vals = arange(m+1)
    periods = 2 * m / l_vals
   
    for l in range(0, m+1):
        sum_val = rtau[0]
        for tau in range(1, m):
            window = 1 + np.cos(pi * tau / m)
            sum_val += rtau[tau] * window * np.cos(pi * l * tau / m)
        if l == 0 or l == m:
            sl[l] = 0.5 * sum_val / m
        else:
            sl[l] = sum_val / m
   

    slm = mean(sl)
    r1 = rtau[1]

    dof = (2 * n - m / 2.0) / m
   
    chi2 = stats.chi2.ppf(1-alpha1, dof)
    print chi2
    sol = slm * (1 - r1**2) / (1 + r1**2 - 2 * r1 * cos(np.pi * l_vals / m))
    sred = sol*chi2/dof

    chi2 = stats.chi2.ppf(1-alpha2, dof)
    swhite = ones(len(l_vals)) * slm*chi2/dof

    return l_vals,periods,sl,sred,swhite,r1

stfn = 'temp.csv'
df= DataFrame.read_table(stfn, delimiter=',', format='%i%3f')
x = df['temp'].values
l_vals,periods,sl,sred,swhite,r1 = spectrumAna(x,9,0.05,0.05,isDetrend=False)

ax1=axes()
plot(l_vals,sl,'black',linewidth=2)   
plot(l_vals, sred, 'r--', linewidth=2)
plot(l_vals, swhite, 'm--', linewidth=2)

ax2 = twiny(ax1)
xaxis(location='top',tickin=False)
xticks(x[1:], ['%.1f' % xx for xx in periods[1:]])