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)%   michael.jachan@tuwien.ac.at 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 factorTDKRn0= [];% \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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%