function [rayamp]=enkelsumma(f,sigma2_,N,t,plotid)
% ENKELSUMMA simulerar ett N st realiseringar av en st. proc. där
%
% X(t) = A0 + sum Ak cos(2 pi fk t + phik)
%
% A0 är normalfördelad med varians sigma2(0)
% Ak är Rayleigh-fördelad med parameter sigma2(k)
% 
% Invariabler,
% f       1-ggr-K vektor med önskade frekvenser i växande ordning
% sigma2_ 1-ggr-K vektor med önskade energier
% N       antal realiserningar
% t       tid
% 
% Utvariabler,
% rayamp  Normalprocessen X(t)
%
% Exempel på användning:
% f = [5 10];
% sigma2_ = [2 2];
% N = 1000;
% dt = 1/(2*max(f)+1);
% t = 0:dt:20;
% [rayamp,normamp]=enkelsumma(f,sigma2_,N,t);
%
% För att använda spekgui
% data.x = rayamp(:,5);
% data.dt = dt;

% Anders Malmberg 20040305


if nargin < 5
     plotid=1;
end

if sigma2_(1)==0
sigma2 = sigma2_(1:end);
else
sigma2 = sigma2_(:);
end


sigma2 = repmat(sigma2,1,N); % K-by-N matrix
%varR = sigma2*(4-pi)/2
R = raylrnd(sqrt(sigma2)); % K-by-N matrix
if sigma2_(1) == 0
   A0 = normrnd(0,sqrt(sigma2_(1)),length(f),1);
   R = [A0 R];
end
%var(R)
phi = 2*pi*rand(length(f),N);

% time interval
t = t(:);
T = repmat(t,1,N);
% signal
X = zeros(size(T));
for i = 1:length(f)
     r = R(i,:);
     r = repmat(r,size(T,1),1);
     p = phi(i,:);
     p = repmat(p,size(T,1),1);
     X = X+r.*cos(2*pi*f(i)*T + p);
end

if plotid
figure(1)
clf
subplot(221)
% spektrum
if f(1) == 0
   stem(f,[sigma2_(1) sigma2_(2:end)./2],'filled')
else
   stem(f,[sigma2_./2],'filled');
end
title('Spektraltäthet');
xlabel('f')
ylabel('R(f)')
axis([min(f)-.5 max(f)+.5 0 max(sigma2_)+1])
grid
subplot(222)
tau = 0:.01:(1/min(f)*4);
r = zeros(size(tau));
for i = 1:length(f)
r = r + sigma2_(i)*cos(2*pi*f(i)*tau);
end
plot(tau,r)
grid on
title('Kovariansfunktion')
     xlabel('\tau')
     ylabel('r(\tau)')
subplot(223)
itend = (1/min(f)*5);
it = t<itend;
plot(t(it),X(it,1),t(it),X(it,2)+5,'r',t(it),X(it,3)-5,'black')
axis([0 t(sum(it)) -15 15])
title('Exempel på realiseringar')
xlabel('t')
ylabel('X(t)')
legend('X_1(t)','X_2(t)+5','X_3(t)-5')
grid on
subplot(224)
it = 45;
t_0 = (t(2)-t(1))*45;
normplot(X(45,:))
tit = sprintf('Normalplot för X(t) vid t=%0.5g i %d realiseringar',t_0,N);
title(tit)
end

% Utdata
rayamp = X;