comfac.m

基于卷积信号的MIMO系统盲信号估计

CC = Ct*Cg;
AA = X*pinv(ppp(BB,CC)).';

if SmallMode == 1
  A=AA;
  B=BB;
  C=CC;
elseif SmallMode == 2 
  A=BB;
  B=AA;
  C=CC;
elseif SmallMode == 3
  A=BB;
  B=CC;
  C=AA;
end


if ~DontShowOutput
  fit = sum(sum(abs(X - AA*ppp(BB,CC).').^2));
  disp([' DTLD fitted raw data with a sum-squared error of ',num2str(fit)])
end


function [A,B,C,G,fit,it]=tucker3(X,DimX,W,maxit);

% [A,B,C,G]=tucker3(X,R,W,maxit);
% This is an SVD-based algorithm for finding the
% parameters of the three-way Tucker3 model when
% the array as well as parameters are complex
% 
% Copyright 1998
% Rasmus Bro & Claus A. Andersson
% KVL, Denmark, [e-mail]

UseNIPALS = 0; % Use NIPALS (svdf.m) instead of SVD. It's cheaper in terms of flops for large arrays. For small there's no big difference
if UseNIPALS == 1 & any(imag(X(:)))
   disp(' Apparently the loadings in this NIPALS are a little oblique. Check that (decrease convergence criterion)')
   error(' NIPALS has not yet been changed to handle complex numbers. Use SVD (set UseNIPALS to 0)')
end

DontShowOutput = 0;

if nargin<4
  maxit=100;
end

%Initialising counters and others
SSX=sum(sum(abs(X).^2));
it=0;
Oldfit=1e100;
Diff=1e100;
I=DimX(1);J=DimX(2);K=DimX(3);
B=orth(rand(J,W(2)));
C=orth(rand(K,W(3)));

while Diff>1e-6&it<maxit

  it=it+1;

  %Updating A
    TA1=C'*reshape(X,I*J,K).';
    TA2=B'*reshape(TA1,W(3)*I,J).';
    TA3=reshape(TA2,W(2)*W(3),I).';
    if UseNIPALS
      [TA4 TA5 TA6]=svdf(TA3,W(:,1));
    else
      [TA4 TA5 TA6]=svd(TA3,0);
    end
    A=TA4(:,1:W(1));

  %Updating C
    TC1=reshape(A'*X,W(1)*J,K).';
    TC2=B'*reshape(TC1,W(1)*K,J).';
    TC3=reshape(TC2,K*W(2),W(1)).';
    TC3=reshape(TC3,W(1)*W(2),K).';
    if UseNIPALS
      [TC4 TC5 TC6]=svdf(TC3,W(3));
    else
      [TC4 TC5 TC6]=svd(TC3,0);
    end
    C=TC4(:,1:W(3));

  %Updating B
    TB1=reshape(C'*TC1,W(1)*W(3),J).';
    if UseNIPALS
       [TB2 TB3 TB4]=svdf(TB1,W(2));
    else
       [TB2 TB3 TB4]=svd(TB1,0);
    end
    B=TB2(:,1:W(2));

  %Calculate core & fit
    G1=reshape(A'*X,W(1)*J,K).';
    G2=reshape(C'*G1,W(1)*W(3),J).';
    G=reshape(B'*G2,W(2)*W(3),W(1)).';
    fit=sum(sum(abs(G.^2)));
    fit=SSX-fit;

    Diff=abs(Oldfit-fit);
    Oldfit=fit;

end

function [u,s,v] = svdf(X,F);

% Rand-reduced SVD based on NIPALS (actually 
% the power-method the way it's implemented here)

maxit = 30; % Use 30
crit = 1e-4; % Use 1e-4
[I,J] = size(X);
u = zeros(I,F);
v = zeros(J,F);

if J > I
   x = X*X';
else
   x = X'*X;
end


for f = 1:F
   p = sum(x)';
   converged=0;
   it = 0;
   while ~converged
      it = it +1;
      pold = p;
      p = x*p;
      p = p/norm(p);
      if norm(p-pold)/norm(pold)<crit | it>maxit
         converged = 1;
      end
   end
  
   if J > I
     u(:,f) = p;
     v(:,f) = X'*p;
     s(f) = norm(v(:,f));     
     v(:,f) = v(:,f)/s(f);
   else
     v(:,f) = p;
     u(:,f) = X*p;
     s(f) = norm(u(:,f));
     u(:,f) = u(:,f)/norm(u(:,f));
   end
   x = x - s(f)^2*p*p';
end


function [X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12,X13,X14,X15,X16,X17]=nshape(X,DimX,f);

% $ Version 1.03 $ Date 18. July 1999 $ Not compiled $
% $ Version 1.031 $ Date 18. July 1999 $ Error in help figure and now outputs new DimX $ Not compiled $
%
% Copyright, 1998 - 
% This M-file and the code in it belongs to the holder of the
% copyrights and is made public under the following constraints:
% It must not be changed or modified and code cannot be added.
% The file must be regarded as read-only. Furthermore, the
% code can not be made part of anything but the 'N-way Toolbox'.
% In case of doubt, contact the holder of the copyrights.
%
% Rasmus Bro
% Chemometrics Group, Food Technology
% Department of Food and Dairy Science
% Royal Veterinary and Agricultutal University
% Rolighedsvej 30, DK-1958 Frederiksberg, Denmark
% Phone  [phone]
% Fax    [phone]
% E-mail [e-mail]
%
%
% [Xf,DimXf] = nshape(X,DimX,f);
%
% Refolds an N-way array so that Xf is X with index
% f as row-index, and the remaining in succesive order. For an 
% I x J x K x L four-way array this means X1 is I x JKL, X2 is
% J x 蘇L, X3 is K x IJL, and X4 is L x IJK
%
%
%    K  _______             
%      /      /|           1      J     2稪    J稫
%     /______/ |         1  _____________________
%    |      |  |           |      |      |      |
%    |      | /    -->     |      |      |      |        f = (Mode) 1 (same as original array)
% I  |______|/          I  |______|______|______|
%           J
%
%                          1      I     2稩    K稩
%                        1  _____________________
%                          |      |      |      |
%                  -->     |      |      |      |        f = (Mode) 2
%                        J |______|______|______|
%
%  
%                          1      I     2稩    I稪
%                        1  _____________________
%                          |      |      |      |
%                  -->     |      |      |      |        f = (Mode) 3
%                        K |______|______|______|
%
%
% If the last input is not given all rearrangements are given.
% For a fourway array this would read
% [X1,X2,X3,X4]=nshape(X,DimX);
%
%	Copyright
%	Rasmus Bro & Claus A. Andersson 1995
%	Denmark
%	E-mail [e-mail]

ord=chkpfdim(X,DimX,NaN);
elemen=prod(DimX);
if nargin==2
  do_it=ones(1,ord);
else
  do_it=zeros(1,ord);
  do_it(f)=1;
end

if do_it(1)==1
  X1=X;
end


% _____Make X2_____

if do_it(2)==1
  X2=X(:,1:DimX(2)).';
  for R2=DimX(2)+1:DimX(2):elemen/DimX(1)
    X2=[X2 X(:,R2:R2+DimX(2)-1).'];
  end
end

if ord>3

% _____Make X3 - Xord-1_____

for comp=3:ord-1
  if do_it(comp)==1
    xx=[];  % Denne kan opbygges til flere
    for R2=1:prod(DimX(2:comp-1)):prod(DimX(2:comp))
      x=[];
      for R3=R2:prod(DimX(2:comp)):prod(DimX(2:ord))
        x=[x reshape(X(:,R3:R3+prod(DimX(2:comp-1))-1),1,prod(DimX(1:comp-1)))];
      end % for
      xx=[xx;x];
    end
    eval(['X',num2str(comp),'=xx;']);
  end % for comp
end
end % if ord>3


% _____Make Xord_____

if do_it(ord)==1
  xx=[];
  for R3=1:elemen/(DimX(1)*DimX(ord)):elemen/DimX(1)
    xx=[xx; reshape(X(:,R3:R3+elemen/(DimX(1)*DimX(ord))-1),1,elemen/DimX(ord))];
  end % for
  eval(['X',num2str(ord),'=xx;']);
end

if nargin==3
   eval(['X1=X',num2str(f),';']);
   X2 = [DimX(f) DimX([1:f-1 f+1:ord])];
end

function ord=chkpfdim(X,DimX,show);

% show == NaN => no text
%
% $ Version 1.02 $ Date 28. July 1998 $ Not compiled $

if nargin<3
   show=0;
end


% Find the order, i.e., number of ways
ord=length(DimX);

% Check if DimX corresponds to size of X
if DimX(1)~=size(X,1)|prod(DimX)/DimX(1)~=size(X,2)
  disp(' ')
  disp(' Size of array does not correspond to dimensions given in DimX')
  error(['disp('' The matrix input must be of size ',num2str(DimX(1)),' x ',num2str(prod(DimX)/DimX(1)),' if DimX is correctly given'')'])
end

if ~isnan(show)
  txt=[];
  for i=1:ord-1
    txt=[txt num2str(DimX(i)) ' x '];
  end
  txt=[txt num2str(DimX(ord))];
  disp([' The array is a ',num2str(ord),'-way array with'])
  disp([' dimensions: ' txt])
end

function AB=ppp(A,B);

% $ Version 1.02 $ Date 28. July 1998 $ Not compiled $
%
% Copyright, 1998 - 
% This M-file and the code in it belongs to the holder of the
% copyrights and is made public under the following constraints:
% It must not be changed or modified and code cannot be added.
% The file must be regarded as read-only. Furthermore, the
% code can not be made part of anything but the 'N-way Toolbox'.
% In case of doubt, contact the holder of the copyrights.
%
% Rasmus Bro
% Chemometrics Group, Food Technology
% Department of Food and Dairy Science
% Royal Veterinary and Agricultutal University
% Rolighedsvej 30, DK-1958 Frederiksberg, Denmark
% Phone  [phone]
% Fax    [phone]
% E-mail [e-mail]
%
% The parallel proportional profiles product - triple-P product
% For two matrices with similar column dimension the triple-P product
% is ppp(A,B) = [kron(B(:,1),A(:,1) .... kron(B(:,F),A(:,F)]
% 
% AB = ppp(A,B);
%
% Copyright 1998
% Rasmus Bro
% KVL,DK
% [e-mail]

[I,F]=size(A);
[J,F1]=size(B);

if F~=F1
   error(' Error in ppp.m - The matrices must have the same number of columns')
end

AB=zeros(I*J,F);
for f=1:F
   ab=A(:,f)*B(:,f).';
   AB(:,f)=ab(:);
end

function [NewA,NewB,NewC,DeltaMin]=linesrch(X,DimX,A,B,C,Ao,Bo,Co,Delta);

dbg=0;

if nargin<5
  Delta=5;
else
  Delta=max(2,Delta);
end

dA=A-Ao;
dB=B-Bo;
dC=C-Co;
Fit1=sum(sum(abs(X-A*ppp(B,C).').^2));
regx=[1 0 0 Fit1];
Fit2=sum(sum(abs(X-(A+Delta*dA)*ppp((B+Delta*dB),(C+Delta*dC)).').^2));
regx=[regx;1 Delta Delta.^2 Fit2];

while Fit2>Fit1
  if dbg
    disp('while Fit2>Fit1')
  end
  Delta=Delta*.6;
  Fit2=sum(sum(abs(X-(A+Delta*dA)*ppp((B+Delta*dB),(C+Delta*dC)).').^2));
  regx=[regx;1 Delta Delta.^2 Fit2];
end

Fit3=sum(sum(abs(X-(A+2*Delta*dA)*ppp((B+2*Delta*dB),(C+2*Delta*dC)).').^2));
regx=[regx;1 2*Delta (2*Delta).^2 Fit3];

while Fit3<Fit2
  if dbg
    disp('while Fit3<Fit2')
  end
  Delta=1.8*Delta;
  Fit2=Fit3;
  Fit3=sum(sum(abs(X-(A+2*Delta*dA)*ppp((B+2*Delta*dB),(C+2*Delta*dC)).').^2));
  regx=[regx;1 2*Delta (2*Delta).^2 Fit2];
end

% Add one point between the two smallest fits
[a,b]=sort(regx(:,4));
regx=regx(b,:);
Delta4=(regx(1,2)+regx(2,2))/2;
Fit4=sum(sum(abs(X-(A+Delta4*dA)*ppp((B+Delta4*dB),(C+Delta4*dC)).').^2));
regx=[regx;1 Delta4 Delta4.^2 Fit4];

%reg=pinv([1 0 0;1 Delta Delta^2;1 2*Delta (2*Delta)^2])*[Fit1;Fit2;Fit3]
reg=pinv(regx(:,1:3))*regx(:,4);
%DeltaMin=2*reg(3);

DeltaMin=-reg(2)/(2*reg(3));

%a*x2 + bx + c = fit
%2ax + b = 0
%x=-b/2a

NewA=A+DeltaMin*dA;
NewB=B+DeltaMin*dB;
NewC=C+DeltaMin*dC;
Fit=sum(sum(abs(X-NewA*ppp(NewB,NewC).').^2));

if dbg
  regx
  plot(regx(:,2),regx(:,4),'o'),
  hold on
  x=linspace(0,max(regx(:,2))*1.2);
  plot(x',[ones(100,1) x' x'.^2]*reg),
  hold off
  drawnow
  [DeltaMin Fit],pause
end

[minfit,number]=min(regx(:,4));
if Fit>minfit
  DeltaMin=regx(number,2);
  NewA=A+DeltaMin*dA;
  NewB=B+DeltaMin*dB;
  NewC=C+DeltaMin*dC;
end