tfarma_est_cepsb.m

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

function [B, REG, PMIN]= tfarma_est_cepsb(Aalpha, Aml, MMA, LMA)
% function [B, REG, PMIN]= tfarma_est_cepsb(Aalpha, Aml, MMA, LMA)
%   This file is part of the TFPM toolbox v0.5 (c)
%   [e-mail] and underlies the GPL.
% 
% Estimates TFMA(MMA; LMA; N; 1/2) models from Aalpha.
% size(Aalpha)= [N, N], Aalpha(N/2+1, N/2+1)= A_y[0, 0]. 
% Uses a cepstral recursion based on IEEESPL vol9(4), apr02, pp130. 

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

% Dimensions
N= size(Aalpha, 1);
[MAR, LAR]= param_dim(Aml);

[Aalpha_hat, REG, PMIN]= tf_ceps2(Aalpha, 1e-4);
Aalpha_hat= Aalpha_hat/N;%<===========================

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(98);cesh(Aalpha_hat)
REG
PMIN
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Recursion in m
% Initialization (m=0)
A0= [Aalpha_hat(N/2+1:N, N/2+1); Aalpha_hat(1:N/2, N/2+1)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(97);clot(A0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rx_hat= real(ifft(A0))*N;
b0n= real(exp(rx_hat/2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(97);clf;hold on
   plot(b0n, 'r')
   plot(real(Bmn(:, 1)), 'g')
hold off;drawnow;

RX_HAT= 2*log(real(Bmn(:, 1)));
figure(98);clf;hold on
   plot(rx_hat, 'r')
   plot(RX_HAT, 'g')
hold off;drawnow;

figure(97);clf;
   subplot(2, 1, 1);hold on;
      plot(real(A0), 'r');
      plot(real(fft(RX_HAT/N)), 'g');
   subplot(2, 1, 2);hold on;
      plot(imag(A0), 'r');
      plot(imag(fft(RX_HAT/N)), 'g');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

b0l= fft(b0n)/N;
B= [[b0l(N-LMA+1:N); b0l(1:LMA+1)] zeros(2*LMA+1, MMA)];

for m= 1:MMA
   for l= -LMA:LMA
      bb= 0;
      for mprime= 0:m-1
         for lprime= -LAR:LAR
            d= 0;
	    for mpp= 0:mprime
               for lpp= -LAR:LAR
                  d= d+param_get(Aml, mpp, lpp)*param_get(B, mprime-mpp, lprime-lpp);
	       end;
	    end;
            bb= bb + (m-mprime)/m*d*Aalpha_hat(N/2+1+l-lprime, N/2+1+m-mprime) ;
	 end;
      end;
      for mprime= 1:m-1
         for lprime= -LAR:LAR
            bb= bb + (m-mprime)/m*(param_get(Aml, m-mprime, l-lprime)*param_get(B, mprime, lprime) - ...
                param_get(Aml, mprime, lprime)*param_get(B, m-mprime, l-lprime));
	 end;
      end;
      for lprime= -LAR:LAR
         bb= bb+param_get(Aml, m, lprime)*param_get(B, 0, l-lprime);
      end;
      B(LMA+1+l, m+1)= bb;
   end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
param_norm(B, 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Recursion in l

% Initialization (l=0)
A0= [Aalpha_hat(N/2+1, N/2+1:N) Aalpha_hat(N/2+1, 1:N/2)];

rx_hat= real(ifft(A0))*N;
bm0= real(exp(rx_hat/2));

bk0= fft(bm0)/N;
B1= [zeros(LMA, MMA+1); bk0(1:MMA+1); zeros(LMA, MMA+1)];
for m= 0:MMA
   for l= -LMA:-1
      bb= 0;
      for mprime= 0:m-1
         for lprime= -LMA:LMA
            bb= bb + (l-lprime)*param_get(B1, mprime, lprime)*...
		Aalpha_hat(N/2+1+l-lprime, N/2+1+m-mprime) ;
	 end;
      end;
      B1(LMA+1+l, m+1)= bb/l;
   end;
   for l= 1:LMA
      bb= 0;
      for mprime= 0:m-1
         for lprime= -LMA:LMA
            bb= bb + (l-lprime)*param_get(B1, mprime, lprime)*...
		Aalpha_hat(N/2+1+l-lprime, N/2+1+m-mprime) ;
	 end;

      end;
      B1(LMA+1+l, m+1)= bb/l;
   end;
end;
param_norm(B, Bml, 1)
param_norm(B1, Bml, 1)
param_norm(B1, B, 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%