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

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')