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