Home > classes > @ao > ngsetup_vpa.m

ngsetup_vpa

PURPOSE ^

ALGONAME = mfilename;

SYNOPSIS ^

function [Tinit,Tprop,E] = ngsetup_vpa(den,fs, ndigits)

DESCRIPTION ^

 ALGONAME = mfilename;
 VERSION  = '$Id: ngsetup_vpa.m,v 1.2 2008/08/01 13:19:42 ingo Exp $';

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [Tinit,Tprop,E] = ngsetup_vpa(den,fs, ndigits)
0002 
0003   % ALGONAME = mfilename;
0004   % VERSION  = '$Id: ngsetup_vpa.m,v 1.2 2008/08/01 13:19:42 ingo Exp $';
0005   digits(ndigits)
0006   fs = sym(fs);
0007   den = sym(den,'d');
0008   den=den';
0009   dt = vpa(1/fs);
0010 
0011   n = length(den)-1;
0012   % digits(d);
0013 
0014   %% setting up matrix Aij
0015   m_a = vpa(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   %% Matrix exponential E
0028   a = m_a*dt;
0029   E = expm(a);
0030 
0031 
0032   %% setting up matrix Bij
0033   B = vpa(zeros(n,n));
0034   for i=1:n
0035     if rem(i,2) ~= 0
0036       j0 = (i+1)/2;
0037       s  = (-1)^(j0+1);
0038       j  = j0;
0039       for k=1:2:(n+1)
0040         d1 = den(k);
0041         d2 = s*d1;
0042         B(i,j) = d2;
0043         s   = -s;
0044         j = j+1;
0045       end
0046     end
0047     if rem(i,2) == 0
0048       j0 = i/2+1;
0049       s  = (-1)^j0;
0050       j  = j0;
0051       for k=2:2:(n+1)
0052         d1 = den(k);
0053         d2 = s*d1;
0054         B(i,j) = d2;
0055         s        = -s;
0056         j        = j+1;
0057       end
0058     end
0059   end
0060 
0061   %% solve B * m = k
0062   m_k = vpa(zeros(n,1));
0063   m_k(n) = 0.5;
0064   m_m = vpa(B\m_k);
0065 
0066   %% filling covariance matrix Cinit
0067   % Cinit = vpa(zeros(n,n));
0068   % for i=1:n
0069   %     for j=1:n
0070   %         if rem((i+j),2) == 0    % even
0071   %             d1 = (-1)^((i-j)/2);
0072   %             d1 = subs(d1)
0073   %             d2 = vpa('m_m((i+j)/2)');
0074   %             d2 = subs(d2)
0075   %             d3 = vpa(ctranspose(d2));
0076   %             d3 = subs(d3)
0077   %             d4 = vpa('d1 * d3');
0078   %             d4 = subs(d4)
0079   %             Cinit(i,j) = d4;
0080   %         else
0081   %             Cinit(i,j) = 0;
0082   %         end
0083   %     end
0084   % end
0085 
0086   Cinit = vpa(zeros(n,n));
0087   for i=1:n
0088     for j=1:n
0089       if rem((i+j),2) == 0    % even
0090         d1 = (-1)^((i-j)/2);
0091         d2 = (i+j)/2;
0092         Cinit(i,j) = d1 * m_m(d2);
0093         %         else
0094         %             Cinit(i,j) = 0;
0095       end
0096     end
0097   end
0098   %cholesky decomposition
0099 
0100   % Tinit = chol(double(Cinit),'lower'); %lower triangular matrix
0101   Tinit = ao.mchol(Cinit);
0102 
0103   %% setting up matrix D
0104   N = n*(n+1)/2;
0105 
0106   m_d = vpa(zeros(N));
0107   g = zeros(n);
0108   for i=1:n
0109     for j=1:n
0110       if i>=j
0111         g(i,j) = (i*i-i)/2+(j);
0112       else
0113         g(i,j) = (j*j-j)/2+(i);
0114       end
0115     end
0116   end
0117 
0118   for i=1:n
0119     for j=i:n
0120       for k=1:n
0121         m_d(g(i,j),g(j,k)) = m_d(g(i,j),g(j,k)) + m_a(i,k);
0122         m_d(g(i,j),g(i,k)) = m_d(g(i,j),g(i,k)) + m_a(j,k);
0123       end
0124     end
0125   end
0126 
0127 
0128   %% setting up q from D * p = q
0129   m_q = vpa(zeros(1,g(n,n)));
0130   for i=1:n
0131     for j=i:n
0132       if i==n
0133         m_q(g(i,j)) = (E(n,n))*(E(n,n))-1;
0134       else
0135         m_q(g(i,j)) = (E(i,n)*E(j,n));
0136       end
0137     end
0138   end
0139 
0140   m_p = m_d\m_q';
0141 
0142   Cprop = vpa(zeros(n));
0143   for i=1:n
0144     for j=1:n
0145       Cprop(i,j) = m_p(g(i,j));
0146     end
0147   end
0148 
0149   Tprop = ao.mchol(Cprop);
0150   Tprop = double(Tprop);
0151   E = double(E);
0152   Tinit = double(Tinit);
0153 
0154   % Tprop = chol(Cprop,'lower');
0155 
0156   % Tprop = chol(double((Cprop)),'lower');
0157   % Tprop = mchol(Cprop);
0158 
0159   % %% writing the generator
0160   % r = randn(n,1);
0161   % y = Tinit * r;
0162   % x = zeros(Nfft,1);
0163   % for i=1:Nfft
0164   %     r = randn(n,1);
0165   %     y = E * y + Tprop * r;
0166   %     x(i) = a*y;
0167   % end
0168 
0169 end
0170

Generated on Mon 08-Sep-2008 13:18:47 by m2html © 2003