clear all
clc

%觀察信號矩陣 X 為 K*N 矩陣
%源信號矩陣 S 為 K*N 矩陣
%混合矩陣 A 為 K*K 矩陣

K=4;
N=157;

%生成 K 個源信號, 存儲到 S
t=1:N;
figure(1);

S(1,:)=rand(1,N);
subplot(2,2,1);
plot(t,S(1,:));
title('隨機信號');
hold on;

S(2,:)=square(2*pi*t/64);  
subplot(2,2,2);
plot(t,S(2,:));
title('方波信號');
hold on;

S(3,:)=sin(2*pi*t/32);  
subplot(2,2,3);
plot(t,S(3,:));
title('正弦信號');
hold on;

S(4,:)=cos(2*pi*t/16); 
subplot(2,2,4);
plot(t,S(4,:));
title('余弦信號');
hold on;

%生成 K 個混合信號,存儲到 X, A 中存儲混合矩陣
A=randn(K);
X=A*S;
figure(2);
plot(t,X);
title('觀察信號');

% X 數據中心化
m=mean(X,2);
for i=1:N
    X(:,i)=X(:,i)-m
end

% X 數據白化
covMat=cov(X');
[E,D]=eig(covMat);
V=E*D^(-0.5)*E';
X=V*X;
figure(3);
plot(t,X);
title('白化后的信號');

%FastIca算法
%求完后 W 是 混合矩陣 A 的近似
W=rand(K);
for p=1:K
    W(:,p)=W(:,p)/norm(W(:,p));
    e=1
    while e<800
        pred=W(:,p);
        
        %定點法
        W(:,p)=1/N*X*((W(:,p)'*X).^3)'-3*W(:,p);

        %采用收縮策略，類似于正交化的過程
        sum=zeros(K,1);
        for i=1:p-1
            sum=sum+W(:,p)'*W(:,i)*W(:,i);
        end
        W(:,p)=W(:,p)-sum;
        W(:,p)=W(:,p)/norm(W(:,p));
        
        %收斂條件
        e=e+1
    end
end

%推算源信號
out=W'*X

figure(4);
for k=1:K  
    subplot(2,2,k);
    plot(t,out(k,:));
    hold on;
end;