Mercurial > hg > ltpda
view m-toolbox/classes/@ao/ngsetup_vpa.m @ 27:29276498ebdb database-connection-manager
Remove LTPDARepositoryManager implementation
* * *
Remove GUI helper
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Mon, 05 Dec 2011 16:20:06 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
line source
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