--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/ngsetup.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,149 @@
+% NGSETUP is called by the function fromPzmodel
+%
+% Inputs calculated by ...
+%  ... NGCONV:
+%            - den:   denominator coefficients
+%
+%  ... USER
+%            - fs:    sampling frequency given as input to LTPDA_NOISEGEN
+%
+%      Outputs:
+%            - Tinit: matrix to calculate initial state vector
+%            - Tprop: matrix to calculate propagation vector
+%            - E:     matrix to calculate propagation vector
+%
+% A Monsky 24-07-07
+%
+% $Id: ngsetup.m,v 1.3 2008/10/20 08:38:29 anneke Exp $
+
+function [Tinit,Tprop,E] = ngsetup(den,fs)
+
+  den=den';
+  dt = 1/fs;
+  length(den);
+
+  n = length(den)-1;
+
+  %% setting up matrix Aij
+
+  m_a = 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
+  E = expm(m_a*dt);
+
+  %% setting up matrix Bij
+  B = 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)
+        B(i,j) = s*den(k);
+        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)
+        B(i,j) = s * den(k);
+        s        = -s;
+        j        = j+1;
+      end
+    end
+  end
+
+  %% solve B * m = k
+  m_k = zeros(n,1);
+  m_k(n) = 0.5;
+  m_m = B\m_k;
+
+  %% filling covariance matrix Cinit
+  Cinit = zeros(n,n);
+  for i=1:n
+    for j=1:n
+      if rem((i+j),2) == 0    % even
+        Cinit(i,j) = (-1)^((i-j)/2) * m_m((i+j)/2);
+      else
+        Cinit(i,j) = 0;
+      end
+    end
+  end
+
+
+  %cholesky decomposition
+
+  Tinit = chol(Cinit,'lower'); %lower triangular matrix
+
+  %% setting up matrix D
+  N = n*(n+1)/2;
+
+  m_d = 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 = 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 = zeros(n);
+  for i=1:n
+    for j=1:n
+      Cprop(i,j) = m_p(g(i,j));
+    end
+  end
+  Tprop = chol(Cprop,'lower');
+
+  % %% 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
+