function [Y,fc,sig] =SSA_decomposition(x,L,fe,Aff)

% [Y,theta,sig] =SSA_filtering(x,L,fe);
% x signal
% L number of components
% fe sampling frequency
% decomposition in the columns of Y
% fc component frequencies in Hz
% sig variances of the components
% Aff =1 display


x = x(:);
y = x;
N = length(x);

X = zeros(N-L+1,L);
for k=1:N-L+1
    X(k,:) = x(k:k+L-1)';
end

[U,S,V] = svd(X,0);
sig = diag(S);

fc = zeros(L,1);
Y = zeros(length(y), L);
for k=1:L
    Z = fliplr(S(k,k)*U(:,k)*V(:,k)');
    for p=-N+L:L-1,
        y(L-p) = mean(diag(Z,p));
    end
    Y(:,k) = y;
    
    % estimate of component frequency
    [P,fr]=pwelch(y-mean(y),round(N/4),[],4000,fe);
    [~,Ind]=max(P);
    fc(k)=fr(Ind);
    if Aff==1,
       figure(1)
       st = 'blue: signal, red: component #';
       st = cat(2,st,num2str(k));
       plot((1:N)/fe,x,(1:N)/fe,y,'r')
       title(st)
       xlabel('seconds')
       disp(cat(2,'frequence ',num2str(fc(k)),' Hz'));
       pause
    end
end
sig = sig.*sig;
if Aff==1,
figure(2)
[fcs,Ind]=sort(fc);
plot(fcs,sig(Ind),'--rs','MarkerFaceColor','r','MarkerSize',8)
xlabel('component frequencies')
ylabel('component variances')
end





    