Home > m > noisegenerator > ngsetup.m

ngsetup

PURPOSE ^

ALGONAME = mfilename;

SYNOPSIS ^

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

DESCRIPTION ^

 ALGONAME = mfilename;
 VERSION  = '$Id: ngsetup.m,v 1.1 2007/07/31 15:49:12 anneke Exp $';

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [Tinit,Tprop,E] = ngsetup(den,fs)
0002 
0003 % ALGONAME = mfilename;
0004 % VERSION  = '$Id: ngsetup.m,v 1.1 2007/07/31 15:49:12 anneke Exp $';
0005 
0006 
0007 den=den';
0008 dt = 1/fs;
0009 length(den);
0010 
0011 n = length(den)-1;
0012 
0013 %% setting up matrix Aij
0014 
0015 m_a = 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 %% Matrix exponential E
0027 E = expm(m_a*dt);
0028 
0029 %% setting up matrix Bij
0030 B = zeros(n,n);
0031 for i=1:n
0032     if rem(i,2) ~= 0
0033         j0 = (i+1)/2;
0034         s  = (-1)^(j0+1);
0035         j  = j0;
0036         for k=1:2:(n+1)
0037             B(i,j) = s*den(k);            
0038             s   = -s;
0039             j = j+1;
0040         end 
0041     end 
0042     if rem(i,2) == 0
0043         j0 = i/2+1;
0044         s  = (-1)^j0;
0045         j  = j0;
0046         for k=2:2:(n+1)
0047             B(i,j) = s * den(k);            
0048             s        = -s;
0049             j        = j+1;
0050         end 
0051     end 
0052 end
0053 
0054 %% solve B * m = k
0055 m_k = zeros(n,1);
0056 m_k(n) = 0.5;
0057 m_m = B\m_k;
0058 
0059 %% filling covariance matrix Cinit
0060 Cinit = zeros(n,n);
0061 for i=1:n
0062     for j=1:n
0063         if rem((i+j),2) == 0    % even
0064             Cinit(i,j) = (-1)^((i-j)/2) * m_m((i+j)/2);
0065         else
0066             Cinit(i,j) = 0;
0067         end
0068     end
0069 end
0070 
0071 
0072 %cholesky decomposition
0073 
0074 Tinit = chol(Cinit,'lower'); %lower triangular matrix
0075 
0076 %% setting up matrix D
0077 N = n*(n+1)/2;
0078 
0079 m_d = zeros(N);
0080 g = zeros(n);
0081 for i=1:n
0082     for j=1:n
0083         if i>=j
0084             g(i,j) = (i*i-i)/2+(j);
0085         else
0086             g(i,j) = (j*j-j)/2+(i);
0087         end
0088     end
0089 end
0090 
0091 for i=1:n
0092     for j=i:n
0093         for k=1:n
0094             m_d(g(i,j),g(j,k)) = m_d(g(i,j),g(j,k)) + m_a(i,k);
0095             m_d(g(i,j),g(i,k)) = m_d(g(i,j),g(i,k)) + m_a(j,k);    
0096         end
0097     end
0098 end
0099 
0100 
0101 %% setting up q from D * p = q
0102 m_q = zeros(1,g(n,n));
0103 for i=1:n
0104     for j=i:n
0105         if i==n
0106             m_q(g(i,j)) = (E(n,n))*(E(n,n))-1;
0107         else
0108             m_q(g(i,j)) = (E(i,n)*E(j,n));
0109         end
0110     end
0111 end
0112 
0113 m_p = m_d\m_q';
0114 
0115 Cprop = zeros(n);
0116 for i=1:n
0117     for j=1:n        
0118             Cprop(i,j) = m_p(g(i,j));           
0119     end
0120 end
0121 Tprop = chol(Cprop,'lower');
0122 
0123 % %% writing the generator
0124 % r = randn(n,1);
0125 % y = Tinit * r;
0126 % x = zeros(Nfft,1);
0127 % for i=1:Nfft
0128 %     r = randn(n,1);
0129 %     y = E * y + Tprop * r;
0130 %     x(i) = a*y;
0131 % end

Generated on Mon 03-Sep-2007 12:12:34 by m2html © 2003