function W=sobi(x,n,p)
%二階盲辨識
%W=sobi(x,n,p)估計分離矩陣
%一般調(diào)用方法W=sobi(x);
%x混合信號
%m傳感器個數(shù)
%N采樣點(diǎn)數(shù)
%n源信號個數(shù)
%p協(xié)方差矩陣個數(shù)
[m,N]=size(x);
if nargin==1
    n=m;
    p=100;
end
if nargin==2
    p=100;
end
xz=x-mean(x')'*ones(1,N);
Rxx=(xz(:,1:N-1)*xz(:,2:N)')/(N-1);
[Ux,Dx,Vx]=svd(Rxx);
Dx=diag(Dx);
if n<m
    Dx=Dx-real((mean(Dx(n+1:m))));
    Q=diag(real(sqrt(1./D(1:n))))*Ux(:,1:n)';
else
    n=max(find(Dx>1e-14));;
    Q=diag(real(sqrt(1./Dx(1:n))))*Ux(:,1:n)';
end
Xb=Q*xz;%prewhitened data
%Estimation of the time delayed covariance matrices
k=1;
pn=p*n;
for u=1:n:pn
    k=k+1;
    Rxp=Xb(:,k:N)*Xb(:,1:N-k+1)'/(N-k+1);
    M(:,u:u+n-1)=norm(Rxp,'fro')*Rxp;
end
%approximate joint daigonalization
eps=1/sqrt(N)/100;
encore=1;
U=eye(n);
while encore
    encore=0;
    for p=1:n-1
        for q=p+1:n
            %givens rotations:
            g=[M(p,p:n:pn)-M(q,q:n:pn);...
                M(p,q:n:pn)+M(q,p:n:pn);...
                i*(M(q,p:n:pn)-M(p,q:n:pn))];
            [Ucp,D]=eig(real(g*g');
            [la,k]=sort(diag(D));
            angles=Ucp(:,k(3));
            angles=sign(angles(1))*angles;
            c=sqrt(0.5+angles(1)/2);
            sr=0.5*(angles(2)-j*angles(3))/c;
            sc=conj(sr);
            encore=encore | asr;
            if asr
                colp=M(:,p:n:pn);
                colq=M(:,q:n:pn);
                M(:,p:n:pn)=c*colp+sr*colq;
                M(:,q:n:pn)=c*colq-sc*colp;
                rowp=M(p,:);
                rowq=M(q,:);
                M(p,:)=c*rowp+sc*rowq;
                M(q,:)=c*rowq-sr(rowp;
                temp=U(:,p);
                U(:,p)=c*U(:,p)+sr*U(:,q);
                U(:,q)=c*U(:,q)-sc*temp;
            end
        end
    end
end
W=U'*Q;
y=W*x;