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