mypca.m

该代码为盲源分离中的相对牛顿法

% MYPCA - principal component analysis
%
%          [Y,B,R0]=mypca(X0,N) 
%
%Input:
%          X0 - input signals in rows
%          N  - number of princip. components
%
%Output:
%          Y  - normalized principal components
%          B  - transformation matrix: Y=B*(X0-mean(X0))
%          R0 - cov(X0);

function    [Y,B,R0]=mypca(X0,N),

  [Nsn,TotT]=size(X0);
  
  if nargin<2 N=Nsn;end

  X=X0-mean(X0,2)*ones(1,TotT);
  R0=cov(X');

  [S0,lam0]=eig(0.5*(R0+R0'));lam0=diag(lam0);
  [lm0,I0]=sort(lam0);%lm(Nsn:-1:1)'
  II0=I0(Nsn:-1:Nsn-N+1);
  lam1=lam0(II0);
  S1=S0(:,II0);
  B=(S1*diag(sqrt(1 ./(lam1+1e-20))))';
  Y=B*X;  % presphearing