--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/ndeigcsd.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,173 @@
+% 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;
+    
+  end
+  
+  
+end
\ No newline at end of file