%Figures produced by Whittle likelihood smoothing will appear differently 
%from the analogous R versions (using simple Gaussian smoothers on the 
%periogograms); often much more smoothed
lfp = load('lfp-ryan.dat');

time = 1:1000;
lfp1 = lfp(1:1000);
f = (0:500)./1000;

knots = [1, 200:200:1000];
splinefit1 = spline(knots, lfp1(knots), time);

lfp1detrend = lfp1-splinefit1';

I1detrend = periodogram(lfp1detrend)';

maxt1 = size(I1detrend, 2);
smoothI1 = zeros(1, maxt1);
bw = 1;

for i = 1:maxt1
    I1ker = normpdf((i-(1:maxt1))/bw);
    smoothI1(i) = sum(I1ker.*I1detrend)/sum(I1ker);
end

sm_I1d_50 = smoothI1(1:50)';

sm_Ivals = zeros(1000, 50);
upperband = zeros(1, 50);
lowerband = zeros(1, 50);

bw = 1;

for i = 1:1000
    rvals = gamrnd(1, 1, 50, 1);
    Ivals = sm_I1d_50.*rvals;

    smoothrep = zeros(1, 50);
    for j = 1:maxt1
        I1kerrep = normpdf((j-(1:50))/bw)';
        sm_Ivals(i, j) = sum(I1kerrep.*Ivals)/sum(I1kerrep);
    end
end

for i = 1:50
	sort_vals = sort(sm_Ivals(:,i));
	upperband(i) = sort_vals(975);
	lowerband(i) = sort_vals(25);
end


plot(f(1:50), sm_I1d_50, '-k', 'LineWidth', 2)
hold on;
plot(f(1:50), lowerband, '-b', 'LineWidth', 2)
plot(f(1:50), upperband, '-b', 'LineWidth', 2)

set(gca, 'Box', 'off', 'FontSize', 20, ...
         'XLim', [-0.001, 0.051], 'YLim', [-100, 6500], ...
         'XTick', 0:0.01:0.05, 'YTick', 0:1500:6000, ...
         'TickDir', 'out')
xlabel('Frequency', 'FontSize', 20)

set(gcf, 'Position', [200, 100, 1200, 800])