fastfixedpoint.m

独立分量分析(Independent Component Analysis

%经典快速固定点算法（基于峭度，梯度下降，超高斯）

clear all;

% 初始化独立源数M和数据长度N
M=3;%4;
sn=3;       %周期数
sl=128;     %每个周期128点
N=sl*sn;    


load EPsignal1r;
% 产生独立信源
signal1=zeros(1,N);         %清零
for i=0:sn-1
    signal1(i*sl+1:(i+1)*sl)=EPsignal1;            %将信号周期延拓，模拟EP信号
end
signal2=((rand(1,N)<0.5)*2-1).*log(rand(1,N));     %模拟超高斯信号 
signal3=randn(1,N);         %产生白高斯信号
%signal4=rand(1,N);         %均匀分布模拟亚高斯信号
%源的均值要求为0，去直流分量
signal1=signal1-mean(signal1);        %去均值
signal1=signal1/max(abs(signal1));    %归一化
signal2=signal2-mean(signal2);
signal3=signal3-mean(signal3);
%signal4=signal4-mean(signal4);
signal=[signal1;signal2;signal3];%;signal4];


% 初始化
v=zeros(M,N);    % 混合后阵矩 v=As
x=zeros(M,N);    % 白化后阵矩 x=Mv=MAs
y=zeros(M,N);    % 输出信号
w=zeros(M,1);    % 解混向量，列
                 % 初值均为零
A=rand(M);       % 混合矩阵，随机产生
%wold=rand(M,1);    % 存放k-1时刻的解混向量，用随机数初始化权向量，放在迭代部分也可以
Wopt=zeros(M,M);    % 权值收敛后的最优权值矩阵
B=zeros(M);         % 由已经求得的解混向量构成的矩阵。为了每次提取出不同的独立源，需要对解混向量进行正交化。
flag=zeros(1,M);       % 判断收敛的标志。收敛，flag=1；不收敛，flag=0。行向量
stepsize=N*ones(1,M);  % 收敛步长初始化，每个源的初始步长均为最大值即总点数


% 信号混合
v=A*signal;


% 白化 
Ev=zeros(M);
eigvalue=zeros(M);    %特征值
eigvector=zeros(M);   %每一列对应一个特征向量
WhiteM=zeros(M);      %白化矩阵
for i=1:N
    Ev=Ev+v(:,i)*v(:,i)';
end
Ev=Ev/N ;             %近似协防差矩阵 E[xx']，平均近似期望
[eigvector,eigvalue]=eig(Ev);
for j=1:M
    eigvalue(j,j)=power(eigvalue(j,j),-0.5);%此时eigvalue为对角阵,对角元为特征值开方
end
WhiteM=eigvector*eigvalue*eigvector';
x=WhiteM*v;


% 权值迭代，产生输出
meanvalue=zeros(M,1);          % 迭代算法中要用到,用来纪录,初始为0
for si=1:M                     %用sn不好==si 表示源
    wold=rand(M,1);            %列，我觉得把他放到循环之内比较好?
    if  si>1
        wold=wold-B*B'*wold;   %投影运算，为满足正交化条件,每个分量之间保持独立
    end
    for i=1:N
        if flag(1,si)==0
            y(si,i)=wold'*x(:,i);
            for j=1:N
                meanvalue=meanvalue+x(:,j)*power(wold'*x(:,j),3);
            end
            meanvalue=meanvalue/N;
            w=meanvalue-3*wold;       %迭代公式实现
            meanvalue=zeros(M,1);     %此步清零很重要           
            if si>1
                w=w-B*B'*w;
            end
            w=w/norm(w);              %归一化,除以范数
            if abs(abs(w'*wold)-1)<1e-6
                flag(1,si)=1;
                stepsize(1,si)=i;     %步长是根据收敛确定的
                Wopt(:,si)=w;
                B(:,si)=w;            %w是矩阵相应的一列，B矩阵用来投影运算
            else
                wold=w; 
            end
        else
            y(si,i)=Wopt(:,si)'*x(:,i);
        end
    end
end 


% 画图
figure(1);                   %源信号
subplot(M,3,1);
plot(signal1);
subplot(M,3,4);
plot(signal2);
subplot(M,3,7);
plot(signal3);
%subplot(M,3,10);
%plot(signal4);

subplot(M,3,2);              %混合信号
plot(v(1,:));
subplot(M,3,5);
plot(v(2,:));
subplot(M,3,8);
plot(v(3,:));
%subplot(M,3,11);
%plot(v(4,:));

subplot(M,3,3);             %输出信号
%plot(x(1,:));
plot(y(1,:));
subplot(M,3,6);
%plot(x(2,:));
plot(y(2,:));
subplot(M,3,9);
%plot(x(3,:));
plot(y(3,:));
%subplot(M,3,12);
%plot(x(4,:));
%plot(y(4,:));



%stepsize=[5	 8	2]


%A=0.11992	0.84077	0.71796
%  0.32012	0.70461	0.13022
%  0.96523	0.14933	0.68027


%Wopt=0.62105	-0.48073	0.61903
%     0.77039	0.51973	   -0.36929
%     -0.1442	0.70625	    0.69313