%% Script file to produce Example 1, ICA

% Create the data:  A noisy circle in 3-d
numpts=400;
t=linspace(0,3*pi,numpts);
S=[cos(t'), sin(t'), randn(numpts,1)]; %Separated (Clean) Signals
A =[ -0.0964 -0.1680 1.6777
      0.4458 0.1795 1.9969
     -0.2958 0.4211 0.6970]; %Mixing Matrix
X=S*A;       % Mixed Data
dX=diff(X);  % Differenced data (derivative)

%% Method 1:  GSVD:

[U,V,B,C,S] = gsvd(X,dX,0); %Clean Signal in U

%% Method 2:  Two SVDs
%The double SVD code, equivalent to the GSVD:
[Ux,Sx,Vx]=svd(X,0);
Y=X*Vx*diag(1./(diag(Sx))); %Stopping here is basic PCA
dY=diff(Y);
[Uy,Sy,Vy]=svd(dY,0);
S2=Y*Vy; %S2 is also the clean signal


%% Plot the results.

% Plotting routines below:
figure(1)  % Plot the original data in R^3
plot3(X(:,1),X(:,2),X(:,3));
view(-133,26)
axis tight

figure(2)
for j=1:3
  subplot(3,1,j)
  plot(U(:,j)); %Clean Signals from GSVD
end

figure(3)
for j=1:3
  subplot(3,1,j);
  plot(S2(:,j)); %Clean Signals from double SVD
end

figure(4)
for j=1:3
  subplot(3,1,j)
  plot(Y(:,j)); %Results of PCA (or KL)
end