0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 function [Tinit,Tprop,E] = ngsetup_vpa(den,fs, ndigits)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
2
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
3 % ALGONAME = mfilename;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
4 % VERSION = '$Id: ngsetup_vpa.m,v 1.2 2008/08/01 13:19:42 ingo Exp $';
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
5 digits(ndigits)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
6 fs = sym(fs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
7 den = sym(den,'d');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
8 den=den';
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
9 dt = vpa(1/fs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
10
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
11 n = length(den)-1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
12 % digits(d);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
13
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
14 %% setting up matrix Aij
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
15 m_a = vpa(zeros(n,n));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
16 for i = 1:n
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
17 for j = 1:n
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
18 if j == i+1
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
19 m_a(i,j) = 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
20 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21 if i == n
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22 m_a(i,j) = -den(j);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27 %% Matrix exponential E
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28 a = m_a*dt;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29 E = expm(a);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32 %% setting up matrix Bij
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33 B = vpa(zeros(n,n));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34 for i=1:n
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
35 if rem(i,2) ~= 0
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
36 j0 = (i+1)/2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
37 s = (-1)^(j0+1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
38 j = j0;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
39 for k=1:2:(n+1)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40 d1 = den(k);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41 d2 = s*d1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42 B(i,j) = d2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43 s = -s;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44 j = j+1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
46 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47 if rem(i,2) == 0
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48 j0 = i/2+1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49 s = (-1)^j0;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50 j = j0;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51 for k=2:2:(n+1)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
52 d1 = den(k);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
53 d2 = s*d1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
54 B(i,j) = d2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
55 s = -s;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
56 j = j+1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
57 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
58 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
59 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
60
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
61 %% solve B * m = k
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
62 m_k = vpa(zeros(n,1));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
63 m_k(n) = 0.5;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
64 m_m = vpa(B\m_k);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
65
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
66 %% filling covariance matrix Cinit
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
67 % Cinit = vpa(zeros(n,n));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
68 % for i=1:n
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
69 % for j=1:n
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
70 % if rem((i+j),2) == 0 % even
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
71 % d1 = (-1)^((i-j)/2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
72 % d1 = subs(d1)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
73 % d2 = vpa('m_m((i+j)/2)');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
74 % d2 = subs(d2)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
75 % d3 = vpa(ctranspose(d2));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
76 % d3 = subs(d3)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
77 % d4 = vpa('d1 * d3');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
78 % d4 = subs(d4)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
79 % Cinit(i,j) = d4;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
80 % else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
81 % Cinit(i,j) = 0;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
82 % end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
83 % end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
84 % end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
85
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
86 Cinit = vpa(zeros(n,n));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
87 for i=1:n
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
88 for j=1:n
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
89 if rem((i+j),2) == 0 % even
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
90 d1 = (-1)^((i-j)/2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
91 d2 = (i+j)/2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
92 Cinit(i,j) = d1 * m_m(d2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
93 % else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
94 % Cinit(i,j) = 0;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
95 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
96 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
97 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
98 %cholesky decomposition
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
99
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
100 % Tinit = chol(double(Cinit),'lower'); %lower triangular matrix
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
101 Tinit = ao.mchol(Cinit);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
102
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
103 %% setting up matrix D
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
104 N = n*(n+1)/2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
105
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
106 m_d = vpa(zeros(N));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
107 g = zeros(n);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
108 for i=1:n
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
109 for j=1:n
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
110 if i>=j
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
111 g(i,j) = (i*i-i)/2+(j);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
112 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
113 g(i,j) = (j*j-j)/2+(i);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
114 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
115 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
116 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
117
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
118 for i=1:n
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
119 for j=i:n
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
120 for k=1:n
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
121 m_d(g(i,j),g(j,k)) = m_d(g(i,j),g(j,k)) + m_a(i,k);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
122 m_d(g(i,j),g(i,k)) = m_d(g(i,j),g(i,k)) + m_a(j,k);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
123 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
124 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
125 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
126
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
127
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
128 %% setting up q from D * p = q
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
129 m_q = vpa(zeros(1,g(n,n)));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
130 for i=1:n
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
131 for j=i:n
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
132 if i==n
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
133 m_q(g(i,j)) = (E(n,n))*(E(n,n))-1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
134 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
135 m_q(g(i,j)) = (E(i,n)*E(j,n));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
136 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
137 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
138 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
139
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
140 m_p = m_d\m_q';
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
141
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
142 Cprop = vpa(zeros(n));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
143 for i=1:n
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
144 for j=1:n
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
145 Cprop(i,j) = m_p(g(i,j));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
146 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
147 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
148
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
149 Tprop = ao.mchol(Cprop);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
150 Tprop = double(Tprop);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
151 E = double(E);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
152 Tinit = double(Tinit);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
153
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
154 % Tprop = chol(Cprop,'lower');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
155
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
156 % Tprop = chol(double((Cprop)),'lower');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
157 % Tprop = mchol(Cprop);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
158
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
159 % %% writing the generator
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
160 % r = randn(n,1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
161 % y = Tinit * r;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
162 % x = zeros(Nfft,1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
163 % for i=1:Nfft
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
164 % r = randn(n,1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
165 % y = E * y + Tprop * r;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
166 % x(i) = a*y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
167 % end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
168
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
169 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
170
|