www.pudn.com > tfarma10.rar > invert_akaike.m, change:2003-06-19,size:2508b

```function AM= invert_akaike(rPM, BM)
% function AM= invert_akaike(rPM, BM)
%   This file is part of the TFPM toolbox v0.5 (c)
%   michael.jachan@tuwien.ac.at and underlies the GPL.
%
% Implements the Akaike algorithm (SIAMAM 24(2), mar73,
% pp234) for solving the block-Toeplitz system
%
% btoep(rPM)*AP = BP; M...block size, P...super size,
% size(rPM)= [M, M(2P-1)].

%% Dimensions:
PP= size(rPM);
M= PP(1);
P= (PP(2)/M+1)/2;

% Permute rPM and transpose BM!!
rrPM= [];
for p= 1:2*P-1
rrPM= [rPM(:, (p-1)*M+1:p*M).' rrPM];
end;
rPM= rrPM;
BM= BM.';

rhs_dim= size(BM, 1);% Dimension of right-hand-side block-vector

%% The output tensor
AM= zeros(P*M, rhs_dim, P);

%% Here is the algorithm for the system [A1..AP] T= [B1..BP]!!
% Init (p= 0):
qinv= rPM(:, (P-1)*M+1:P*M);
sinv= qinv;
am  = BM(:, 1:M)*inv(qinv);
AM(:, :, 1)= [am.'; zeros(M*(P-1), rhs_dim)];

if(P==1)
return;
end;

Xp= -rPM(:, (P-2)*M+1:(P-1)*M) * inv(qinv);%X1
Yp= -rPM(:, (P  )*M+1:(P+1)*M) * inv(qinv);%Y1

ap= rPM(:, (P-2)*M+1:(P-1)*M);             %a1
rp= rPM(:, (P  )*M+1:(P+1)*M);             %r1

qinv_new= rPM(:, (P-1)*M+1:P*M) - rPM(:, (P-2)*M+1:(P-1)*M)*inv(qinv)*rPM(:, (P  )*M+1:(P+1)*M);
sinv_new= rPM(:, (P-1)*M+1:P*M) - rPM(:, (P  )*M+1:(P+1)*M)*inv(qinv)*rPM(:, (P-2)*M+1:(P-1)*M);

for p= 1:P-2
Ep= block_exchange(M, p);
qinv= qinv_new;
sinv= sinv_new;

App= zeros(rhs_dim, M);
for m= 1:p
App= App + am(:, (m-1)*M+1:m*M)*rPM(:, (P+p-m)*M+1:(P+p+1-m)*M);
end;
App= (BM(:, p*M+1:(p+1)*M)-App)*inv(qinv);
am_new= [];
for m= 1:p
am_new= [am_new am(:, (m-1)*M+1:m*M) + App*Xp(:, (m-1)*M+1:m*M)];
end;
am= [am_new App];
AM(:, :, p+1)= [am.'; zeros(M*(P-1-p), rhs_dim)];

Xp_new= [zeros(M) Xp] - (rPM(:, (P-2-p)*M+1:(P-1-p)*M) + Xp*   ap)*inv(sinv)*[eye(M) Yp];
Yp_new= [Yp zeros(M)] - (rPM(:, (P  +p)*M+1:(P+1+p)*M) + Yp*Ep*rp)*inv(qinv)*[Xp eye(M)];

ap= [ap; rPM(:, (P-2-p)*M+1:(P-1-p)*M)];
rp= [rp; rPM(:, (P  +p)*M+1:(P+1+p)*M)];

qinv_new= qinv - Xp_new(:, 1:M)          *Yp_new(:, p*M+1:(p+1)*M) * qinv;
sinv_new= sinv - Yp_new(:, p*M+1:(p+1)*M)*Xp_new(:, 1:M)           * sinv;

Xp= Xp_new;
Yp= Yp_new;
end;

p= P-1;
qinv= qinv_new;
sinv= sinv_new;
App= zeros(rhs_dim, M);
for m= 1:p
App= App + am(:, (m-1)*M+1:m*M)*rPM(:, (P+p-m)*M+1:(P+p+1-m)*M);
end;
App= (BM(:, p*M+1:(p+1)*M)-App)*inv(qinv);
am_new= [];
for m= 1:p
am_new= [am_new am(:, (m-1)*M+1:m*M) + App*Xp(:, (m-1)*M+1:m*M)];
end;
am= [am_new App];
AM(:, :, P)= am.';
```