% NDEIGCSD calculates TFs from ND cross-correlated spectra.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DESCRIPTION:%% Calculates TFs or WFs from Ndim cross correlated spectra. Input% elemnts of a cross-spectral matrix in 2D are assumed to be:%% / csd11(f) csd12(f) \% CSD(f) = | |% \ csd21(f) csd22(f) /%%% CALL: h = eigcsd(csd,varargin)%% INPUT:%% csd are the elements of the cross spectra matrix. It is a (n,n,m)% matrix where n is the dimensionality of the system and m is the% number of frequency samples% % Input also the parameters specifying calculation options% % 'OTP' define the output type. Allowed values are 'TF' output the% transfer functions or 'WF' output the whitening filters frequency% responses. Default 'TF'% 'MTD' define the method for the calculation of the csd matrix of a% multichannel system. Admitted values are 'PAP' referring to Papoulis% [1] style calculation in which csd = TF*I*TF' and 'KAY' referring to% Kay [2] style calculation in which csd = conj(TF)*I*TF.'.% Default 'PAP'%% OUTPUT:%% h are the TFs or WFs frequency responses. It is a (n,n,m) matrix in% which n is the dimensionality of the system and m is the number of% frequency samples.% % REFERENCES:% % [1] A. Papoulis, Probability Random Variable and Stochastic Processes,% McGraw-Hill, third edition, 1991.% [2] S. M. Kay, Modern Spectral Estimation, Prentice Hall, 1988.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% HISTORY: 23-04-2009 L Ferraioli% Creation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% VERSION: $Id: ndeigcsd.m,v 1.4 2009/11/06 16:55:51 luigi Exp $%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function h = ndeigcsd(csd,varargin) [l,m,npts] = size(csd); if m~=l error('!!! The first two dimensions of csd must be equal. csd must be a square matrix frequency by frequency.') end % Finding parameters % default otp = 'TF'; mtd = 'PAP'; if ~isempty(varargin) for j=1:length(varargin) if strcmp(varargin{j},'OTP') otp = upper(varargin{j+1}); end if strcmp(varargin{j},'MTD') mtd = upper(varargin{j+1}); end end end % Finding suppression suppr = ones(l,l); for ii = 2:l k = ones(l,1); for jj = ii-1:-1:1 k(jj) = min(sqrt(csd(jj,jj,:)./csd(ii,ii,:))); if k(jj)>=1 suppr(jj,ii) = floor(k(jj)); else n=0; while k(jj)<1 k(jj)=k(jj)*10; n=n+1; end k(jj) = floor(k(jj)); suppr(jj,ii) = k(jj)*10^(-n); end end % csuppr(ii) = prod(suppr(:,ii)); end csuppr = prod(suppr,2); supmat = diag(csuppr); ssup = supmat*supmat.'; supmat = rot90(rot90(supmat)); isupmat = inv(supmat); % Core Calculation % initializing output dat h = ones(l,m,npts); for phi = 1:npts % Appliing suppression PP = csd(:,:,phi); PP = supmat*PP*supmat;% [V,D] = eig(PP,ssup); [V,D,U] = svd(PP,0);% [V,D] = eig(PP); % Correcting the output of eig% V = fliplr(V);% D = rot90(rot90(D));% % Correcting the output of eig% Vp = fliplr(V);% Lp = rot90(rot90(D));% % % Correcting the phase% [a,b] = size(PP);% for ii=1:b% Vp(:,ii) = Vp(:,ii).*(cos(angle(Vp(ii,ii)))-1i*sin(angle(Vp(ii,ii))));% Vp(ii,ii) = real(Vp(ii,ii));% end% % V = Vp;% D = Lp; % Definition of the transfer functions switch otp case 'TF' switch mtd case 'PAP'% HH = ssup*V*sqrt(D); HH = isupmat*V*sqrt(D);% HH = V*sqrt(D); case 'KAY'% HH = conj(ssup*V*sqrt(D)); HH = conj(isupmat*V*sqrt(D));% HH = conj(V*sqrt(D)); end case 'WF' switch mtd case 'PAP' %HH = inv(ssup*V*sqrt(D)); HH = inv(isupmat*V*sqrt(D)); case 'KAY' %HH = inv(conj(ssup*V*sqrt(D))); HH = inv(conj(isupmat*V*sqrt(D))); end end h(:,:,phi) = HH; endend