1 | function [LD, D] = ld2ld(LD, strs, strt) |
---|
2 | % convert L'DL decomposition described by a source structure strs |
---|
3 | % to another corresponding to a target structure strt |
---|
4 | % LD = ld2ld(LD, strs, strt) |
---|
5 | % LD = ld2ld(L,D) composes the factors L and D into the matrix LD |
---|
6 | % [L,D] = ld2ld(LD) decomposes the matrix LD into the factors L and D |
---|
7 | % |
---|
8 | % LD : L'DL decomposition of a positive definite matrix |
---|
9 | % LD ... D placed instead of the unit diagonal of L |
---|
10 | % strs: source structure |
---|
11 | % strt: target structure |
---|
12 | % |
---|
13 | % Designed : M. Karny, P. Nedoma |
---|
14 | % Updated : April 2000, MK July 2004 |
---|
15 | % Project : post-ProDaCTool |
---|
16 | |
---|
17 | nPsi = length(LD); |
---|
18 | |
---|
19 | % call [L,D] = ld2ld(LD) |
---|
20 | if nargin==1, D = diag(diag(LD)); |
---|
21 | LD = LD - D + eye(nPsi); |
---|
22 | return; |
---|
23 | end |
---|
24 | |
---|
25 | % call LD = ld2ld(L,D) |
---|
26 | if nargin==2, LD = LD + strs - eye(nPsi); |
---|
27 | return; |
---|
28 | end |
---|
29 | |
---|
30 | % call LD = ld2ld(LD, strs, strt) |
---|
31 | if isempty(strt), LD=LD(1,1); return; end |
---|
32 | dim = size(strt,2); |
---|
33 | lens = size(strs,2); |
---|
34 | if lens<length(LD)-1 | lens>length(LD) |
---|
35 | error('incompatible source structure and LD'); |
---|
36 | end |
---|
37 | % add a fictive output channel (1000) |
---|
38 | if lens<length(LD) |
---|
39 | strs = [ [1000;0], strs ]; |
---|
40 | strt = [ [1000;0], strt ]; |
---|
41 | lens = lens + 1; |
---|
42 | dim = dim + 1; |
---|
43 | end |
---|
44 | for i=dim:-1:1 |
---|
45 | chn = strt(1,i); |
---|
46 | td = strt(2,i); |
---|
47 | j = find(strs(1,:)==chn & strs(2,:)==td); |
---|
48 | %chn,td,strs, strt, j, dim, |
---|
49 | if isempty (j), error('target structure item is not in source'); end |
---|
50 | LD = ldperm(LD,j); % column j to 1th position |
---|
51 | strs = [strs(:,j),strs(:,1:j-1),strs(:,j+1:lens)]; % rebuild strs |
---|
52 | %strs, pause |
---|
53 | end |
---|
54 | % cut and transform L'DL if target is shorter |
---|
55 | if lens>dim |
---|
56 | for i=2:dim |
---|
57 | LD = ldperm(LD,dim); |
---|
58 | end |
---|
59 | for i=dim+1:lens |
---|
60 | LD = ldperm(LD,lens); |
---|
61 | end |
---|
62 | LD = LD(lens-dim+1:lens, lens-dim+1:lens); |
---|
63 | LD = ldperm(LD, dim); |
---|
64 | end |
---|
65 | |
---|
66 | if nargout==2 |
---|
67 | [LD, D] = ld2ld(LD); |
---|
68 | end |
---|