--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/eigcsd.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,249 @@
+% EIGCSD calculates TFs from 2D cross-correlated spectra.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION:
+% 
+%     Calculates TFs or WFs from 2dim cross correlated spectra. Input the
+%     elemnts of a cross-spectral matrix
+% 
+%                / csd11(f)  csd12(f) \
+%      CSD(f) =  |                    |
+%                \ csd21(f)  csd22(f) /
+% 
+%     and output the frequency response of four innovation transfer
+%     functions or four whitening filters that can be used to color white
+%     noise or to whitening colored noise.
+% 
+% CALL: [h11,h12,h21,h22] = eigcsd(csd11,csd12,csd21,csd22,varargin)
+%       [h11,h12,h21,h22] = eigcsd(csd11,csd12,[],csd22,varargin)
+%       [h11,h12,h21,h22] = eigcsd(csd11,[],csd21,csd22,varargin)
+% 
+% INPUT:
+% 
+%     csd11, csd12, csd21 and csd22 are the elements of the cross spectral
+%     matrix
+%     Input also the parameters specifying calculation options
+%     'USESYM' define the calculation method, allowed values are 0
+%     (double precision arithmetic) and 1 (symbolic toolbox variable
+%     precision arithmetic)
+%     'DIG' define the digits used in VPA calculation
+%     'OTP' define the output type. Allowed values are 'TF' output the
+%     transfer functions or 'WF' output the whitening filters frequency
+%     responses
+%     'KEEPVAR' get a whitening filter that preserve the variance of
+%     the input data. Values are true or false
+%     'VARS' first and second channels variance
+% 
+% OUTPUT:
+% 
+%     h11, h12, h21 and h22 are the four innovation TFs frequency responses
+%     or the four whitening filters frequency responses
+% 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% HISTORY:     02-10-2008 L Ferraioli
+%                 Creation
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% VERSION: $Id: eigcsd.m,v 1.8 2009/06/16 16:06:50 luigi Exp $
+%
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [h11,h12,h21,h22] = eigcsd(csd11,csd12,csd21,csd22,varargin)
+
+  % Init
+
+  S11 = csd11;
+  S22 = csd22;
+  
+  if isempty(csd12) && isempty(csd21)
+    error(' You must input csd12 or csd21 ')
+  end
+  
+  if isempty(csd12)
+    S12 = conj(csd21);
+  else
+    S12 = csd12;
+  end
+  if isempty(csd21)
+    S21 = conj(csd12);
+  else
+    S21 = csd21;
+  end
+  
+  npts = length(S11);
+
+  % Finding parameters
+
+  % default
+  usesym = 0;
+  dig = 50;
+  otp = 'TF';
+  kv = false;
+  vars = [1 1];
+%   fs = 10;
+
+  if ~isempty(varargin)
+    for j=1:length(varargin)
+      if strcmp(varargin{j},'USESYM')
+        usesym = varargin{j+1};
+      end
+      if strcmp(varargin{j},'DIG')
+        dig = varargin{j+1};
+      end
+      if strcmp(varargin{j},'OTP')
+        otp = varargin{j+1};
+      end
+      if strcmp(varargin{j},'KEEPVAR')
+        kv = varargin{j+1};
+      end
+      if strcmp(varargin{j},'VARS')
+        vars = varargin{j+1};
+      end
+    end
+  end
+
+  % Defining calculation method
+  if usesym == 1
+    method = 'VPA';
+  else
+    method = 'NUM';
+  end
+
+  % Finding suppression
+
+  k = min(sqrt(S11./S22));
+  if k>=1
+    suppr = floor(k);
+  else
+    n=0;
+    while k<1
+      k=k*10;
+      n=n+1;
+    end
+    k = floor(k);
+    suppr = k*10^(-n);
+  end
+
+  supmat = [1 0;0 suppr];
+  % isupmat = [1 0;0 1/suppr];
+  isupmat = inv(supmat);
+  
+  % check for keep variance option
+  if isequal(otp,'WF') && kv
+    % get mean of spectra
+    s1 = sqrt(vars(1));
+    s2 = sqrt(vars(2));
+  end
+
+  % Core Calculation
+
+  switch method
+    case 'NUM'
+
+      % initializing output data
+      h11 = ones(npts,1);
+      h12 = ones(npts,1);
+      h21 = ones(npts,1);
+      h22 = ones(npts,1);
+      
+      % T = zeros(2,2,npts);
+      % D = zeros(2,2,npts);
+
+      for phi = 1:npts
+
+        % Appliing suppression
+        PP = supmat*[S11(phi) S12(phi);S21(phi) S22(phi)]*supmat;
+        
+        [V,D] = eig(PP);
+        
+        % 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
+                
+        % Definition of the transfer functions
+        switch otp
+          case 'TF'
+            HH = isupmat*[Vp(:,1) Vp(:,2)]*[sqrt(Lp(1,1)) 0;0 sqrt(Lp(2,2))];
+          case 'WF'
+            HH = [1/sqrt(Lp(1,1)) 0;0 1/sqrt(Lp(2,2))]*inv(Vp)*supmat;
+            % HH = [1/sqrt(Lp(1,1)) 0;0 1/sqrt(Lp(2,2))]*(Vp')*supmat;
+        end
+
+        if isequal(otp,'WF') && kv
+          h11(phi,1) = s1*HH(1,1);
+          h12(phi,1) = s1*HH(1,2);
+          h21(phi,1) = s2*HH(2,1);
+          h22(phi,1) = s2*HH(2,2);
+        else
+          h11(phi,1) = HH(1,1);
+          h12(phi,1) = HH(1,2);
+          h21(phi,1) = HH(2,1);
+          h22(phi,1) = HH(2,2);
+        end
+        
+        % T(:,:,phi) = conj(HH)*[S11(phi) S12(phi);S21(phi) S22(phi)]*HH.';
+        % D(:,:,phi) = conj(HH)*HH.' - [S11(phi) S12(phi);S21(phi) S22(phi)];
+
+      end
+
+    case 'VPA'
+      % Define the numerical precision
+      digits(dig)
+
+      % initializing output data
+      h11 = vpa(ones(npts,1));
+      h12 = vpa(ones(npts,1));
+      h21 = vpa(ones(npts,1));
+      h22 = vpa(ones(npts,1));
+
+      SS11 = vpa(S11);
+      SS12 = vpa(S12);
+      SS21 = vpa(S21);
+      SS22 = vpa(S22);
+      
+      for phi = 1:npts
+
+        % Appliing suppression
+        PP = supmat*[SS11(phi) SS12(phi);SS21(phi) SS22(phi)]*supmat;
+        
+        [V,D] = eig(PP);
+        Vp1 = V(:,1);
+        Vp2 = -1.*V(:,2);
+        Lp1 = D(1,1);
+        Lp2 = D(2,2);
+
+        % Definition of the transfer functions
+        switch otp
+          case 'TF'
+            HH = isupmat*[Vp1 Vp2]*[sqrt(Lp1) 0;0 sqrt(Lp2)];
+          case 'WF'
+            % HH = [1/sqrt(Lp1) 0;0 1/sqrt(Lp2)]*inv([Vp1 Vp2])*supmat;
+            HH = inv(isupmat*[Vp1 Vp2]*[sqrt(Lp1) 0;0 sqrt(Lp2)]);
+        end
+
+        if isequal(otp,'WF') && kv
+          h11(phi,1) = s1*HH(1,1);
+          h12(phi,1) = s1*HH(1,2);
+          h21(phi,1) = s2*HH(2,1);
+          h22(phi,1) = s2*HH(2,2);
+        else
+          h11(phi,1) = HH(1,1);
+          h12(phi,1) = HH(1,2);
+          h21(phi,1) = HH(2,1);
+          h22(phi,1) = HH(2,2);
+        end
+        
+      end
+      h11 = double(h11);
+      h12 = double(h12);
+      h21 = double(h21);
+      h22 = double(h22);
+  end
+