param_tdkr.m

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

function [TDKRnm, TDKRn0]= param_tdkr(Cnm, Initval)
% function [TDKRnm, TDKRn0]= param_tdkr(Cnm, Initval)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Computes the Kamen roots [Kamen, LINALG, vol98, pp263, 1988] of
% a TV polynomial given through Cnm. If no initial value is 
% given (Initval==[]), the TDIR 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  =   2;
N    = 256;
re_im= 'r';
mo_no= 'm';
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
TDKRn0= [];% \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;
   TDKRn0= [TDKRn0; cn0];
end;

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

if(~M)%% Multiplier has no roots!
   TDKRnm= zeros(N, M);% \PORO{0:0}[0:N-1]
else
   if(M==1)%% Order-1 system has very simple roots!
      TDKRnm= -Cnm(:, 2);% \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)
         PM= 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) )/PM(n+M-MM))];
            end;
            Enm= [Enm; [1; flipud(eMM(2:end))].'];
            PM= [PM; eMM(M)-Cnm(n+1, 2)];
         end;
         PM= PM(M:end);
         TDKRnm= [param_tdkr(Enm, Initval) PM];
      else%% The final step, M==2
         TDKRnm= [];
         for n= 0:N-1
            p1= Cnm(n+1, 3)/Initval;
            p2= -p1-Cnm(n+1, 2);
            Initval= p2;
            TDKRnm= [TDKRnm; [p1 p2]];
         end;
      end;%if(M>2)%% The recursion
   end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M= 3;
L= 2;
Cml= param_rand(M, L, N, 'i', 'n');
Cnm= param_ml_to_nm(Cml, N);

Initval= Cnm(N-M+2:N, end)
[TDKRnm, TDKRn0]= param_tdkr(Cnm, Initval);
figure(1);clf;tfpm_plot_roots(TDKRnm, max(1, max(max(abs(real(TDKRnm))))), ...
			              max(1, max(max(abs(imag(TDKRnm))))), 'x')
figure(2);clf;plot(real(TDKRnm8))
mesh(abs(TDKRnm2(N+1:end, :) - TDKRnm8(7*N+1:end, :)))

%[TDKRnm2, TDKRn02]= param_tdfk([Cnm;Cnm]);



M= 7;
L= 5;
Cml= param_rand(M, L, N, 'i', 'n');
Cnm= param_ml_to_nm(Cml, N);
[TDKRnm81, TDKRn08]= param_tdkr([Cnm;Cnm;Cnm], j*randn);
[TDKRnm82, TDKRn08]= param_tdkr([Cnm;Cnm;Cnm], j*randn);
figure(1);clf;plot(abs(TDKRnm81))
figure(2);clf;plot(abs(TDKRnm82))
figure(3);clf;plot(abs(TDKRnm82(2*N+1:end, :)-TDKRnm81(2*N+1:end, :)))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%