0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 function [Tinit,Tprop,E] = ngsetup(den,fs)
0020
0021 den=den';
0022 dt = 1/fs;
0023 length(den);
0024
0025 n = length(den)-1;
0026
0027
0028
0029 m_a = zeros(n,n);
0030 for i = 1:n
0031 for j = 1:n
0032 if j == i+1
0033 m_a(i,j) = 1;
0034 end
0035 if i == n
0036 m_a(i,j) = -den(j);
0037 end
0038 end
0039 end
0040
0041
0042 E = expm(m_a*dt);
0043
0044
0045 B = zeros(n,n);
0046 for i=1:n
0047 if rem(i,2) ~= 0
0048 j0 = (i+1)/2;
0049 s = (-1)^(j0+1);
0050 j = j0;
0051 for k=1:2:(n+1)
0052 B(i,j) = s*den(k);
0053 s = -s;
0054 j = j+1;
0055 end
0056 end
0057 if rem(i,2) == 0
0058 j0 = i/2+1;
0059 s = (-1)^j0;
0060 j = j0;
0061 for k=2:2:(n+1)
0062 B(i,j) = s * den(k);
0063 s = -s;
0064 j = j+1;
0065 end
0066 end
0067 end
0068
0069
0070 m_k = zeros(n,1);
0071 m_k(n) = 0.5;
0072 m_m = B\m_k;
0073
0074
0075 Cinit = zeros(n,n);
0076 for i=1:n
0077 for j=1:n
0078 if rem((i+j),2) == 0
0079 Cinit(i,j) = (-1)^((i-j)/2) * m_m((i+j)/2);
0080 else
0081 Cinit(i,j) = 0;
0082 end
0083 end
0084 end
0085
0086
0087
0088
0089 Tinit = chol(Cinit,'lower');
0090
0091
0092 N = n*(n+1)/2;
0093
0094 m_d = zeros(N);
0095 g = zeros(n);
0096 for i=1:n
0097 for j=1:n
0098 if i>=j
0099 g(i,j) = (i*i-i)/2+(j);
0100 else
0101 g(i,j) = (j*j-j)/2+(i);
0102 end
0103 end
0104 end
0105
0106 for i=1:n
0107 for j=i:n
0108 for k=1:n
0109 m_d(g(i,j),g(j,k)) = m_d(g(i,j),g(j,k)) + m_a(i,k);
0110 m_d(g(i,j),g(i,k)) = m_d(g(i,j),g(i,k)) + m_a(j,k);
0111 end
0112 end
0113 end
0114
0115
0116
0117 m_q = zeros(1,g(n,n));
0118 for i=1:n
0119 for j=i:n
0120 if i==n
0121 m_q(g(i,j)) = (E(n,n))*(E(n,n))-1;
0122 else
0123 m_q(g(i,j)) = (E(i,n)*E(j,n));
0124 end
0125 end
0126 end
0127
0128 m_p = m_d\m_q';
0129
0130 Cprop = zeros(n);
0131 for i=1:n
0132 for j=1:n
0133 Cprop(i,j) = m_p(g(i,j));
0134 end
0135 end
0136 Tprop = chol(Cprop,'lower');
0137
0138
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148 end
0149