Code covered by the BSD License  

Highlights from
Blind Modal Identification with Hankel Matrices (BMIDHM) toolbox

image thumbnail

Blind Modal Identification with Hankel Matrices (BMIDHM) toolbox

by

 

03 Jun 2013 (Updated )

This toolbox is for performing modal identification from ambient random data.

[W,R,V,L,ord,cost]=hmsobi(x,nhank,dec,ord,nbrc)
function [W,R,V,L,ord,cost]=hmsobi(x,nhank,dec,ord,nbrc) 
% 
% Type: [W,R,V,L,ord]=hmsobi(x,nhank,dec,ord,nbrc);
% 
% Inputs:
%
% x        := nt x m complex-valued input data matrix (no whitening needed) (complex)
% nhank    := 1 x 1 the number of Hankel matrices for JAD (default=10) 
% dec      := 1 x 1 decimation order between Hankel matrices (default=1)
% ord      := 1 x 1 order (number of sources)
%             ord = 0 => choose order interactivly (default=0)
% nbrc     := 1 x 1 number of block rows and columns in Hankel matrices (default=1)
%
%
% Outputs:
%
% W   := m x m Inverse of right diagonalizer. Estimate of the inverse of mixing matrix in first block row. (complex)
% R   := m x m Right diagonalizer. Estimate of the mixing matrix in first block column. (complex) 
% V   := m x m Inverse of left diagonalizer (complex)
% L   := m x m Left diagonalizer (complex) 
% ord := 1 x 1 order (number of sources) used
%
% Function for SOBI algorithm applied to complex-valued data. The JAD routine employed is
% wsf_AB. Find decomp: Hxx(tau) = L*Rqq(tau)*R', where Rqq(tau) is diagonal, Hxx(tau) is the Hankel matrix of Rxx(tau).
%
% This functon calls: cstdcov, wsf_AB

% Scot McNeill, January 2013
%
%
error(nargchk(1,5,nargin));
if nargin<5 | isempty(nbrc)
 nbrc=1;
end
if nargin<4 | isempty(ord)
 ord=0;
end
if nargin<3 | isempty(dec)
 dec=1;
end 
if nargin<2 | isempty(nhank)
 nhank=10;
end
%
[nt,m]=size(x);
x=x.';
%
% Compute time-delayed complex-valued correlation matrices
%
nth=dec*nhank+2*nbrc-2;
Rx=zeros(m,nth*m);
for i=1:nth
 ii=(i-1)*m+1:i*m;
 Rx(:,ii)=cstdcov(x,i);
end
%
% Build Hankel Matrices
%
icol=1:nbrc*m;
Hx=zeros(nbrc*m,nbrc*m,nhank);
for i=1:nhank
 for k=1:nbrc
  irow=[(k-1)*m+1:k*m];
  i1=(i-1)*dec*m+(k-1)*m+[1:nbrc*m];
  Hx(irow,icol,i)=Rx(:,i1);
 end
end
%
% Order selection
%
if ord == 0
 svl=svd(squeeze(Hx(:,:,1)));
 figure(81); 
 semilogy(svl,'*');
 xlabel(' Number');ylabel('SV Magnitude');
 title(' Hankel Singular Values (HSV)');
 disp([' ']);
 disp(['The HSV plot (See Figure) allows you to select the # of modes'])
 disp(['by choosing the desired number of singular values.']);
 disp(['# of singular values should be equal to the number of modes.']);
 disp([' ']);
 %
 ord=input(['Desired Number of Singular Values (positive integer): ']);
 if ( isempty(ord)==1 | ord==0 )
  ord=nmax;
 end
 sigkpt=sprintf('%g',100*sum(svl(1:ord))/sum(svl));
 disp([' ']);
 disp([' # HSV set to: ',sprintf('%i',ord)]);
 disp([' Model Describes, ',sigkpt,' (%) of Data']);
 disp([' ']);
else
 ord=abs(ord);
end
%
% joint diagonalization
%
[L,R,cost]=wsf_AB(Hx,ord); % L - left, R - right
%
W=pinv(R);
V=pinv(L);
%
%%

Contact us