tfar_est_tfywu.m

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

function [AA, BB]= tfar_est_tfywu(Ayy, N)
% function [AA, BB]= tfar_est_tfywu(Ayy, N)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Estimates TFAR(MAR, LAR; N; 1/2) models from Ayy. 
% size(Ayy)= [6*LAR+1, 2*MAR+1], Ayy(3*LAR+1, MAR+1)= A_y[0, 0]. 
% Uses the underspread TFYW estimator order-recursively in m. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    = 256;
MAR  =   3;
LAR  =   2;
MMA  =   0;
LMA  = LAR;
re_im= 'i';
mo_no= 'n';
tfpm_file_gen;
%-------------
alpha= 1/2;
beta= 1/2;
%Aalpha= tfarma_ambi(Aml, Bml, N, alpha);
%Ayy= Aalpha(N/2-3*LAR+1:N/2+3*LAR+1, N/2-MAR+1:N/2+MAR+1);

x= tfarma_gen(randn(N, 1), Aml, Bml, beta);
y= tfarma_gen(randn(N, 1), Aml, Bml, beta);
qyy= corr_est(x, x, MAR, 1/2);
ayy= fft(qyy);
ayy= [ayy(N/2+1:N, :); ayy(1:N/2, :)];
ayy= ayy(N/2+1-3*LAR:N/2+1+3*LAR, :);
Ayy= ayy;
%HAR= tvarma_impr(Amn, Bmn);
%AAR= nm_to_ml(ker_to_lag(HAR*HAR', 0, alpha));
%Ayy= AAR   (N/2-3*LAR+1:N/2+3*LAR+1, N/2-MAR+1:N/2+MAR+1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

% The output
AA= zeros(LL, MARmax+1, MARmax);
BB= zeros(LL,        1, MARmax);

% Fill the Toeplitz/block-Toeplitz generator rPM and the RHS BM
rPM= [];
for mu= MARmax-1:-1:-MARmax+1%This is the right one;)
   rPM= [rPM; fliplr(Ayy(3*LAR+2-LL:3*LAR+LL, MARmax+mu+1).')];
end;

BM = [];
for m= 1:MARmax
   BM= [BM; -Ayy(2*LAR+1:4*LAR+1, MARmax+m+1)];
end;

% Solve the system
AM= invert_waxkailath(rPM, BM);

% Reshape AM to A and Compute B
for k= 1:MARmax
   Aml_e= [[zeros(LAR, 1); 1; zeros(LAR, 1)] reshape(AM(1:LL*k, 1, k), LL, k)];
   AA(:, :, k)= [Aml_e zeros(LL, MARmax-k)];
   B0l= tfar_est_b0l(Ayy(:, MARmax+1-k:MARmax+1+k), Aml_e, N);
   BB(:, :, k)= B0l;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
param_norm(AA(:, :, MAR), Aml, 0)
param_norm(BB(:, :, MAR), Bml, 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Inverse filtering:
e= tfarma_inv(y, AA(:, :, end), BB(:, :, end));
var(e)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%