Mercurial > hg > ltpda

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