Part c¶

In [13]:

hi = 1250

data = load('/home/bhardin2/ece458/2015.01.08.14.02.24.MST_10min.npz')
dat = data['MSTdata']

M = 150
n = 32

F2 = zeros(n)
for j in range(M):
    Estub = dat[n*j:n*(j+1),hi]
    F = fft.fft(Estub)/n
    F2 += abs(F)**2
F2 = F2/M
F2s = fft.fftshift(F2)
x = arange(n)
xs = fftshift(x)
xs[xs>=n/2] = xs[xs>=n/2] - n
plot(xs,F2s,'k')
xlabel('FFT bin')
ylabel('Magnitude-squared')
axis('tight')

Out[13]:

(-16.0, 15.0, 15165.447642751795, 21829.398021851859)

In [26]:

# Use an unbiased estimator of the grass level from the 32 observations (std(F2) is a biased estimator)
sqrt(1./(N-1) * sum(dF**2))

Out[26]:

1616.798152045453

The estimated grass level $\langle \delta S^2\rangle^{1/2}$ is 1617, which is relatively close to the expected value of 1530.

Part d¶

In [30]:

Mvec = [10,50,100,150]
fmt  = ['k','k.-','k--','k-.']
grass = zeros(len(Mvec))

for k in range(len(Mvec)):
    M = Mvec[k]
    
    F2 = zeros(n)
    for j in range(M):
        Estub = dat[n*j:n*(j+1),hi]
        F = fft.fft(Estub)/n
        F2 += abs(F)**2
    F2 = F2/M
    F2s = fft.fftshift(F2)
    x = arange(n)
    xs = fftshift(x)
    xs[xs>=n/2] = xs[xs>=n/2] - n
    plot(xs,F2s,fmt[k],label='M=%i' % M)
    
    dF = F2 - mean(F2)
    grass[k] = sqrt(1./(N-1) * sum(dF**2))
    
xlabel('FFT bin')
ylabel('Magnitude-squared')
axis('tight')
legend(loc='best')

Out[30]:

<matplotlib.legend.Legend at 0x7ff608bc5bd0>

In [40]:

figure(figsize=(4,4))
loglog(Mvec,grass,'ko-')
grid()
axis('equal')
xlabel('I')
ylabel('Grass level')

Out[40]:

<matplotlib.text.Text at 0x7ff6081dffd0>

It looks like the grass levels go down as expected, as a square-root relationship. In a log-log plot, the slope should be $1/2$, and it looks like the slope is indeed close to $1/2$.