tfarma_sid_gkmu.m

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

function [A, B]= tfarma_sid_gkmu(Salpha, MAR, LAR, MMA, LMA, N)
% function [A, B]= tfarma_sid_gkmu(Salpha, MAR, LAR, MMA, LMA, N)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
MAR  =   3;
LAR  =   2;
MMA  =   2;
LMA  = LAR;
N    = 256;
re_im= 'i';
mo_no= 'n';
tfpm_file_gen;
%-------------
beta = 1/2;
alpha= 1/2;

Salpha= tfarma_sprf(Aml, Bml, N, alpha);
figure(1);tf_show(Salpha)

%N= 16;
%Salpha= (-N/2:N/2-1).'*ones(1, N) + ones(N, 1)*(j*(-N/2:N/2-1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

alpha= 1/2;

% Dimensions
N= size(Salpha, 1);
LL   = 2*LMA+1;

% Fill the Toeplitz/block-Toeplitz generator rPM
rPM= [];
for mu= MAR-1:-1:-MAR+1
   rPM= [rPM; fliplr(Salpha(N/2-2*LAR+1:N/2+2*LAR+1, N/2-2+MMA+MAR+mu).')];
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
toeptoep(rPM)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

RHS= [];
for m= 1:MAR
   RHS= [RHS; -Salpha(N/2-LAR+1:N/2+LAR+1, N/2+m+1)];
end;

% Solve the system
AAstack= invert_waxkailath(rPM, RHS);

A= [[zeros(LAR, 1); 1; zeros(LAR, 1)] reshape(AAstack(:, :, end), LL, MAR)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reshape to BB
BB= zeros(2*LMA+1, 1+MMA, MMA);
for m= 1:MMA
   BB(:, :, m)= [[zeros(LMA, 1); param_get(D0l, 0, 0); zeros(LMA, 1)] reshape(BBstack(1:LL*m, 1, m), LL, m) zeros(LL, MMA-m)];
end;
param_norm(BB(:, :, end), Bml, 1)
param_norm(B, Bml, 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The TFAR part
B1= [ [zeros(LC-LMA, 1+MMA); B/param_get(D0l, 0, 0); zeros(LC-LMA, 1+MMA)] zeros(2*LC+1, MC-MMA) ];
A= zeros(2*LAR+1, 1+MAR);
A(LAR+1, 1)= 1;
for m= 1:MAR
   for l= -LAR:LAR
      aml= Cml(LC+1+l, 1+m);
      for mprime= 1:m
         for lprime= -LAR:LAR
            aml= aml + Cml(LC+1+l-lprime, 1+m-mprime)*...
                       B1(LC+1+lprime, 1+mprime);
         end;
      end;
      A(LAR+1+l, 1+m)= aml;
   end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
param_norm(A, Aml, 0)
param_norm(B, Bml, 0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
MAR  =   3;
LAR  =   2;
MMA  =   2;
LMA  =   2;
N    = 256;
re_im= 'i';
mo_no= 'm';
tfpm_file_gen;
%-------------
alpha= 1/2;
beta = 1/2;
mm   =   1;

KC=  2;
LC= KC*(LAR);
MC= KC*(MAR);

MM= 1000;
%Psi= tf_multiwin(N, MC, LC, 5, 1, 1);
Psi= ones(N);
psi= Psi(N/2+1-3*LC:N/2+1+3*LC, N/2+1-MC:N/2+1+MC);

AA= zeros(size(Aml));
BB= zeros(size(Bml));
for mm= 1:MM
   y= tfarma_gen(randn(N, 1), Aml, Bml, beta);
   qyy= corr_est(y, y, MC, 1/2);
   ayy= fft(qyy);
   ayy= [ayy(N/2+1:N, :); ayy(1:N/2, :)];
   Ayy= ayy(N/2+1-3*LC:N/2+1+3*LC, :).*conj(psi);
% Estimate the long TFAR model
   [CA, CB]= tfar_est_tfywu(Ayy, N);
   Cml= CA(:, :, end);
   D0l= CB(:, :, end);
   [A, B]= tfarma_est_gkmu(Cml, D0l, MAR, LAR, MMA, LMA, N);
   [mm param_norm(A, Aml, 0) param_norm(B, Bml, 0)]
   AA= AA+A;
   BB= BB+B;
end;
AA/MM
Aml
BB/MM
Bml
[param_norm(AA/MM, Aml, 0) param_norm(BB/MM, Bml, 0)]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%