newicar.m

我最近写的一篇论文的源码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  多维参考信号的ICA算法               by Jan Woo 20060409
%%%%  out1 输出结果；out2分离矩阵W;out3是去除伪差后的结果；
%%%%  说明：

function [out1,out2,out3] = newicar(x,ref);
[M,N]=size(x);
[LM,LN] = size(ref);
disp('--------------------------------------------------------');
disp('This algorithm is written by Jan Woo,April.2006.')
disp('--------------------------------------------------------');
if N~=LN
    disp('错误01：参考信号的长度与观测信号长度不一致！');
else
% introduction
disp('参数说明:');
disp(['观测信号的维数:',int2str(M)]);
disp(['参考信号的维数:',int2str(LM)]);
disp(['输出信号的维数:',int2str(LM)]);
disp('--------------------------------------------------------');
% whiten algorithm
xremean=x-(mean(x'))'*ones(1,N);
RX = cov(xremean',1);
[V,D]=eig(RX);
u=D^(-0.5)*V'*xremean;
% prime algorithm
Wwiener=ref*u'/N;     
WW = zeros(LM,M);

for l = 1:LM    
W(:,1)=(Wwiener(l,:)/(norm(Wwiener(l,:),2)))'; % W is a column; 
yita=norm(Wwiener(l,:),2);
temp=zeros(M,1);

for i=1:1:1000
    for j=1:1:N
        temp=u(:,j)*((W(:,i)'*u(:,j)).^3)+temp;
    end
    W(:,1+i)=temp/N-3*W(:,i);
    W(:,1+i)=W(:,1+i)/norm(W(:,1+i),2);
    if abs(norm(W(:,1+i)'*W(:,i),2)-1)<0.0001
        disp(['第',int2str(l),'个输出经过迭代的次数为：',int2str(i)]);
        break;
    end
    if norm((W(:,i+1)-W(:,1)),2)>yita
        W(:,i+1)=W(:,1)+0.05*randn(M,1);
    end
    W(:,1+i)=W(:,1+i)/norm(W(:,1+i),2);
end %% for i = 1:1000
if i==1000
    disp(['第',int2str(l),'个量经过1000次迭代不收敛。']);
else
    WW(l,:)=W(:,i+1)';
end %% if i == 1000
end %% for l=1:LM
disp('--------------------------------------------------------');
y = WW*u;
Yout = u;
for l = 1:LM
    b = u*y(l,:)'/N;
    Yout = Yout - b*y(l,:);
end
out1 = y;
out2 = WW;
out3 = Yout;

end %% if N~=LN else....