tvarma_est_cepsb.m

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

function [Bnm, B, REG, PMIN]= tvarma_est_cepsb(Aalpha, MMA, LMA, Aml)
% function [Bnm, B, REG, PMIN]= tvarma_est_cepsb(Aalpha, MMA, LMA, Aml)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Estimates TVARMA models from Ralpha. size(Ralpha)= [N, N]. 
% Uses the complex TV cepstral recursion in IEEESPL vol9(4), apr02,
% pp130. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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= 'n';
tfpm_file_gen;
%-------------
alpha= 1/2;
beta = 1/2;
%Psi= ones(N);
Aalpha= tfarma_ambi(Aml, Bml, N, alpha);
Bnmref=  param_expand(Bml, N);
n=0
m= 1
mprime= 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ralpha= ml_to_nm(Aalpha);

% Dimensions
N= size(Ralpha, 1);

Anm=  param_expand(Aml, N);
Bnm= [];
for n= 0:N-1
   [r_hat, REG, PMIN]= ti_ceps2(Ralpha(n+1, :).', 1e-4);

%% Recursion in m
% Initialization (m=0)
%   BB= [exp(.5*r_hat(N/2+1)) zeros(1, MMA)];
  BB= 1;
  AA= Anm(n+1, :);
% m<0
   for m= 1:MMA
      bb= m*(r_hat(N/2+m+1) + Anm(n+1, m+1));
      if(m>1)
         for mprime= 1:m-1
            bb= bb + mprime*( AA(mprime+1)*BB(m-mprime+1) - BB(mprime+1)*AA(m-mprime+1) );
         end;
         bbb= 0;
         for mprime= 1:m-1
            for mpp= 0:m-mprime
               bbb= bbb + AA(mpp+1)*BB(m-mprime-mpp+1);	 
            end;
            bb= bb + mprime*r_hat(N/2+mprime+1)*bbb;
         end;
      end;%if(m>1)
      BB(m+1)= bb/m;
   end;
   Bnm= [Bnm; BB];
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
max(max(abs(Bnm-Bnmref)))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

HMAmonic= [zeros(N, N/2) Bnm zeros(N, N/2-MMA-1)];
HAR     = [zeros(N, N/2) Anm zeros(N, N/2-MMA-1)];
Wmonic= abs(nm_to_nk(HMAmonic)./nm_to_nk(HAR)).^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);mesh(abs(Wmonic))
figure(2);mesh(abs(tfarma_wvsp(Aml, Bml, N, 1/2)))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

B0= sqrt(real(mean( (nm_to_nk(Ralpha) ./ Wmonic).' ).'));
%figure(1);clf;plot(real(B0));
%figure(2);clf;plot(real(param_expand(Bml, N)))
Bnm= Bnm.*(B0*ones(1, MMA+1));
B= param_compress(Bnm, LMA);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
max(max(abs(Bnm-Bnmref)))
max(max(abs(Bml-B)))
param_norm(param_expand(Bml, N), Bnm, 1)
param_norm(Bml, B, 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%