Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/xCovmat.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 %
3 % Compute cross covariance matrix of two data series
4 %
5 % CALL
6 %
7 % Get autocovariance
8 % cmat = xCovmat(x)
9 % cmat = xCovmat(x,[])
10 % cmat = xCovmat(x,[],'cutbefore',10,'cutafter',10)
11 %
12 % Get crosscovariance
13 % cmat = xCovmat(x,y)
14 % cmat = xCovmat(x,y,'cutbefore',10,'cutafter',10)
15 %
16 % INPUT
17 %
18 % - x, y, data series
19 % - cutbefore, followed by the data samples to cut at the starting of the
20 % data series
21 % - cutafter, followed by the data samples to cut at the ending of the
22 % data series
23 %
24 %
25 % L Ferraioli 10-10-2010
26 %
27 % $Id: xCovmat.m,v 1.1 2010/11/16 16:41:37 luigi Exp $
28 %
29 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30
31 function cmat = xCovmat(x,y,varargin)
32
33 acov = true;
34 % willing to work with columns
35 if size(x,1)<size(x,2)
36 x = x.';
37 end
38 % subtract the mean
39 x = x - mean(x);
40 if nargin > 1
41 if ~isempty(y) % crosscovariance
42 acov = false;
43 if size(y,1)<size(y,2)
44 y = y.';
45 end
46 y = y - mean(y);
47 end
48 end
49
50 nx = size(x,1);
51
52 cutbefore = [];
53 cutafter = [];
54 if ~isempty(varargin)
55 for j=1:length(varargin)
56 if strcmp(varargin{j},'cutbefore')
57 cutbefore = varargin{j+1};
58 end
59 if strcmp(varargin{j},'cutafter')
60 cutafter = varargin{j+1};
61 end
62 end
63 end
64
65 if acov % autocovariance
66
67 % init trim matrix
68 trim = zeros(nx,nx);
69 % fillin the trim matrix
70 for ii=1:nx
71 trim(1:nx-ii+1,ii) = x(ii:nx,1);
72 end
73 %trim = sparse(trim);
74 cmat = trim*conj(trim);
75 % normalization
76 %nmat = sparse(triu(ones(nx,nx),0));
77 nmat = triu(ones(nx,nx),0);
78 nmat = rot90(nmat);
79 normat = nmat*nmat;
80 cmat = cmat./normat;
81 %cmat = full(cmat);
82
83 if ~isempty(cutbefore)
84 % cut rows
85 cmat(1:cutbefore,:)=[];
86 % cut columns
87 cmat(:,1:cutbefore)=[];
88 end
89 [nn,mm]=size(cmat);
90 if ~isempty(cutafter)
91 % cut rows
92 cmat(nn-cutafter:nn,:)=[];
93 % cut columns
94 cmat(:,mm-cutafter:mm)=[];
95 end
96
97 else % cross covariance
98
99 % init trim matrix
100 trimx = zeros(nx,nx);
101 trimy = zeros(nx,nx);
102 % fillin the trim matrix
103 for ii=1:nx
104 trimx(1:nx-ii+1,ii) = x(ii:nx,1);
105 trimy(1:nx-ii+1,ii) = y(ii:nx,1);
106 end
107 cmat = trimx*conj(trimy);
108 % normalization
109 nmat = triu(ones(nx,nx),0);
110 nmat = rot90(nmat);
111 normat = nmat*nmat;
112 cmat = cmat./normat;
113
114 if ~isempty(cutbefore)
115 % cut rows
116 cmat(1:cutbefore,:)=[];
117 % cut columns
118 cmat(:,1:cutbefore)=[];
119 end
120 [nn,mm]=size(cmat);
121 if ~isempty(cutafter)
122 % cut rows
123 cmat(nn-cutafter:nn,:)=[];
124 % cut columns
125 cmat(:,mm-cutafter:mm)=[];
126 end
127
128 end
129
130
131
132
133
134
135 end