www.pudn.com > newicar.rar > 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....