--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@matrix/mchNoisegen.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,560 @@
+% MCHNOISEGEN Generates multichannel noise data series given a model
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% FUNCTION:    mchNoisegen generates multichannel noise data series given a
+% mutichannel noise generating filter.
+%
+% DESCRIPTION: 
+% 
+%
+% CALL:        out = mchNoisegen(mod, pl)
+%
+% INPUT:
+%         mod: is a matrix containing a multichannel noise generating
+%         filter aiming to generate colored noise from unitary variance
+%         independent white data series. Each element of the multichannel
+%         filter must be a parallel bank of filters as that produced by
+%         matrix/mchNoisegenFilter or ao/Noisegen2D. The filter matrix must
+%         be square.
+% 
+% OUTPUT:
+%         out: is a matrix containg a multichannel colored noise time
+%         series which csd matrix is defined by mod'*mod.
+%
+% HISTORY:     19-10-2009 L Ferraioli
+%              Creation
+%
+% ------------------------------------------------------------------------
+%
+% <a href="matlab:utils.helper.displayMethodInfo('matrix', 'mchNoisegen')">Parameters Description</a>
+%
+% VERSION:     $Id: mchNoisegen.m,v 1.6 2011/04/08 08:56:31 hewitson Exp $
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function varargout = mchNoisegen(varargin)
+
+  % Check if this is a call for parameters
+  if utils.helper.isinfocall(varargin{:})
+    varargout{1} = getInfo(varargin{3});
+    return
+  end
+  
+  import utils.const.*
+  utils.helper.msg(msg.OMNAME, 'running %s/%s', mfilename('class'), mfilename);
+  
+  % Collect input variable names
+  in_names = cell(size(varargin));
+  for ii = 1:nargin,in_names{ii} = inputname(ii);end
+  
+  % Collect all ltpdauoh objects
+  [mtxs, mtxs_invars] = utils.helper.collect_objects(varargin(:), 'matrix', in_names);
+  [pl, invars] = utils.helper.collect_objects(varargin(:), 'plist');
+  
+  inhists = mtxs.hist;
+  
+  % combine plists
+  pl = parse(pl, getDefaultPlist());
+  pl.getSetRandState();
+  
+  outm = copy(mtxs,nargout);
+  
+  % Get parameters and set params for fit
+  Nsecs       = find(pl, 'nsecs');
+  fs          = find(pl, 'fs');
+  yunit       = find(pl, 'yunits');
+  
+  % total number of data in the series
+  Ntot = round(Nsecs*fs);
+  
+  % chose case for input filter
+  if isa(outm,'matrix') && isa(outm.objs(1),'filterbank')
+    % discrete system
+    sys = 'z';
+    mfil = outm.objs;
+    [nn,mm] = size(mfil);
+    if nn~=mm
+      error('!!! Filter matrix must be square')
+    end
+    % check if filter is a matrix built with noisegen2D
+    fromnsg2D = false;
+    if nn == 2
+      nn11 = numel(mfil(1,1).filters);
+      nn21 = numel(mfil(2,1).filters);
+      if nn11 == nn21
+        % get poles out of the filters
+        pl11 = zeros(nn11,1);
+        pl21 = zeros(nn21,1);
+        for ii=1:nn11
+          pl11(ii) = -1*mfil(1,1).filters(ii).b(2);
+          pl21(ii) = -1*mfil(2,1).filters(ii).b(2);
+        end
+        % check if poles are equal that means they are produced with
+        % noisegen2D
+        if all((pl11-pl21)<eps)
+          fromnsg2D = true;
+        end
+      end  
+    end
+  elseif isa(outm,'matrix') && isa(outm.objs(1),'parfrac')
+    % continuous system
+    sys = 's';
+    mfil = outm.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
+      
+      % init output data
+      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;
+        
+        cns = getinitz(res,pls,filtsz,fromnsg2D);
+        
+        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
+      
+      outm = 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
+      
+      outm = matrix(out);
+      
+  end
+  
+  % set output
+  varargout{1} = outm;
+  
+
+end
+
+%--------------------------------------------------------------------------
+% Get Info Object
+%--------------------------------------------------------------------------
+function ii = getInfo(varargin)
+  if nargin == 1 && strcmpi(varargin{1}, 'None')
+    sets = {};
+    pl   = [];
+  else
+    sets = {'Default'};
+    pl   = getDefaultPlist;
+  end
+  % Build info object
+  ii = minfo(mfilename, 'matrix', 'ltpda', utils.const.categories.sigproc, '$Id: mchNoisegen.m,v 1.6 2011/04/08 08:56:31 hewitson Exp $', sets, pl);
+  ii.setArgsmin(1);
+  ii.setOutmin(1);
+  ii.setOutmax(1);
+end
+
+%--------------------------------------------------------------------------
+% Get Default Plist
+%--------------------------------------------------------------------------
+function plout = getDefaultPlist()
+  persistent pl;  
+  if exist('pl', 'var')==0 || isempty(pl)
+    pl = buildplist();
+  end
+  plout = pl;  
+end
+
+function pl = buildplist()
+  
+  pl = plist();
+  
+  % Nsecs
+  p = param({'nsecs', 'Number of seconds in the desired noise data series.'}, {1, {1}, paramValue.OPTIONAL});
+  pl.append(p);
+  
+  % Fs
+  p = param({'fs', 'The sampling frequency of the noise data series.'}, {1, {1}, paramValue.OPTIONAL});
+  pl.append(p);
+  
+  % Yunits
+  p = param({'yunits','Unit on Y axis.'},  paramValue.STRING_VALUE(''));
+  pl.append(p);
+  
+end
+
+%--------------------------------------------------------------------------
+% Local function
+% Estimate init values for the z domain case
+%--------------------------------------------------------------------------
+function cns = getinitz(res,pls,filtsz,fromnsg2D)
+
+  if fromnsg2D % init in the case of filters produced with noisegen2D
+    % divide in 2 the problem
+    res1 = res(1:filtsz(1));
+    res2 = res(filtsz(1)+1:filtsz(2));
+    pls1 = pls(1:filtsz(1));
+    pls2 = pls(filtsz(1)+1:filtsz(2));
+    
+    %%% the first series %%%
+    Nrs = numel(res1);
+    % get covariance matrix
+    R = zeros(Nrs);
+    for aa=1:Nrs
+      for bb=1:Nrs
+        R(aa,bb) = (res1(aa)*conj(res1(bb)))/(1-pls1(aa)*conj(pls1(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
+    cns1 = A*ns;
+
+    % need to correct for roundoff errors
+    % cleaning up results for numerical approximations
+    idx = imag(pls1(:,1))==0;
+    cns1(idx) = real(cns1(idx));
+
+    % cleaning up results for numerical roundoff errors
+    % states associated to complex conjugate poles must be complex
+    % conjugate
+    idxi = imag(pls1(:,1))~=0;
+    icns1 = cns1(idxi);
+    for jj = 1:2:numel(icns1)
+      icns1(jj+1,1) = conj(icns1(jj,1));
+    end
+    cns1(idxi) = icns1;
+    
+    %%% the second series %%%
+    Nrs = numel(res2);
+    % get covariance matrix
+    R = zeros(Nrs);
+    for aa=1:Nrs
+      for bb=1:Nrs
+        R(aa,bb) = (res2(aa)*conj(res2(bb)))/(1-pls2(aa)*conj(pls2(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
+    cns2 = A*ns;
+
+    % need to correct for roundoff errors
+    % cleaning up results for numerical approximations
+    idx = imag(pls2(:,1))==0;
+    cns2(idx) = real(cns2(idx));
+
+    % cleaning up results for numerical roundoff errors
+    % states associated to complex conjugate poles must be complex
+    % conjugate
+    idxi = imag(pls2(:,1))~=0;
+    icns2 = cns2(idxi);
+    for jj = 1:2:numel(icns2)
+      icns2(jj+1,1) = conj(icns2(jj,1));
+    end
+    cns2(idxi) = icns2;
+    
+    %%% combine results %%%
+    cns = [cns1;cns2];
+    
+  else % init in the case of filters produced with matrix constructor
+    
+    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;
+    
+  end
+
+end