% This program generates a distribution approximating a random walk, but
% without the humungus flicker component produced by integrating a
% Gaussian distribution. To make it more realistic, we start with an
% exponential distribution, then integrate over a window of 1000
% samples.
%
% The time series at each lag should be almost identical for lags up to
% the fileter length.
%
clear all;
max = 1000000;
lng = 1000;
phase = 1;
%
% Generate and filter the distribution.
%
expon = random('Exponential', phase, max, 1);
lng = 1000;
x1 = linspace(1, max - 1, max);
a = linspace(0, 0, lng);
a(1) = 1;
b = linspace(1, 1, lng) / lng;
%for i = 1:lng
%    b(i) = i .^ -0.1;
%end
%b = b / sum(b);
y = filter(b, a, expon);
%
% Plot the signal for each lag.
%
y1 = y; i = 1; d = 1;
while length(y1) >= 100    
	z1(i) = var(y1);
	m1(i) = d;
    clf reset; h = newplot;
    set(h, 'FontSize', 12);
    set(h, 'Position', [.12 .15 .85 .8]);
    plot(x1 / 1000, y1, '-k')
    axis([0 1000 .8 1.2]);
    lag = strcat('Time (ks) Interval=', num2str(d));
    xlabel(lag);
    ylabel('Delay (ms)');
    lag = strcat('test_', num2str(d));
    print('-dtiff', '-r600', lag)
    x1 = (x1(2:2:length(x1)));
    y1 = (y1(2:2:length(y1)) + y1(1:2:length(y1) - 1)) / 2;
    i = i + 1; d = 2 * d;
end
%
% Plot the variance-time graph.
%
e = m1 / m1(1);
e = 1 ./ e;
loglog(m1, z1, '-k', m1, e * z1(1))
%axis([1 1e5 1e-4 100]);
xlabel('Time Interval (s)');
ylabel('Variance (s^2)');
print -dtiff -r600 test_var
%
% Plot the periodogram.
%
f = fft(expon);
r = f .* conj(f);
loglog(r(length(r) / 2:length(r)), '.')
%axis([1 1e5 1e-4 100]);
xlabel('Frequency (Hz)');
ylabel('Power');
print -dtiff -r600 test_per
%
% fft
%
f = fft(y);
r = f .* conj(f);
loglog(r(length(r) / 2:length(r)), '-k')
%axis([1 1e5 1e-4 100]);
xlabel('Frequency (Hz)');
ylabel('Power');
print -dtiff -r600 test_per