gkmu.m

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

function [A, B, INSTAB]= gkmu(Ayy, N, MAR, LAR, MMA, LMA)
% function [A, B, INSTAB]= gkmu(Ayy, N, MAR, LAR, MMA, LMA)
%   This file is part of the TFPM toolbox v0.3 (c)
%   [e-mail] and underlies the GPL.
% 
% Estimates TFARMA(MAR, LAR; MMA, LMA; N; alpha) models WITH MONIC
% TFMA PART 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_tfywu.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. 
%
% IS FUCKING UNSTABLE!!!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    =  64;
MAR  =   3;
LAR  =   3;
MMA  =   2;
LMA  = LAR;
re_im= 'r';
mo_no= 'n';
tfpm_file_gen;

%Bml(LMA+1, 1)= 6;

%-------------
beta = 1/2;
alpha= 1/2;

Aalpha= tfarma_ambi(Aml, Bml, N, alpha);

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

Ayy= Aalpha(N/2-3*LC+1:N/2+3*LC+1, N/2-MC+1:N/2+MC+1);
%max(max(abs(Ayy)))

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
[CA, CB]= tfar_est_tfywu(Ayy, N);
Cml= reshape(CA(:, :, MCmax), LL, MCmax+1);

% A pure TFAR model?
if(~MMA)
   A= Cml(LC+1-LAR:LC+1+LAR, 1:MAR+1);
   B= CB(LC+1-LMA:LC+1+LMA, 1, end);
   return;
end;

% A pure TFMA model?
%if(~MAR)
%   A= nk_to_nm(nm_to_ml(1./param_weyl(Cml, N, 1/2) ))
%end;

% The TFMA part
CCC= Cml(LC-3*LMA+1:LC+3*LMA+1, MAR-MMA+1:MAR+MMA+1);

[BML, B0L]= tfar_est_tfywu(CCC, N);



B= reshape(BML(:, :, MMA), 2*LMA+1, MMA+1)*CB(LC+1, 1, end);

% 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)*...
                       B1(LC+1+lprime, 1+mprime);
         end;
      end;
      A(LAR+1+l, 1+m)= aml;
   end;
end;

%Stable??
TDIRA= param_tdir(param_expand(A, N));
TDIRB= param_tdir(param_expand(B, N));
if(max(max(abs([TDIRA TDIRB])))>.95)
   INSTAB= 1;
end;

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

%Bml(LMA+1, 1)= 2;

%-------------
beta = 1/2;
alpha= 1/2;

KC=  2;
LC= KC*(LAR+LMA);
MC= KC*(MAR+MMA);
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);

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]= gkmu(Ayy, N, MAR, LAR, MMA, LMA)
param_norm(A, Aml, 1)
param_norm(B, Bml, 1)
param_plot('exp', N, 1, 'r', 'g', A, 1, 0, 0, 1, 0, 0, B, 1, 0, 0, 2, 0, 0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%