tfarma_est_gkme.m

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

function [A, B, INSTAB]= tfarma_est_gkme(Ayy, N, MAR, LAR, MMA, LMA)
% function [A, B, INSTAB]= tfarma_est_gkme(Ayy, N, MAR, LAR, MMA, LMA)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Estimates TFARMA(MAR, LAR; MMA, LMA; N; alpha) models using the
% method of Graupe/Krause/Moore (IEEEAC 20, feb75, pp104). In this
% method: MAR>MMA, LAR>LMA! The ambiguity function Ayy needs to be
% given up to MCmax, the maximum M-dimension of the long underlying
% TFAR model. Its L-dimension is LC. The estimator tfar_est_tfywe.m
% is used in here. INSTAB is zero, if the long TFAR model in stable
% (no pole magnitude > 1) otherwise the maximum value of the impulse
% response. 
% MC, LC >= ??

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

KC=  2;
LC= KC*(LAR+LMA);
MC= KC*(MAR+MMA);

Aalpha= tfarma_ambi(Aml, Bml, N, alpha);
Ayy= Aalpha(N/2-3*LC+1:N/2+3*LC+1, N/2-MC+1:N/2+MC+1);

%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, :);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

INSTAB= 0;
alpha= 1/2;

% Dimensions
MCmax= (size(Ayy, 2)-1)/2;
LC   = (size(Ayy, 1)-1)/6;
LL   = 2*LC+1;

% Estimate the long TFAR model
[Cml, CB]= tfar_est_tfywe(Ayy, N);

%param_plot('exp', N, 0, 'r', 'g', Cml, 1, 1, 2, 3, 4, 5, 1, 0, 0, 0, 0, 0, 0);
[TDIRnm, TDIRn0]= param_tdir(param_expand(Cml, N));
if(max(max(abs(TDIRnm))) > 1)
   HAR= tfarma_impr(Cml, 1, N, alpha);
   INSTAB= max(max(abs(HAR)));
end;
   
% The TFMA part
CCC= [
      Cml(LC-3*LMA+1:LC+3*LMA+1, MAR-MMA+1:MAR+MMA+1)
     ];

[B, B0L]= tfar_est_tfywe(CCC, N);

% The TFAR part
B1= [ [zeros(LC-LMA, 1+MMA); B; zeros(LC-LMA, 1+MMA)] zeros(2*LC+1, MCmax-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)*...
                       exp(2*pi*j/N*mprime*(l-lprime))*...
                       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)
INSTAB
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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  =   1;
LAR  =   1;
MMA  =   1;
LMA  =   1;
N    = 256;
re_im= 'i';
mo_no= 'm';
tfpm_file_gen;
%-------------
alpha= 1/2;
beta = 1/2;

KC=  2;
LC= KC*(LAR+LMA);
MC= KC*(MAR+MMA);

MM= 100;
Psi= tf_multiwin(N, MC, LC, 5, 1, 1);
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);
   [A, B, INSTAB]= tfarma_est_gkme(Ayy, N, MAR, LAR, MMA, LMA);
   [mm param_norm(A, Aml, 0) param_norm(B, Bml, 0) INSTAB]
   AA= AA+A;
   BB= BB+B;
end;
AA/MM
BB/MM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%