param_fact.m

用于模拟时变非平稳的ARMA过程

function [CML11, Pre, Post, NrIter]= param_fact(Cml, N, re_im, mo_no, beta)
% function [CML11, Pre, Post, NrIter]= param_fact(Cml, N, re_im, mo_no, beta)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Computes the approximate factorization of the system Cml of order
% (M, L) into (1, 1) and LTI or LFI systems. Pre/Post: GM-phd,
% p43. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    = 128;
MAR  =   2;
LAR  =   2;
MMA  =   1;
LMA  =   1;
re_im= 'r';
mo_no= 'm';
tfpm_file_gen;
%-------------
alpha= 1/2;
beta = 1/2;
Cml= param_tconv(Aml, Bml, N, beta);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Dimensions
[M, L]= param_dim(Cml);
thetaC= param_stack_ml(Cml);

CML11= zeros(3, 2, min(M, L));

M1= min(1, M);
L1= min(1, L);

M2= M-M1;
L2= L-L1;

tfeps= .125;
% Accuracy:
if(M2>3)
   tfeps= .075;
end;
if(M2>4)
   tfeps= .050;
end;
if(M2>5)
   tfeps= .025;
end;


% Step0: Init
NrIter= 0;
Pre= [];
Post= [];
%Error= [1000 100 1000 100 1000 100];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%RANDOM METHOD:
Mask1= 1./(1+(abs(-L1:L1)+1)'*(1:M1+1)).^2;
Mask2= 1./(1+(abs(-L2:L2)+1)'*(1:M2+1)).^2;

Bml1k= param_rand(M1, L1, N, 'i', 'n');
%Bml1k(L1+1, 1)= Cml(L+1, 1);
Bml2k= param_rand(M2, L2, N, 'i', 'n');
%Bml2k(L2+1, 1)= 1;
Bml1k= Bml1k.*Mask1;
Bmk2k= Bml2k.*Mask2;
if(re_im=='r')
   Bml1k= param_hermite(Bml1k);
   Bml2k= param_hermite(Bml2k);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

thetaB1= param_stack_ml(Bml1k);
thetaB2= param_stack_ml(Bml2k);

Bml= param_tconv(Bml1k, Bml2k, N, beta);
B1= param_tconv_b1(Bml1k, Bml2k, N, beta);
B2= param_tconv_b2(Bml1k, Bml2k, N, beta);

% Step1: 1st Iteration
   Bml2k= param_destack_ml(inv(B1'*B1)*B1'*thetaC, M2, L2);
   Bml2k(:, 1)= param_hermite(Bml2k(:, 1));
   B2= param_tconv_b2(Bml1k, Bml2k, N, beta);
   Bml1k= param_destack_ml(inv(B2'*B2)*B2'*thetaC, M1, L1);
   Bml1k(:, 1)= param_hermite(Bml1k(:, 1));
   B1= param_tconv_b1(Bml1k, Bml2k, N, beta);
   Bml= param_tconv(Bml1k, Bml2k, N, beta);
   Error= norm(Cml-Bml)/norm(Cml);
%   Error= [Error norm(Cml-Bml)/norm(Cml) ];Error= Error(end-5:end);
   NrIter= NrIter + 1;
% Correction
   corr= Bml1k(L1+1, 1);
   Bml1k= Bml1k/corr;
   Bml2k= Bml2k*corr;

% Step2: 2nd Iteration
   Bml2k= param_destack_ml(inv(B1'*B1)*B1'*thetaC, M2, L2);
   Bml2k(:, 1)= param_hermite(Bml2k(:, 1));
   B2= param_tconv_b2(Bml1k, Bml2k, N, beta);
   Bml1k= param_destack_ml(inv(B2'*B2)*B2'*thetaC, M1, L1);
   Bml1k(:, 1)= param_hermite(Bml1k(:, 1));
   B1= param_tconv_b1(Bml1k, Bml2k, N, beta);
   Bml= param_tconv(Bml1k, Bml2k, N, beta);
   Error= norm(Cml-Bml)/norm(Cml);
%   Error= [Error norm(Cml-Bml)/norm(Cml) ]; Error= Error(end-5:end);
   NrIter= NrIter + 1;
% Correction
   corr= Bml1k(L1+1, 1);
   Bml1k= Bml1k/corr;
   Bml2k= Bml2k*corr;

%while( (Error(end)>tfeps) & (NrIter<100) )
while(Error>tfeps)
% Iteration:
   Bml2k= param_destack_ml(inv(B1'*B1)*B1'*thetaC, M2, L2);
   Bml2k(:, 1)= param_hermite(Bml2k(:, 1));
   B2= param_tconv_b2(Bml1k, Bml2k, N, beta);
   Bml1k= param_destack_ml(inv(B2'*B2)*B2'*thetaC, M1, L1);
   Bml1k(:, 1)= param_hermite(Bml1k(:, 1));
   B1= param_tconv_b1(Bml1k, Bml2k, N, beta);
   Bml= param_tconv(Bml1k, Bml2k, N, beta);
%   Error= [Error norm(Cml-Bml)/norm(Cml)];Error= Error(end-5:end);
   Error= norm(Cml-Bml)/norm(Cml);
   [var(Error)/mean(Error) Error(end) tfeps M2 L2]

%   if( (var(Error)/mean(Error)<1e-3) & (M2>1) )%BL諨SINN!!!!!!1
%      Bml1k= param_rand(M1, L1, N, 'i', 'n').*Mask1;
%      Bml2k= param_rand(M2, L2, N, 'i', 'n').*Mask2;
%      display('REINIT!!!!!')
%      Error= [1000 100 1000 100 1000 100];
%   end;
   
   NrIter= NrIter + 1;
% Correction
   corr= Bml1k(L1+1, 1);
   Bml1k= Bml1k/corr;
   Bml2k= Bml2k*corr;
   figure(98);clf;mesh(abs(Bml1k))
   figure(99);clf;mesh(abs(Bml2k))
   drawnow
end;


CML11(:, :, 1)= Bml1k;
if( (M2==1) & (L2==1) )
   CML11(:, :, 2)= Bml2k;
   return
else
   [CML11(:, :, 2:end), xxx, yyy, nriter]= param_fact(Bml2k, N, beta);
   NrIter= NrIter + nriter
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
beta= 1/2;
Cml= param_rand(5, 5, N, 'i', 'n');
figure(1);clf;mesh(abs(Cml))

[CML11, Pre, Post, NrIter]= param_fact(Cml, N, beta)
Cml1= param_itconv(CML11, N, beta);
norm(Cml1-Cml)/norm(Cml)
figure(1);clf;mesh(abs(Cml))
figure(2);clf;mesh(abs(Cml1))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%