?? invert_waxkailath.m
字號:
function AM= invert_waxkailath(rPM, BM)% function AM= invert_waxkailath(rPM, BM)% This file is part of the TFPM toolbox v0.5 (c)% michael.jachan@tuwien.ac.at and underlies the GPL.% % Implements the Wax-Kailath Algorithm (IEEEASSP 31(5), oct83,% pp1218) for solving the Toeplitz/Block-Toeplitz system %% toeptoep(rPM)*AM= BM; P...block size, M...super size,% size(rPM)= [2M-1, 2P-1].%% ATTENTION: (JP*R).' ~= JP*(R.')!!% Dimensions:PM= size(rPM);P= (PM(2)+1)/2;M= (PM(1)+1)/2;rhs_dim= size(BM, 2);% Dimension of right-hand-side block-vector% The output tensorAM= zeros(P*M, rhs_dim, M);% Definitions:NP = zeros(P);IP = eye(P);JP = rot90(IP);JPn= JP;% Init (n=0):Alpha= toep(rPM(M, :));%Alpha_0Delta= BM(1:P, :); %Delta_0if(M==1) AM= inv(Alpha)*Delta; return;end;Beta = JP*toep(rPM(M+1, :)); %Beta_0Gamma= (toep(rPM(M-1, :))*JP).';%Gamma_0am= inv(Alpha)*Delta; %AM_0AM(:, :, 1)= [am; zeros(P*(M-1), rhs_dim)];% Init (n=1):Vn= -JP*inv(toep(rPM(M, :))) *Gamma;%V_1Wn= -JP*inv(toep(rPM(M, :))).'*Beta ;%W_1Rp= toep(rPM(M+1, :)); %R_{+1}Tm= toep(rPM(M-1, :)); %R^T_{-1}for n= 1:M-2 Alpha= Alpha - Gamma.'*inv(Alpha).'*Beta; Beta = Wn.'*JPn*Rp + JP*toep(rPM(M+1+n, :)); Gamma= Vn.'*(Tm*JPn).' + (JP*toep(rPM(M-1-n, :)).'); Delta= Vn.'*BM(1:P*n, :) + BM(P*n+1:P*(n+1), :); Vhat_new= [JPn*Vn; NP] - [Wn; IP]*inv(Alpha) *Gamma; What_new= [JPn*Wn; NP] - [Vn; IP]*inv(Alpha).'*Beta ; am= [am; zeros(P, rhs_dim)] + [Wn; IP]*inv(Alpha)*Delta; AM(:, :, n+1)= [am; zeros(P*(M-1-n), rhs_dim)]; JPn= rot90(eye((n+1)*P)); Vn= JPn*Vhat_new; Wn= JPn*What_new; Rp= [toep(rPM(M+1+n, :)); Rp]; Tm= [toep(rPM(M-1-n, :)) Tm];end;Alpha= Alpha-Gamma.'*inv(Alpha).'*Beta;Delta= Vn.'*BM(1:P*(M-1), :) + BM(P*(M-1)+1:P*M, :);am= [am; zeros(P, rhs_dim)] + [Wn; IP]*inv(Alpha)*Delta;AM(:, :, M)= am;?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -