Code covered by the BSD License

# Blind Modal Identification with Hankel Matrices (BMIDHM) toolbox

### Scot McNeill (view profile)

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);
%
%%```