0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 % TEST_MATRIX_LINLSQSVD tests the linlsqsvd method of the AO class.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
3 % L Ferraioli 10-11-2010
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
4 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
5 % $Id: test_matrix_linlsqsvd.m,v 1.1 2011/02/18 17:07:35 luigi Exp $
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
6 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
8
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
9 %% 1) Determine the coefficients of a linear combination of noises and
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
10 %% comapre with lscov:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
11 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
12 % Make some data
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
13 fs = 10;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
14 nsecs = 10;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
15 B1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
16 B1.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
17 B2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
18 B2.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
19 B3 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
20 B3.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21 B4 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22 B4.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24 C1 = matrix(B1,B2,plist('shape',[2,1]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25 C1.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26 C2 = matrix(B3,B4,plist('shape',[2,1]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27 C2.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29 C = matrix([B1 B3;B2 B4]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30 C.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32 n1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33 n2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
35 n = matrix(n1,n2,plist('shape',[2,1]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
36 n.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
37
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
38 a = [ao(1,plist('yunits','m/T')) ao(2,plist('yunits','m/T'))];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
39 A = matrix(a,plist('shape',[2,1]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41 % assign output values
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42 y = C*A + n;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44 %%% Get a fit with linlsqsvd
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45 pobj1 = linlsqsvd(C1, C2, y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
46
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47 % combine results
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48 for ii=1:numel(pobj1.y)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49 prs(ii) = ao(cdata(pobj1.y(ii)));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50 prs(ii).setYunits(pobj1.yunits(ii));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
52 Pars = matrix(prs,plist('shape',[numel(prs),1]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
53 yfit1 = C*Pars;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
54
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
55 %%% do linear combination: using eval
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
56 yfit2 = pobj1.eval;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
57
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
58 % Plot (compare data with fit)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
59 iplot(y.objs(1), yfit1.objs(1), yfit2.objs(1))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
60 iplot(y.objs(2), yfit1.objs(2), yfit2.objs(2))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
61
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
62 %% 2) Determine the coefficients of a linear combination of noises:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
63 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
64 % Make some data
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
65 fs = 10;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
66 nsecs = 10;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
67 x1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
68 x1.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
69 x2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
70 x2.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
71 x3 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
72 x3.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
73 x4 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
74 x4.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
75
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
76 C1 = matrix(x1,x3,plist('shape',[2,1]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
77 C1.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
78 C2 = matrix(x2,x4,plist('shape',[2,1]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
79 C2.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
80
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
81 C = matrix([x1 x2;x3 x4]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
82 C.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
83
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
84 n1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
85 n2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
86 n = matrix(n1,n2,plist('shape',[2,1]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
87 n.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
88
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
89 a = [ao(1,plist('yunits','m/T')) ao(2,plist('yunits','m/m'))];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
90 A = matrix(a,plist('shape',[2,1]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
91 A.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
92
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
93 y = C*A + n;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
94
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
95 %%% Get a fit with linlsqsvd
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
96 pobj1 = linlsqsvd(C1, C2, y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
97
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
98
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
99 % combine results
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
100 for ii=1:numel(pobj1.y)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
101 prs(ii) = ao(cdata(pobj1.y(ii)));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
102 prs(ii).setYunits(pobj1.yunits(ii));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
103 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
104 Pars = matrix(prs,plist('shape',[numel(prs),1]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
105 yfit1 = C*Pars;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
106
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
107 %%% do linear combination: using eval
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
108 yfit2 = pobj1.eval;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
109
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
110 % Plot (compare data with fit)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
111 iplot(y.objs(1), yfit1.objs(1), yfit2.objs(1))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
112 iplot(y.objs(2), yfit1.objs(2), yfit2.objs(2))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
113
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
114
|