%this is a matlab program to calculate and plot spectra of red

% noise.

% It shows how samples from red noise gradually converge from

% spikey spectra to smooth as the sample size or number of realizations

% is increased.

% For the basic demonstration, start the program and just press the return button

% using the default values, and watch the graphic change. Then play with the parameters.

% You need to subroutines inputwd.m and fstat.m

n = inputwd('What is length of window? (256)',256)

number = inputwd('Maximum number of realizations? (40)',40)

nplot = inputwd('Number of realizations between plots? (1)',1)

% Generate red-noise autocorrelation

alpha = inputwd('What is one-lag autocorrelation of red noise? (0.5)',0.5)

frac = inputwd('What fraction of frequencies to plot? (1.0)',1.0)

seeder = inputwd('What is seed for random number generator? (3)',3.0)

height = inputwd('By what factor should the vertical scale exceed the null hypothesis? (2.0)',2.0)

randn('seed',seeder)

n2=n/2;

f=(0:n2)/(n);

% OK, now construct expected red noise spectrum

for i=1:n2+1

rspec(i)=(1-alpha*alpha)/(1-2*alpha*cos(pi*(i-1)/n2)+alpha*alpha);

end

nn=fix(n2*frac)+1;

factor = sqrt(1-alpha^2);

ymax=rspec(1)*height;

xmax=f(nn)

x=zeros(n,1);

pnum=0;

% loop on realizations

for num = 1:number

pnum=pnum+1;

x(1) = x(n)*alpha + factor*randn;

for i=2:n

x(i)=x(i-1)*alpha + factor*randn;

end

p = spectrum(x,n,0);

if (num == 1 )

psum=p;

else

psum=p+psum;

end

if (pnum == nplot)

pnum=0;

pave=psum/num;

%normalize spectrum

pmmn=squeeze(mean(pave));

pave=pave/pmmn;

%specplot(p,1.)

% find F-statistic

dof = 2* num;

fst = fstat(dof);

spec99=fst*rspec;

plot(f(1:nn),pave(1:nn,1),f(1:nn),rspec(1:nn),f(1:nn),spec99(1:nn),'--','LineWidth',2.0);

%%h = axes('FontSize',14)

axis([0.0,xmax,0.0,ymax])

title(['Red Noise: ','n=',num2str(n),', a=',num2str(alpha), ', realizations=',num2str(num)],'FontSize',18)

xlabel('Frequency (cycles per time step)','FontSize', 16)

ylabel('Power','FontSize', 16)

pause

end

end