+
−% 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
+
−
+
−