ti_ceps2.m

来自「用于模拟时变非平稳的ARMA过程」· M 代码 · 共 96 行

function [r_hat, REGULARIZED, RMIN]= ti_ceps2(r, ABSMIN)
% function [r_hat, REGULARIZED, RMIN]= ti_ceps2(r, ABSMIN)
%   This file is part of the TFPM toolbox v0.5 (c)
%   [e-mail] and underlies the GPL.
% 
% Computes the complex cepstrum of a cplx-valued, non-causal
% signal (correlation sequence) r[m], m= -N/2:N/2-1. ABSMIN is used
% to regularize before taking the log. The switch <UNWRAP>
% determines the computation of the CTFC:
% +1: unwrapping
%  0: no unwrapping
% -1: zero phase
% -2: Pabs= real(Palpha)!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
MAR  =   0;
LAR  =   0;
MMA  =   5;
LMA  =   0;
N    =  64;
re_im= 'i';
mo_no= 'n';
tfpm_file_gen;
%-------------
alpha= 1/2;
beta = 1/2;
ABSMIN= 1e-3;

Bml= [-6 -5*j -4 -3 -2 -1];



Ralpha= tfarma_corr(Aml, Bml, N, alpha);
r= Ralpha(1, :).';

y= tfarma_gen(randn(N, 1), Aml, Bml, beta);
qYY= corr_est(y, y,   0, alpha);
r= sum(qYY)/N;


figure(97);clot(r)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Dimensions
N= length(r);

R= fft(fftshift(r));

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

Rabs= abs(R);

% Regularization: min(Rabs)->ABSMIN!
RMIN= min(Rabs);
mask= Rabs<ABSMIN;
Rabs= (1-mask).*Rabs + ABSMIN*mask;
REGULARIZED= sum(mask);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(98);
subplot(2, 1, 1);plot(Rang)
subplot(2, 1, 2);plot(Rang1)
max(abs(Rang-Rang1))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The complex log:
R_hat= log(Rabs);

r_hat= fftshift(ifft(R_hat));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r_hat1= tiarma_ceps(Aml, Bml, N);
figure(97);clot(r_hat)
figure(98);clot(r_hat1)
figure(99);clot(r_hat-r_hat1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%