Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/ngsetup.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 % NGSETUP is called by the function fromPzmodel
2 %
3 % Inputs calculated by ...
4 % ... NGCONV:
5 % - den: denominator coefficients
6 %
7 % ... USER
8 % - fs: sampling frequency given as input to LTPDA_NOISEGEN
9 %
10 % Outputs:
11 % - Tinit: matrix to calculate initial state vector
12 % - Tprop: matrix to calculate propagation vector
13 % - E: matrix to calculate propagation vector
14 %
15 % A Monsky 24-07-07
16 %
17 % $Id: ngsetup.m,v 1.3 2008/10/20 08:38:29 anneke Exp $
18
19 function [Tinit,Tprop,E] = ngsetup(den,fs)
20
21 den=den';
22 dt = 1/fs;
23 length(den);
24
25 n = length(den)-1;
26
27 %% setting up matrix Aij
28
29 m_a = zeros(n,n);
30 for i = 1:n
31 for j = 1:n
32 if j == i+1
33 m_a(i,j) = 1;
34 end
35 if i == n
36 m_a(i,j) = -den(j);
37 end
38 end
39 end
40
41 %% Matrix exponential E
42 E = expm(m_a*dt);
43
44 %% setting up matrix Bij
45 B = zeros(n,n);
46 for i=1:n
47 if rem(i,2) ~= 0
48 j0 = (i+1)/2;
49 s = (-1)^(j0+1);
50 j = j0;
51 for k=1:2:(n+1)
52 B(i,j) = s*den(k);
53 s = -s;
54 j = j+1;
55 end
56 end
57 if rem(i,2) == 0
58 j0 = i/2+1;
59 s = (-1)^j0;
60 j = j0;
61 for k=2:2:(n+1)
62 B(i,j) = s * den(k);
63 s = -s;
64 j = j+1;
65 end
66 end
67 end
68
69 %% solve B * m = k
70 m_k = zeros(n,1);
71 m_k(n) = 0.5;
72 m_m = B\m_k;
73
74 %% filling covariance matrix Cinit
75 Cinit = zeros(n,n);
76 for i=1:n
77 for j=1:n
78 if rem((i+j),2) == 0 % even
79 Cinit(i,j) = (-1)^((i-j)/2) * m_m((i+j)/2);
80 else
81 Cinit(i,j) = 0;
82 end
83 end
84 end
85
86
87 %cholesky decomposition
88
89 Tinit = chol(Cinit,'lower'); %lower triangular matrix
90
91 %% setting up matrix D
92 N = n*(n+1)/2;
93
94 m_d = zeros(N);
95 g = zeros(n);
96 for i=1:n
97 for j=1:n
98 if i>=j
99 g(i,j) = (i*i-i)/2+(j);
100 else
101 g(i,j) = (j*j-j)/2+(i);
102 end
103 end
104 end
105
106 for i=1:n
107 for j=i:n
108 for k=1:n
109 m_d(g(i,j),g(j,k)) = m_d(g(i,j),g(j,k)) + m_a(i,k);
110 m_d(g(i,j),g(i,k)) = m_d(g(i,j),g(i,k)) + m_a(j,k);
111 end
112 end
113 end
114
115
116 %% setting up q from D * p = q
117 m_q = zeros(1,g(n,n));
118 for i=1:n
119 for j=i:n
120 if i==n
121 m_q(g(i,j)) = (E(n,n))*(E(n,n))-1;
122 else
123 m_q(g(i,j)) = (E(i,n)*E(j,n));
124 end
125 end
126 end
127
128 m_p = m_d\m_q';
129
130 Cprop = zeros(n);
131 for i=1:n
132 for j=1:n
133 Cprop(i,j) = m_p(g(i,j));
134 end
135 end
136 Tprop = chol(Cprop,'lower');
137
138 % %% writing the generator
139 % r = randn(n,1);
140 % y = Tinit * r;
141 % x = zeros(Nfft,1);
142 % for i=1:Nfft
143 % r = randn(n,1);
144 % y = E * y + Tprop * r;
145 % x(i) = a*y;
146 % end
147
148 end
149