Home > m > noisegenerator > 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/01/23 17:23:56 anneke Exp $

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Mon 31-Mar-2008 13:54:54 by m2html © 2003