param_tdfkm.m

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

function [TDFKnM, TDFKn0, Enm]= param_tdfkM(Cnm, Initval)
% function [TDFKnM, TDFKn0, Enm]= param_tdfkM(Cnm, Initval)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Computes the right outer Kamen root [Kamen, LINALG, vol98, pp263,
% 1988] of a cyclic TV polynomial given through Cnm. If no initial
% value is  given (Initval==[]), the TDFR with minimum magnitude
% and phase is taken?. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
MAR  =   0;
LAR  =   0;
MMA  =   3;
LMA  =   3;
N    =  16;
re_im= 'r';
mo_no= 'n';
tfpm_file_gen;
%-------------
Cml= Bml;
Cnm= param_ml_to_nm(Cml, N);
%Cnm= ones(N, 1)*BLTI;
Initval= [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

% TV Amplification factor
TDFKn0= [];% \parc{0}[0:N-1]
for n= 0:N-1
   cn0= Cnm(n+1, 1);% \parc{0}[n]
   Cnm(n+1, :)= Cnm(n+1, :)/cn0;
   TDFKn0= [TDFKn0; cn0];
end;

if(nargin==1)
   Initval= [];
end;

if(~M)%% Multiplier has no roots!
   TDFKnm= zeros(N, M);% \PORO{0:0}[0:N-1]
   Enm= Cnm;
else
   if(M==1)%% Order-1 system has very simple roots!
      TDFKnm= -Cnm(:, 2)./Cnm(:, 1);% \PORO{1:1}[0:N-1]= -\parc{1}[0:N-1]
   else
      if(isempty(Initval))
% Systen is TI prior to n= 0, takes min-root in abs and angle!!
%         Initval= sort(roots(Cnm(1, :)));Initval= Initval(end);%*ones(M-1, 1);
         Initval= j;
      end;

% The recursion
      if(M>2)
         TDFKnM= Initval.*ones(M-1, 1);
         Enm= [];
         for n= 0:N-1
            eMM= 0;
            for MM= M-1:-1:1
               eMM= [eMM; double(( eMM(M-MM)-Cnm(n+1, 2+MM) )/TDFKnM(n+M-MM))];
            end;
            Enm= [Enm; [1; flipud(eMM(2:end))].'];
            TDFKnM= [TDFKnM; eMM(M)-Cnm(n+1, 2)];
         end;
         TDFKnM= TDFKnM(M:end);
%         TDFKnm= [param_tdfk(Enm, Initval) TDFKnM];
      else%% The final step, M==2
         TDFKnM= [];
	 Enm= [];
         for n= 0:N-1
            p1= Cnm(n+1, 3)/Initval;
            p2= -p1-Cnm(n+1, 2);
            Initval= p2;
            TDFKnM= [TDFKnM; p2];
	    Enm= [Enm; -p1*Cnm(n+1, 1)];
         end;
      end;%if(M>2)%% The recursion
   end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Cml= param_rand(3, 1, N, 'r', 'n');
Cnm= param_ml_to_nm(Cml, N);

Initval= [];M= MMA;
[TDFKnM, Enm]= param_tdfkM(Cnm, []);
Initval= TDFKnM(N-M+1:N)
TDFKnM

[TDFKnm, TDFKn0]= param_tdfk(Cnm, []);

figure(1);clf;tfpm_plot_roots(TDFKnm, max(1, max(max(abs(real(TDFKnm))))), ...
			              max(1, max(max(abs(imag(TDFKnm))))), 'x')
figure(2);clf;plot(real(TDFKnm8))
mesh(abs(TDFKnm2(N+1:end, :) - TDFKnm8(7*N+1:end, :)))

%[TDFKnm2, TDFKn02]= param_tdfk([Cnm;Cnm]);
[TDFKnm8, TDFKn08]= param_tdfk([Cnm;Cnm;Cnm;Cnm;Cnm;Cnm;Cnm;Cnm]);


[TDFKnm, TDFKn0]= param_tdfk(Cnm);
clf
tfpm_plot_roots(TDFKnm, max(max(abs(real(TDFKnm)))), max(max(abs(imag(TDFKnm)))), 'o')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%