--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@matrix/MultiChannelNoise.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,303 @@
+% MULTICHANNELNOISE
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% FUNCTION:    MultiChannelNoise
+%
+% DESCRIPTION: 
+% 
+%
+% CALL:        a = MultiChannelNoise(a, pl)
+%
+% PARAMETER:
+%              pl:       plist containing Nsecs, fs
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function a = MultiChannelNoise(a, pli)
+
+  VERSION = '$Id: MultiChannelNoise.m,v 1.10 2010/10/29 16:09:13 ingo Exp $';
+
+  % get AO info
+  iobj = matrix.getInfo('matrix', 'Multichannel Noise');
+
+  % Set the method version string in the minfo object
+  iobj.setMversion([VERSION '-->' iobj.mversion]);
+
+  % Add default values
+  pl = parse(pli, iobj.plists);
+  
+  % Get parameters and set params for fit
+  Nsecs       = find(pl, 'nsecs');
+  fs          = find(pl, 'fs');
+  filter      = find(pl, 'model');
+  yunit       = find(pl, 'yunits');
+  
+  % total number of data in the series
+  Ntot = round(Nsecs*fs);
+  
+  % chose case for input filter
+  if isa(filter,'matrix') && isa(filter.objs(1),'filterbank')
+    % discrete system
+    sys = 'z';
+    mfil = filter.objs;
+    [nn,mm] = size(mfil);
+    if nn~=mm
+      error('!!! Filter matrix must be square')
+    end
+  elseif isa(filter,'matrix') && isa(filter.objs(1),'parfrac')
+    % continuous system
+    sys = 's';
+    mfil = filter.objs;
+    [nn,mm] = size(mfil);
+    if nn~=mm
+      error('!!! Filter matrix must be square')
+    end
+  else
+    error('!!! Input filter must be a ''matrix'' of ''filterbank'' or ''parfrac'' objects.')
+  end
+  
+  % init output
+  out(nn,1) = ao;
+  
+  % switch between input filter type
+  switch sys
+    case 'z' % discrete input filters
+      
+      o = zeros(nn,Ntot);
+      
+      for zz=1:nn % moving along system dimension
+        
+        % extract residues and poles from input objects
+        % run over matrix dimensions
+        res = [];
+        pls = [];
+        filtsz = [];
+        for jj=1:nn % collect filters coefficients along the columns zz
+          bfil = mfil(jj,zz).filters;
+          filtsz = [filtsz; numel(bfil)];
+          for kk=1:numel(bfil)
+            num = bfil(kk).a;
+            den = bfil(kk).b;
+            res = [res; num(1)];
+            pls = [pls; -1*den(2)];
+          end
+        end
+        
+        % rescaling residues to get the correct result for univariate psd
+        res = res.*sqrt(fs/2);
+        
+        Nrs = numel(res);
+        % get covariance matrix
+        R = zeros(Nrs);
+        for aa=1:Nrs
+          for bb=1:Nrs
+            R(aa,bb) = (res(aa)*conj(res(bb)))/(1-pls(aa)*conj(pls(bb)));
+          end
+        end
+        
+        % avoiding problems caused by roundoff errors
+        HH = triu(R,0); % get upper triangular part of R
+        HH1 = triu(R,1); % get upper triangular part of R above principal diagonal
+        HH2 = HH1'; % do transpose conjugate
+        R = HH + HH2; % reconstruct R in order to be really hermitian
+        
+        % get singular value decomposition of R
+        [U,S,V] = svd(R,0);
+
+        % conversion matrix
+        A = V*sqrt(S);
+
+        % generate unitary variance gaussian random noise
+        %ns = randn(Nrs,Ntot);
+        ns = randn(Nrs,1);
+
+        % get correlated starting data points
+        cns = A*ns;
+        
+        % need to correct for roundoff errors
+        % cleaning up results for numerical approximations
+        idx = imag(pls(:,1))==0;
+        cns(idx) = real(cns(idx));
+        
+        % cleaning up results for numerical roundoff errors
+        % states associated to complex conjugate poles must be complex
+        % conjugate
+        idxi = imag(pls(:,1))~=0;
+        icns = cns(idxi);
+        for jj = 1:2:numel(icns)
+          icns(jj+1,1) = conj(icns(jj,1));
+        end
+        cns(idxi) = icns;
+        
+        y = zeros(sum(filtsz),2);
+        rnoise = zeros(sum(filtsz),1);
+        rns = randn(1,Ntot);
+        %rns = utils.math.blwhitenoise(Ntot,fs,1/Nsecs,fs/2);
+        %rns = rns.'; % willing to work with raw
+        
+        y(:,1) = cns;
+        
+        % start recurrence
+        for xx = 2:Ntot+1
+          rnoise(:,1) = rns(xx-1);
+          y(:,2) = pls.*y(:,1) + res.*rnoise;
+          idxr = 0;
+          for kk=1:nn
+            o(kk,xx-1) = o(kk,xx-1) + sum(y(idxr+1:idxr+filtsz(kk),2));
+            idxr = idxr+filtsz(kk);
+          end
+          y(:,1) = y(:,2);
+        end
+        
+      end
+      
+      clear rns
+      
+      % build output ao
+      for dd=1:nn
+        out(dd,1) = ao(tsdata(o(dd,:),fs));
+        out(dd,1).setYunits(unit(yunit));
+      end
+      
+      a = matrix(out);
+      
+    case 's' % continuous input filters
+      
+      o = zeros(nn,Ntot);
+      
+      T = 1/fs; % sampling period
+      
+      for zz=1:nn % moving along system dimension
+        
+        % extract residues and poles from input objects
+        % run over matrix dimensions
+        res = [];
+        pls = [];
+        filtsz = [];
+        for jj=1:nn % collect filters coefficients along the columns zz
+          bfil = mfil(jj,zz);
+          filtsz = [filtsz; numel(bfil.res)];
+          res = [res; bfil.res];
+          pls = [pls; bfil.poles];
+        end
+        
+        % rescaling residues to get the correct result for univariate psd
+        res = res.*sqrt(fs/2);
+        
+        Nrs = numel(res);
+        
+        % get covariance matrix for innovation
+        Ri = zeros(Nrs);
+        for aa=1:Nrs
+          for bb=1:Nrs
+            Ri(aa,bb) = (res(aa)*conj(res(bb)))*(exp((pls(aa) + conj(pls(bb)))*T)-1)/(pls(aa) + conj(pls(bb)));
+          end
+        end
+        
+        % avoiding problems caused by roundoff errors
+        HH = triu(Ri,0); % get upper triangular part of R
+        HH1 = triu(Ri,1); % get upper triangular part of R above principal diagonal
+        HH2 = HH1'; % do transpose conjugate
+        Ri = HH + HH2; % reconstruct R in order to be really hermitian
+        
+        % get singular value decomposition of R
+        [Ui,Si,Vi] = svd(Ri,0);
+
+        % conversion matrix for innovation
+        Ai = Vi*sqrt(Si);
+        
+        % get covariance matrix for initial state
+        Rx = zeros(Nrs);
+        for aa=1:Nrs
+          for bb=1:Nrs
+            Rx(aa,bb) = -1*(res(aa)*conj(res(bb)))/(pls(aa) + conj(pls(bb)));
+          end
+        end
+        
+        % avoiding problems caused by roundoff errors
+        HH = triu(Rx,0); % get upper triangular part of R
+        HH1 = triu(Rx,1); % get upper triangular part of R above principal diagonal
+        HH2 = HH1'; % do transpose conjugate
+        Rx = HH + HH2; % reconstruct R in order to be really hermitian
+        
+        % get singular value decomposition of R
+        [Ux,Sx,Vx] = svd(Rx,0);
+
+        % conversion matrix for initial state
+        Ax = Vx*sqrt(Sx);
+        
+        % generate unitary variance gaussian random noise
+        %ns = randn(Nrs,Ntot);
+        ns = randn(Nrs,1);
+
+        % get correlated starting data points
+        cns = Ax*ns;
+        
+        % need to correct for roundoff errors
+        % cleaning up results for numerical approximations
+        idx = imag(pls(:,1))==0;
+        cns(idx) = real(cns(idx));
+        
+        % cleaning up results for numerical roundoff errors
+        % states associated to complex conjugate poles must be complex
+        % conjugate
+        idxi = imag(pls(:,1))~=0;
+        icns = cns(idxi);
+        for jj = 1:2:numel(icns)
+          icns(jj+1,1) = conj(icns(jj,1));
+        end
+        cns(idxi) = icns;
+        
+        y = zeros(sum(filtsz),2);
+        rnoise = zeros(sum(filtsz),1);
+        rns = randn(1,Ntot);
+        %rns = utils.math.blwhitenoise(Ntot,fs,1/Nsecs,fs/2);
+        %rns = rns.'; % willing to work with raw
+        
+        y(:,1) = cns;
+        
+        % start recurrence
+        for xx = 2:Ntot+1
+%           innov = Ai*randn(sum(filtsz),1);
+          rnoise(:,1) = rns(xx-1);
+          innov = Ai*rnoise;
+          % need to correct for roundoff errors
+          % cleaning up results for numerical approximations
+          innov(idx) = real(innov(idx));
+
+          % cleaning up results for numerical roundoff errors
+          % states associated to complex conjugate poles must be complex
+          % conjugate
+          iinnov = innov(idxi);
+          for jj = 1:2:numel(iinnov)
+            iinnov(jj+1,1) = conj(iinnov(jj,1));
+          end
+          innov(idxi) = iinnov;
+                   
+          y(:,2) = diag(exp(pls.*T))*y(:,1) + innov;
+          
+          idxr = 0;
+          for kk=1:nn
+            o(kk,xx-1) = o(kk,xx-1) + sum(y(idxr+1:idxr+filtsz(kk),2));
+            idxr = idxr+filtsz(kk);
+          end
+          y(:,1) = y(:,2);
+        end
+        
+      end
+      
+%       clear rns
+      
+      % build output ao
+      for dd=1:nn
+        out(dd,1) = ao(tsdata(o(dd,:),fs));
+        out(dd,1).setYunits(unit(yunit));
+      end
+      
+      a = matrix(out);
+      
+  end
+  
+
+end
+
+