Home > classes > @ao > ngsetup.m

ngsetup

PURPOSE ^

NGSETUP is called by the function LTPDA_NOISEGEN

SYNOPSIS ^

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

DESCRIPTION ^

 NGSETUP is called by the function LTPDA_NOISEGEN

 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.2 2008/08/01 13:19:42 ingo Exp $

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 % NGSETUP is called by the function LTPDA_NOISEGEN
0002 %
0003 % Inputs calculated by ...
0004 %  ... NGCONV:
0005 %            - den:   denominator coefficients
0006 %
0007 %  ... USER
0008 %            - fs:    sampling frequency given as input to LTPDA_NOISEGEN
0009 %
0010 %      Outputs:
0011 %            - Tinit: matrix to calculate initial state vector
0012 %            - Tprop: matrix to calculate propagation vector
0013 %            - E:     matrix to calculate propagation vector
0014 %
0015 % A Monsky 24-07-07
0016 %
0017 % $Id: ngsetup.m,v 1.2 2008/08/01 13:19:42 ingo Exp $
0018 
0019 function [Tinit,Tprop,E] = ngsetup(den,fs)
0020 
0021   den=den';
0022   dt = 1/fs;
0023   length(den);
0024 
0025   n = length(den)-1;
0026 
0027   %% setting up matrix Aij
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   %% Matrix exponential E
0042   E = expm(m_a*dt);
0043 
0044   %% setting up matrix Bij
0045   B = zeros(n,n);
0046   for i=1:n
0047     if rem(i,2) ~= 0
0048       j0 = (i+1)/2;
0049       s  = (-1)^(j0+1);
0050       j  = j0;
0051       for k=1:2:(n+1)
0052         B(i,j) = s*den(k);
0053         s   = -s;
0054         j = j+1;
0055       end
0056     end
0057     if rem(i,2) == 0
0058       j0 = i/2+1;
0059       s  = (-1)^j0;
0060       j  = j0;
0061       for k=2:2:(n+1)
0062         B(i,j) = s * den(k);
0063         s        = -s;
0064         j        = j+1;
0065       end
0066     end
0067   end
0068 
0069   %% solve B * m = k
0070   m_k = zeros(n,1);
0071   m_k(n) = 0.5;
0072   m_m = B\m_k;
0073 
0074   %% filling covariance matrix Cinit
0075   Cinit = zeros(n,n);
0076   for i=1:n
0077     for j=1:n
0078       if rem((i+j),2) == 0    % even
0079         Cinit(i,j) = (-1)^((i-j)/2) * m_m((i+j)/2);
0080       else
0081         Cinit(i,j) = 0;
0082       end
0083     end
0084   end
0085 
0086 
0087   %cholesky decomposition
0088 
0089   Tinit = chol(Cinit,'lower'); %lower triangular matrix
0090 
0091   %% setting up matrix D
0092   N = n*(n+1)/2;
0093 
0094   m_d = zeros(N);
0095   g = zeros(n);
0096   for i=1:n
0097     for j=1:n
0098       if i>=j
0099         g(i,j) = (i*i-i)/2+(j);
0100       else
0101         g(i,j) = (j*j-j)/2+(i);
0102       end
0103     end
0104   end
0105 
0106   for i=1:n
0107     for j=i:n
0108       for k=1:n
0109         m_d(g(i,j),g(j,k)) = m_d(g(i,j),g(j,k)) + m_a(i,k);
0110         m_d(g(i,j),g(i,k)) = m_d(g(i,j),g(i,k)) + m_a(j,k);
0111       end
0112     end
0113   end
0114 
0115 
0116   %% setting up q from D * p = q
0117   m_q = zeros(1,g(n,n));
0118   for i=1:n
0119     for j=i:n
0120       if i==n
0121         m_q(g(i,j)) = (E(n,n))*(E(n,n))-1;
0122       else
0123         m_q(g(i,j)) = (E(i,n)*E(j,n));
0124       end
0125     end
0126   end
0127 
0128   m_p = m_d\m_q';
0129 
0130   Cprop = zeros(n);
0131   for i=1:n
0132     for j=1:n
0133       Cprop(i,j) = m_p(g(i,j));
0134     end
0135   end
0136   Tprop = chol(Cprop,'lower');
0137 
0138   % %% writing the generator
0139   % r = randn(n,1);
0140   % y = Tinit * r;
0141   % x = zeros(Nfft,1);
0142   % for i=1:Nfft
0143   %     r = randn(n,1);
0144   %     y = E * y + Tprop * r;
0145   %     x(i) = a*y;
0146   % end
0147 
0148 end
0149

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