tfar_est_cepsb.m

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

function [A, B, REG, PMIN]= tfar_est_cepsb(Aalpha, MAR, LAR)
% function [A, B, REG, PMIN]= tfar_est_cepsb(Aalpha, MAR, LAR)
%   This file is part of the TFPM toolbox v0.5 (c)
%   [e-mail] and underlies the GPL.
% 
% Estimates TFAR(MAR; LAR; 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;
MAR  =   3;
LAR  =   0;
MMA  =   0;
LMA  = LAR;
N    = 256;
re_im= 'i';
mo_no= 'n';
tfpm_file_gen;
%-------------
alpha= 1/2;
beta = 1/2;

Aalpha= tfarma_ambi(Aml, Bml, N, alpha);
%Halpha= tfarma_impr(Aml, Bml, N, alpha);

Psi= tf_multiwin(N, 3*MMA, 4*LMA, 0, 1, 1);

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 (m=0)
B0= [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(B0))*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-LAR+1:N); b0l(1:LAR+1)];
A= [zeros(LAR, 1); 1; zeros(LAR, 1)];

% m<0
for m= 1:MAR
   for l= -LAR:LAR
      aa= 0;
      for mprime= 0:m-1
         for lprime= -LAR:LAR
            aa= aa + (m-mprime)*param_get(A, mprime, lprime)*...
		Aalpha_hat(N/2+1+l-lprime, N/2+1+m-mprime) ;
	 end;
      end;
      A(LAR+1+l, m+1)= -aa/m;
   end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
param_norm(A, Aml, 1)
param_norm(B, Bml, 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%