function kap8_3c()
% find warmup time with Welch method
k = 10;
n = 1000;
tEnd = n;
t = (1:tEnd);
X = zeros(k, n);

load_system('singleserver4B');

% collect data from k runs
for I=1:k
  NQ = runSingleserver4B(2*I, tEnd);
  fprintf('  Runde %3d\n', I)
  tIn = NQ.time;
  nQIn = NQ.signals.values;
  nQ = interp1(tIn, nQIn, t, 'previous');
  X(I,:) = nQ;
end

% average over runs
Xbar = mean(X);

% smoothing
w = [100, 250, 500];
Xsm1 = movmean2(Xbar, w(1)); 
Xsm2 = movmean2(Xbar, w(2)); 
Xsm3 = movmean2(Xbar, w(3)); 

% Plots
set(gcf(), 'Position',[1 1 500 900], 'Color',[1 1 1]);
subplot(4,1,1)
plot(t, X(1,:), t, X(2,:), t, X(3,:))
title('Beispielkurven')
xlim([0,tEnd])

subplot(4,1,2)
plot(t, Xbar)
title('Mittelwert über Läufe')
xlim([0,tEnd])

subplot(4,1,3)
plot(t, Xsm1, t, Xsm2, t, Xsm3)
title('geglättete Mittelwerte')
xlim([0,tEnd])
legend(sprintf('w=%d',w(1)), sprintf('w=%d',w(2)), ...
       sprintf('w=%d',w(3)), 'Location', 'Best')

subplot(4,1,4)
plot(t, Xsm2)
title('ausgewählte Kurve')
xlim([0,tEnd])

F = getframe(gcf(), [40, 75, 435, 780]);
imwrite(F.cdata, 'bild242.png');

%----------------------------------------------------------------
function y = movmean2(x, b)
% substitute for Matlabs movmean with shrinking window width
% at the beginning
k = 2*b+1;
yr = movmean(x, k);
y = [zeros(1,b), yr((b+1):end)];

for I=1:b
  y(I) = sum(x(1:(2*I-1)))/(2*I-1);
end
