line source
+ − function [Tinit,Tprop,E] = ngsetup_vpa(den,fs, ndigits)
+ −
+ − % ALGONAME = mfilename;
+ − % VERSION = '$Id: ngsetup_vpa.m,v 1.2 2008/08/01 13:19:42 ingo Exp $';
+ − digits(ndigits)
+ − fs = sym(fs);
+ − den = sym(den,'d');
+ − den=den';
+ − dt = vpa(1/fs);
+ −
+ − n = length(den)-1;
+ − % digits(d);
+ −
+ − %% setting up matrix Aij
+ − m_a = vpa(zeros(n,n));
+ − for i = 1:n
+ − for j = 1:n
+ − if j == i+1
+ − m_a(i,j) = 1;
+ − end
+ − if i == n
+ − m_a(i,j) = -den(j);
+ − end
+ − end
+ − end
+ −
+ − %% Matrix exponential E
+ − a = m_a*dt;
+ − E = expm(a);
+ −
+ −
+ − %% setting up matrix Bij
+ − B = vpa(zeros(n,n));
+ − for i=1:n
+ − if rem(i,2) ~= 0
+ − j0 = (i+1)/2;
+ − s = (-1)^(j0+1);
+ − j = j0;
+ − for k=1:2:(n+1)
+ − d1 = den(k);
+ − d2 = s*d1;
+ − B(i,j) = d2;
+ − s = -s;
+ − j = j+1;
+ − end
+ − end
+ − if rem(i,2) == 0
+ − j0 = i/2+1;
+ − s = (-1)^j0;
+ − j = j0;
+ − for k=2:2:(n+1)
+ − d1 = den(k);
+ − d2 = s*d1;
+ − B(i,j) = d2;
+ − s = -s;
+ − j = j+1;
+ − end
+ − end
+ − end
+ −
+ − %% solve B * m = k
+ − m_k = vpa(zeros(n,1));
+ − m_k(n) = 0.5;
+ − m_m = vpa(B\m_k);
+ −
+ − %% filling covariance matrix Cinit
+ − % Cinit = vpa(zeros(n,n));
+ − % for i=1:n
+ − % for j=1:n
+ − % if rem((i+j),2) == 0 % even
+ − % d1 = (-1)^((i-j)/2);
+ − % d1 = subs(d1)
+ − % d2 = vpa('m_m((i+j)/2)');
+ − % d2 = subs(d2)
+ − % d3 = vpa(ctranspose(d2));
+ − % d3 = subs(d3)
+ − % d4 = vpa('d1 * d3');
+ − % d4 = subs(d4)
+ − % Cinit(i,j) = d4;
+ − % else
+ − % Cinit(i,j) = 0;
+ − % end
+ − % end
+ − % end
+ −
+ − Cinit = vpa(zeros(n,n));
+ − for i=1:n
+ − for j=1:n
+ − if rem((i+j),2) == 0 % even
+ − d1 = (-1)^((i-j)/2);
+ − d2 = (i+j)/2;
+ − Cinit(i,j) = d1 * m_m(d2);
+ − % else
+ − % Cinit(i,j) = 0;
+ − end
+ − end
+ − end
+ − %cholesky decomposition
+ −
+ − % Tinit = chol(double(Cinit),'lower'); %lower triangular matrix
+ − Tinit = ao.mchol(Cinit);
+ −
+ − %% setting up matrix D
+ − N = n*(n+1)/2;
+ −
+ − m_d = vpa(zeros(N));
+ − g = zeros(n);
+ − for i=1:n
+ − for j=1:n
+ − if i>=j
+ − g(i,j) = (i*i-i)/2+(j);
+ − else
+ − g(i,j) = (j*j-j)/2+(i);
+ − end
+ − end
+ − end
+ −
+ − for i=1:n
+ − for j=i:n
+ − for k=1:n
+ − m_d(g(i,j),g(j,k)) = m_d(g(i,j),g(j,k)) + m_a(i,k);
+ − m_d(g(i,j),g(i,k)) = m_d(g(i,j),g(i,k)) + m_a(j,k);
+ − end
+ − end
+ − end
+ −
+ −
+ − %% setting up q from D * p = q
+ − m_q = vpa(zeros(1,g(n,n)));
+ − for i=1:n
+ − for j=i:n
+ − if i==n
+ − m_q(g(i,j)) = (E(n,n))*(E(n,n))-1;
+ − else
+ − m_q(g(i,j)) = (E(i,n)*E(j,n));
+ − end
+ − end
+ − end
+ −
+ − m_p = m_d\m_q';
+ −
+ − Cprop = vpa(zeros(n));
+ − for i=1:n
+ − for j=1:n
+ − Cprop(i,j) = m_p(g(i,j));
+ − end
+ − end
+ −
+ − Tprop = ao.mchol(Cprop);
+ − Tprop = double(Tprop);
+ − E = double(E);
+ − Tinit = double(Tinit);
+ −
+ − % Tprop = chol(Cprop,'lower');
+ −
+ − % Tprop = chol(double((Cprop)),'lower');
+ − % Tprop = mchol(Cprop);
+ −
+ − % %% writing the generator
+ − % r = randn(n,1);
+ − % y = Tinit * r;
+ − % x = zeros(Nfft,1);
+ − % for i=1:Nfft
+ − % r = randn(n,1);
+ − % y = E * y + Tprop * r;
+ − % x(i) = a*y;
+ − % end
+ −
+ − end
+ −