tfma_est_mepsb.m

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

function [B, sigma, REG, PMIN]= tfma_est_mepsb(Aalpha, MMA, LMA)
% function [B, sigma, REG, PMIN]= tfma_est_mepsb(Aalpha, MMA, LMA)
%   This file is part of the TFPM toolbox v0.5 (c)
%   [e-mail] and underlies the GPL.
% 
% Estimates monic TFMA(MMA; LMA; N; 1/2) models and B0 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;
MAR  =   0;
LAR  =   0;
MMA  =   3;
LMA  =   2;
N    =  64;
re_im= 'r';
mo_no= 'm';
tfpm_file_gen;
%-------------
alpha= 1/2;
beta = 1/2;
%Bml(1:LMA, 1)= zeros(LMA, 1);
%Bml(LMA+2:end, 1)= zeros(LMA, 1);
sigma0= 1.4;
Aalpha= tfarma_ambi(Aml, sigma0*Bml, N, alpha);

Psi= tf_multiwin(N, MMA, 2*LMA, 5, min(3, N/32), 2);

y= tfarma_gen(randn(N, 1), Aml, Bml, beta);
aalpha= ambi_est_cyc(y);
Aalpha= aalpha.*conj(Psi);

figure(97);cesh(Aalpha)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Dimensions
N= size(Aalpha, 1);

[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 (B0)
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;
sigma= b0l(1);
B= [[zeros(LMA, 1); 1; zeros(LMA, 1)] zeros(2*LMA+1, MMA)];

% m<0
for m= 1:MMA
   for l= -LMA:LMA
      bb= 0;
      for mprime= 0:m-1
         for lprime= -LMA:LMA
            bb= bb + (m-mprime)*param_get(B, mprime, lprime)*Aalpha_hat(N/2+1+l-lprime, N/2+1+m-mprime) ;
	 end;
      end;
      B(LMA+1+l, m+1)= bb/m;
   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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%