% 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