function [y,m,v] = marco_ica(x,epsilon)

% Perform ICA with precision epsilon on data in x rows and return original 
% signals y and mixing matrix m
% Written by Marco F. Duarte for ECE 631 miniproject
% March 14, 2005

% Set epsilon default
%if isempty(epsilon),
    epsilon = 1e-3;
%end
% Set number of mixes and samples
m = size(x,1);
s = size(x,2);

% Start by whitening the data x into z, with matrix v
% Estimate covariance matrix
cxx = zeros(m);
for i = 1:s,
    cxx = cxx+x(:,i)*x(:,i)';
end
cxx = cxx/s;
% Do eigenvalue decomposition
[e,d]=eig(cxx);
% Whitening matrix
di = diag(diag(d).^(-0.5));
v=e*di*e';
% Whiten data
z=v*x;
% Initialize inversion matrix
wm = [];
% For each separate signal we want to obtain
for i=1:m,
    % Initialize w weight vector
    w = rand(m,1);
    % Orthogonalize wwith previous weight vectors
    for j=1:size(wm,1),
        wc = wm(j,:)*w;
        w = w - wc*wm(j,:)'/(norm(wm(j,:))^2);
    end
    % Normalize w
    w = w/norm(w);
    % Run fixed point algorithm
    dw = 1;
    k = 0;
    while dw > epsilon,
        wold = w;
        % Do fixed point iteration
        % CASE 1: Maximize kurtosis
        t = zeros(m,1);
        for j = 1:s,
            t = t+z(:,j)*((w'*z(:,j))^3);
        end
        % t = E(z*(w'*z)^3)
        t = t/s;
        w = t-3*wold;
        % CASE 2: Maximize negentropy
%         t1 = zeros(m,1);
%         for j = 1:s,
%             t1 = t1+z(:,j)*tanh((w'*z(:,j)));
%         end
%         % t1 = E(z*g(w'*z))
%         t1 = t1/s;
%         t2 = 0;
%         for j = 1:s,
%             t2 = t2+1-(tanh(w'*z(:,j))^2);
%         end
%         % t2 = E(g'(w'*z))
%         t2 = t2/s;
%         w = t1-t2*wold;
        % Orthogonalize with previous weight vectors
        for j=1:size(wm,1),
            wc = wm(j,:)*w;
            w = w - wc*wm(j,:)'/(norm(wm(j,:))^2);
        end
        w = w/norm(w);
        dw = 1-abs(wold'*w);
        k = k + 1;
    end
    disp(['Vector ' num2str(i) ': convergence in ' num2str(k) ' iterations, dw = ' num2str(dw)]);
    % Concatenate new weight vector with previous ones
    wm = [wm; w'];
end
% Separate signals
y = wm*z;
for i=1:size(y,1),
    y(i,:) = y(i,:)/max(abs(y(i,:)));
end
% Get mixing matrix from inverse
m = inv(wm);