--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/ngsetup_vpa.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,170 @@
+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
+