jader.m

这是一个关于盲信号处理的联合对角化的程序希望对大家有用

%%%%%%%%%%%%%%%%%%%%%%%%联合对角化实信号盲分离算法%%%%%%%%%%%%%%%%%%%%%%%%%%
function [W]=jadeR(X)
verbose=0;%信息显示控制变量
%%%%%%%%%%%%%%确定原信号个数%%%%%%%%%%%
[n,T]=size(X);
if nargin==1
    m=n;%原信号数据默认为观测信号个数
end
if m>n
    fprintf('jade->原信号个数不能大于传感器数目！\n')
return
end
if verbose 
    fprintf('jade->Looking for %d sources\n',m);
end
%%%%%%%%%%%%%%%%%%%%%%%%%预处理%%%%%%%%%%%%%%%%%%%%%%%%%%%
if verbose
    fprintf('jdae->去均值\n');
end
X=X-mean(X')'*ones(1,T);
%%%%%%%%%%%%%%%%%%%%%%%%%白化确定抽取信号个数%%%%%%%%%%%%%%%%%%%%
if verbose
    fprintf('jade->Whitening the data\n');
end
[U,D]=eig((X*X')/T);
[puiss,k]=sort(diag(D));
rangeW=n-m+1:n;%indices to the m most significant dirctions
scales=sqrt(puiss(rangeW));%scales
W=diag(1./scales)*U(1:n,k(rangeW))';%白化
iW=U(1:n,k(rangeW))*diag(scales);%its pseudo-inverse
X=W*X;
%%%%%%%%%%%%%%%%%%%%%%%%%%%累积量矩阵估计%%%%%%%%%%%%%%%
if verbose
    fprintf('jade->估计累积量矩阵\n');
end
dimsymm=(m*(m+1))/2;%是对称矩阵的维数
nbcm=dimsymm;%累积量矩阵个数
CM=zeros(m,m*nbcm);%存储累积量矩阵
R=eye(m);
Qij=zeros(m);%累积量矩阵的临时变量
Xim=zeros(1,m);%临时
Xjm=zeros(1,m);%临时
scale=ones(m,1)/T;
Range=1:m;%累积量矩阵存储的指数
for im=1:m
    Xim=X(im,:);
    Qij=((scale*(Xim.*Xim)).*X)*X'-R-2*R(:,im)*R(:,im)';
    CM(:,Range)=Qij;
    Range=Range+m;
    for jm=1:im-1
        Xjm=X(jm,:);
        Qij=((scale*(Xim.*Xjm)).*X)*X'-R(:,im)*R(:,jm)'-R(:,jm)*R(:,im)';
        CM(:,Range)=sqrt(2)*Qij;
        Range=Range+m;
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%累积量矩阵的联合对角化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%
if 1 %%%%%%%%%%%%%%%%%%%%对角化信号累积量初始化，节约时间
    if verbose
        fprintf('jade->对角初始化\n');
    end
    [V,D]=eig(CM(:,1:m));
    for u=1:m:m*nbcm
        CM(:,u:u+m-1)=CM(:,u:u+m-1)*V;
    end
    CM=V'*CM;
else V=eye(m);
end


seuil=1/sqrt(T)/100;%统计门限
encore=1;
sweep=1;
updates=0;
g=zeros(2,nbcm);
gg=zeros(2,2);
c=0;
s=0;
ton=0;
toff=0;
theta=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%联合对角化过程%%%%%%%%%%%%%%%%%%%%%%%%%%%
if verbose
    fprintf('jade->利用联合对角优化对比函数\n');
end
while encore
    encore=0;
    if verbose
        fprintf('jade->Sweep#%d\n',sweep);
    end
    sweep=sweep+1;
    for p=1:m-1
        for q=p+1:m
            Ip=p:m:m*nbcm;
            Iq=q:m:m*nbcm;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算Givens变换角度%%%%%%%%%%%%%%%%
            g=[CM(p,Ip)-CM(q,Iq);CM(p,Iq)+CM(q,Ip)];
            gg=g*g';
            ton=gg(1,1)-gg(2,2);
            toff=gg(1,2)+gg(2,1);
            theta=0.5*atan2(toff,ton+sqrt(ton*ton+toff*toff));
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Givens变换%%%%%%%%%%%%%%%%%%%%%
            if abs(theta)>seuil
                encore=1;
                updates=updates+1;
                c=cos(theta);
                s=sin(theta);
                G=[c -s;s c];
                pair=[p;q];
                V(:,pair)=V(:,pair)*G;
                CM(pair,:)=G'*CM(pair,:);
                CM(:,Ip )=[c*CM(:,Ip)+s*CM(:,Iq)-s*CM(:,Ip)+c*CM(:,Iq)];
            end
        end
    end
end
if verbose
    fprintf('jade->Total of %d Givens rotations\n',updates);
end


W=V'*V;
%%%%%%%%%%%%%%%%%%%%%%%%%%行置换求第一分量，这个信号被规格化为单位协方差%%%%%%%%%%%%%%%%%%%%%%%%%%
if verbose
    fprintf('jdae->分量排序\n',updates);
end
A=iW*V;
[vars,keys]=sort(sum(A.*A));
W=W(keys,:);
W=W(m:-1:1,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%固定符号，使W的第一行非负%%%%%%%%%%%%%%%%%%%%%%%%%%
if verbose
    fprintf('shubbs->固定符号\n',uodates);
end
b=W(:,1);
signs=sign(sign(b)+0.1);
W=diag(signs)*W;