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