param_kpoly.m

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

function Cmn= param_kpoly(CKmn, CK0, initval)
% function Cmn= param_kpoly(CKmn, CK0, initval)
%   This file is part of the TFPM toolbox v0.9 (c)
%   [e-mail] and underlies the GPL.
% 
% Computes model parameters from Kamen roots given in CKmn. Initval
% is always 's'. CHANGE!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;%tfpm;
MAR  =   7;
LAR  =   5;
MMA  =   5;
LMA  =   3;
N    = 128;
re_im= 'r';
mo_no= 'n';
tfpm_file_gen;
%-------------
C= Bmn;
[CKmn, CK0]= param_kroots(C, 's');
initval= 's';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Dimensions:
NM= size(CKmn);
N= NM(1);
M= NM(2);

if(nargin==1)
   CK0= ones(N, 1);
   initval= 'c';
end;

if(nargin==2)
   initval= 'c';
end;

if(~M)%% Multiplier has no parameters!
   Cmn= ones(N, 1);
else
   if(M==1)%% Order-1 system has very simple parameters!
      Cmn= [ones(N, 1) -CKmn];
   else

% Set initval value (STATIONARY!!):CHANGE!!!!!!!!!!!
%      if(initval=='s')
%         Initval= roots(Cmn(1, :));
%%         Initval= Initval(1);
%         Initval= Initval(M)*ones(M-1, 1);
%      else

% Set initval value (CYCLIC!!):CHANGE!!!!!!!!!!!
%         if(initval=='c')
%            Initval= [];
%            for mm= 1:M-1
%               Init= roots(Cmn(N+1-mm, :));
%%               Init= Init(1);
%               Init= Init(M);
%               Initval= [Init; Initval];
%            end;
%         else
%            initval= initval(:);
%            Initval= initval(1:M-1);
%         end;%if(initval=='c')
%      end;%if(initval=='s')

      if(M>2)%% The recursion
         Emn= param_kpoly(CKmn(:, 2:end), ones(N, 1), initval);
	 Cmn= [Emn zeros(N, 1)] - ...
	      [ zeros(N, 1) [ Emn(1, :); Emn(1:N-1, :)] ...
	      ].*(CKmn(:, 1)*ones(1, M+1));
      else%% The final step
         Cmn= [];
         for n= 0:N-1
            c1= -CKmn(n+1, 1)-CKmn(n+1, 2);
	    c2= CKmn(n+1, 1)*CKmn(max(n, 1), 2);
	    Cmn= [Cmn; [1 c1 c2]];
	 end;
      end;%if(M>2)%% The recursion
   end;
end;

for n= 0:N-1
   Cmn(n+1, :)= Cmn(n+1, :)*CK0(n+1);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
max(max(abs(Cmn-C)))
figure(99);mesh(abs(Cmn-C))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%