function M = triag(N);% Householder triangularization, % N is rectangular matrix, M is square matrix,% M is determined so that the equality  M*M' = N*N' holds% References:% G.W. Steward (1998): Matrix Algorithms, Vol. I: Basic Decomposition.% M.S. Grewal, A.P. Andrews (2001): Kalamn Filtering: Theory and Practise Using MATLAB.% M. Noorgard, N.K. Poulsen, O. Ravn (2000): KALMTOOL.% G.H. Golub, C.F. van Loan (1996): Matrix Computation, 3rd Edition. % Nonlinear Filtering Toolbox version 2.0-rc4% Copyright (c) 1995 - 2007 NFT developement Team,%              Department of Cybernetics,%              University of West Bohemia in Pilsen%initialization[a,b]=size(N);dif = b-a;v = zeros(a,1);u = zeros(1,b);M = zeros(a,a);%Householder triangularizationfor i = a:-1:1,   m(1:dif+i) = N(i,1:dif+i);   ny = norm(m(1:dif+i));   if ny == 0,     m(dif+i) = sqrt(2);   else      m(1:dif+i) = m(1:dif+i)/ny;      if m(dif+i) >= 0,          m(dif+i) = m(dif+i) + 1;        ny = -ny;      else        m(dif+i) = m(dif+i) - 1;      end      m(1:dif+i) = m(1:dif+i)/sqrt(abs(m(dif+i)));   end   v(1:i-1,1) = N(1:i-1,1:dif+i)*m(1,1:dif+i)';   N(1:i-1,1:dif+i) = N(1:i-1,1:dif+i) - N(1:i-1,1:dif+i)*m(1,1:dif+i)'*m(1,1:dif+i);   M(1:i-1,i) = N(1:i-1,dif+i);   M(i,i) = ny;end