In [1]:
%matplotlib inline
from numpy import linspace,exp
from numpy.random import randn
from scipy.interpolate import LSQUnivariateSpline,UnivariateSpline
from scipy.optimize import bisect
import matplotlib
from pylab import*
from scipy.fftpack import rfft
#-----------------------------------------------------------------------
fig=figure(1)
ax=fig.add_subplot(111)
M=10
Q=50
n = 2**16
nf = 2**7
ratio = n/nf
kc = 10.
pi = 3.141592
xn = linspace(0,1,n)*2.*pi
xnf = linspace(0,1,nf+1)
xnf2 = linspace(0,2,2*nf+1)
k0 = linspace(0,n/2,n/2+1)
k1 = linspace(0,nf/2,nf/2+1)*(n/nf)
k2 = linspace(0,nf/2,nf/2+1)*(n/nf)
k3 = linspace(0,nf,nf)*(n/(2*nf))
k4 = linspace(0,nf/2,nf/2+1)*(n/nf)
spect = 1./(1.+(k0[:]/kc)**(5./3.))*exp(-0.01*(k0[:]/kc)**1.)
spect[0] = 0.
varo = 2.*trapz(y=spect[:],x=k0[:])
L = spect[1]/varo
print('Integral scale :',L)
print('Size of original sample :',2.*pi/L,' Int. Scale')
print('Size of sub-sampling :',2.*pi*float(nf)/float(n)/L,' Int. Scale')
fftfx = zeros(n/2+1,dtype=complex)
rfx1 = zeros(nf+1,dtype=float)
rfx2 = zeros(nf+1,dtype=float)
rfx3 = zeros(2*nf+1,dtype=float)
rfx4 = zeros(nf+1,dtype=float)
nspect0 = zeros(n/2+1,dtype=float)
nspect1 = zeros(nf/2+1,dtype=float)
nspect2 = zeros(nf/2+1,dtype=float)
nspect3 = zeros(nf,dtype=float)
nspect4 = zeros(nf/2+1,dtype=float)
var0 = 0.
var1 = 0.
var2 = 0.
var3 = 0.
var4 = 0.
for ns in range(0,int(Q)):
ang = rand(n/2)*2.*pi
for i in range(0,int(n/2)):
fftfx[i]=sqrt(spect[i])*complex(cos(ang[i]),sin(ang[i]))
rfx = irfft(fftfx)*float(n)
fftfx0 = rfft(rfx)/float(n)
nspect0[0:n/2] = nspect0[0:n/2] + fftfx0[0:n:2]**2 + fftfx0[1:n:2]**2
var0=var0+var(rfx)
#=================================
for j in range(0,M):
first = int(rand(1)*n/2)
#--------------------------------
rfx1[0:nf+1] = rfx[first:first+nf+1]
var1=var1+var(rfx1)
#--------------------------------
rfx2[0:nf+1] = rfx1[0:nf+1]
for k in range(0,nf+1):
rfx2[k] = rfx1[k] - (rfx1[nf] - rfx1[0])*k/nf - rfx1[0]
var2=var2+var(rfx2)
#--------------------------------
rfx3[0:nf+1] = rfx2[:]
for k in range(0,nf+1):
rfx3[nf+k] = - rfx2[nf-k]
var3=var3+var(rfx3)
#--------------------------------
# window = np.parsen(nf+1)
window = np.hanning(nf+1)
# window = np.blackman(nf+1)
rfx4[0:nf+1] = rfx1[0:nf+1]* window[:]
var4=var4+var(rfx4)
#--------------------------------
# plot(xnf,rfx1,'k-')
# plot(xnf,rfx2,'r-')
# plot(xnf2,rfx3,'b-',marker='.')
# show()
fftfx1 = fft(rfx1[0:nf])/float(nf-1)
fftfx2 = fft(rfx2[0:nf])/float(nf-1)
fftfx3 = fft(rfx3[0:2*nf])/float(2*nf-1)
fftfx4 = fft(rfx4[0:nf])/float(nf-1)
nspect1 = nspect1 + absolute(fftfx1[0:nf/2+1])**2
nspect2 = nspect2 + absolute(fftfx2[0:nf/2+1])**2
nspect3 = nspect3 + absolute(fftfx3[0:nf])**2
nspect4 = nspect4 + absolute(fftfx4[0:nf/2+1])**2
nspect0 = nspect0/float(Q)
nspect1 = nspect1/float(M)/float(Q)/float(ratio)
nspect2 = nspect2/float(M)/float(Q)/float(ratio)
nspect3 = nspect3/float(M)/float(Q)/float(0.5*ratio)
nspect4 = nspect4/float(M)/float(Q)/float(ratio)*2
vars1 = 2.*trapz(y=nspect1[:],x=k1[:])
print('Var 1 spectre (black):', vars1)
var0 = var0/float(Q)
var1 = var1/float(Q)/float(M)
var2 = var2/float(Q)/float(M)
var3 = var3/float(Q)/float(M)
var4 = var4/float(Q)/float(M)
print('Var (green) =',varo)
print('Var 0 (yellow) =',var0)
print('Var 1 (black) =',var1)
print('Var 2 (red) =',var2)
print('Var 3 (blue) =',var3)
print('Var 4 (cyan) =',var4)
plot(k0,nspect0,'y-',label=('Full sample (size '+r'$L$'+')'))
plot(k1,nspect1,'k-',label=('Sub-sample (size '+r'$L_{SUB}$'+')'))
plot(k2,nspect2,'r-',label=('Sub-sample - lin. func.'))
plot(k3,nspect3,'b-',label=('Sub-sample - lin. func. + sym.'))
plot(k4,nspect4,'c-',label=('Sub-sample + Parsen win'))
#plot(k4,nspect4,'c-',label=('Sub-sample + Blackman win'))
#plot(k4,nspect4,'c-',label=('Sub-sample + Hanning win'))
plot(k0,spect,'g-',label=('Theoritical'))
ax.set_xlim(1,n)
ax.set_ylim(1.E-15,1.E+02 )
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel(r"$k$")
ax.set_ylabel( r'$E(k)$')
kfit = array([2.e+01,1.E+03])
yfit = 1.E-00*(kfit/kfit[0])**(-5./3.)
text(1.e+02,2.0E-01,r'$k^{-5/3}$',fontsize=14)
plot(kfit,yfit,'k:')
txt1 = (r"${L/l}=$")+'{:.4f}'.format(2.*pi/L)
text(3.e+03, 1.E+1, txt1 , fontsize=14, fontweight='heavy', color='black', bbox=dict(boxstyle="round,pad=.5",facecolor='white', alpha=1.0))
txt2 = (r"${L_{SUB}/l}=$")+'{:.4f}'.format(2.*pi*float(nf)/float(n)/L)
text(3.e+03, 1.E-1, txt2 , fontsize=14, fontweight='heavy', color='black', bbox=dict(boxstyle="round,pad=.5",facecolor='white', alpha=1.0))
legend(loc='lower left', labelspacing=0, borderpad=0.5,handlelength=3,handletextpad=1,ncol=1)
savefig("test_spectra.eps")
In [ ]: