tfpm_generate_stable.m

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

function Cml= tfpm_generate_stable(M, L, N, re_im, mo_no)
% function Cml= tfpm_generate_stable(M, L, N, re_im, mo_no)
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
M    =   2;
L    =   4;
N    = 256;
re_im= 'i';
mo_no= 'n';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(~M)
   C= (randn(L+1, 1) + j*randn(L+1, 1))/5;
   Cml= conv(C, flipud(C'.'));
   if(mo_no=='m')
      Cml(:, 1)= [zeros(L, 1); 1; zeros(L, 1)];
   end;
   return;
end;

win= ([1:L L+1 L:-1:1].'*(M+1:-1:1)/(M+1)/(L+1)).^(1/4);

Cml= win.*(randn(2*L+1, M+1) + j*randn(2*L+1, M+1))/sqrt(1+M*L);
plot(real(param_expand(Cml(:, 1), N)))
if(re_im=='r')
   Cml= param_hermite(Cml);
end;
if(mo_no=='m')
   Cml(:, 1)= [zeros(L, 1); 1; zeros(L, 1)];
else
   C= (randn(L+1, 1) + j*randn(L+1, 1))/10;
   C= conv(C, flipud(C'.'));
   C= C./C(L+1);
   Cml(:, 1)= C;
end;
maxmax= max(max(abs(param_tdir(Cml, N))))

while( (maxmax <1.1) || (maxmax>3.0) )
   Cml= win.*(randn(2*L+1, M+1) + j*randn(2*L+1, M+1))/sqrt(1+M*L);
   Cml(L+1, 1)= 2*mean(abs(Cml(:, 1)));
   if(re_im=='r')
      Cml= param_hermite(Cml);
   end;
   if(mo_no=='m')
      Cml(:, 1)= [zeros(L, 1); 1; zeros(L, 1)];
   else
      C= (randn(L+1, 1) + j*randn(L+1, 1))/5;
      C= conv(C, flipud(C'.'));
      C= C./C(L+1);
      Cml(:, 1)= C;
   end;
   maxmax= max(max(abs(param_tdir(Cml, N))))
end;

figure(3);mesh(abs(tfarma_impr(Cml, 1, N, 1/2)))

while(maxmax > 1)
   Dml= win.*(randn(2*L+1, M+1) +j*randn(2*L+1, M+1))/10/M/(L+1)*maxmax;
   if(re_im=='r')
      Dml= param_hermite(Dml);
   end;
   if(mo_no=='m')
      Dml(:, 1)= zeros(2*L+1, 1);
   else
      Dml(:, 1)= param_hermite(Dml(:, 1));
   end;
   maxmax1= max(max(abs(param_tdir(Cml + Dml, N))));
   [maxmax maxmax1]
   if(maxmax1<maxmax)
      Cml= Cml + Dml;
      maxmax= maxmax1;
%      param_plot('exp', N, 0, 'r', 'g', Cml, 1, 1, 0, 2, 0, 0, 1, 0, 0, 0, 0, 0, 0)
      figure(3);mesh(fliplr(abs(tfarma_impr(Cml, 1, N, 1/2))))
   end;
   drawnow
end;

H= tfarma_impr(Cml, 1, N, 1/2);
nornor= norm(H(:, 1:N/2))
maxmax1= max(max(abs(param_tdir(Cml + Dml, N))));

while( (nornor > 1/4/sqrt(N)) || (maxmax > 1) )
   Dml= win.*(randn(2*L+1, M+1) +j*randn(2*L+1, M+1))/20/M/(L+1)*min(nornor, 1);
   if(re_im=='r')
      Dml= param_hermite(Dml);
   end;
   if(mo_no=='m')
      Dml(:, 1)= zeros(2*L+1, 1);
   else
      Dml(:, 1)= param_hermite(Dml(:, 1));
   end;
   H= tfarma_impr(Cml + Dml, 1, N, 1/2);
   nornor1= norm(H(:, 1:N/2));
   maxmax1= max(max(abs(param_tdir(Cml + Dml, N))));
   [nornor nornor1 1/4/sqrt(N)]
   if( (nornor1<nornor) && (maxmax1 < 1) )
      Cml= Cml + Dml;
      nornor= nornor1;
%      param_plot('exp', N, 0, 'r', 'g', Cml, 1, 1, 0, 2, 0, 0, 1, 0, 0, 0, 0, 0, 0)
      figure(3);mesh(fliplr(abs(tfarma_impr(Cml, 1, N, 1/2))))
   end;
   drawnow
end;
Cml
abs(Cml)
figure(3);mesh(fliplr(abs(tfarma_impr(Cml, 1, N, 1/2))))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;

for n= 8:8
N= 2^n;

for MAR= 4:4
   for LAR= 2:2
      tfpm_sysfile= sprintf('data/%04d/AR%d%d.mat', N, MAR, LAR)
      Amlr= tfpm_generate_stable(MAR, LAR, N, 'r', 'm');
      Amli= tfpm_generate_stable(MAR, LAR, N, 'i', 'm');
	    
      Anmr= param_expand(Amlr, N);
      Anmi= param_expand(Amli, N);

      save(tfpm_sysfile, ...
         'Amli' , 'Anmi' , ...
         'Amlr' , 'Anmr')
   end
end
      
for MMA= 5:5
   for LMA= 2:2
      tfpm_sysfile= sprintf('data/%04d/MA%d%d.mat', N, MMA, LMA)
      Bmlrm= tfpm_generate_stable(MMA, LMA, N, 'r', 'm');
      Bmlim= tfpm_generate_stable(MMA, LMA, N, 'i', 'm');
	    
      Bmlrn= tfpm_generate_stable(MMA, LMA, N, 'r', 'n');
      Bmlin= tfpm_generate_stable(MMA, LMA, N, 'i', 'n');

      Bnmrm= param_expand(Bmlrm, N);
      Bnmim= param_expand(Bmlim, N);

      Bnmrn= param_expand(Bmlrn, N);
      Bnmin= param_expand(Bmlin, N);

      save(tfpm_sysfile, ...
         'Bmlin', 'Bnmin', ...
         'Bmlrn', 'Bnmrn', ...
         'Bmlim', 'Bnmim', ...
         'Bmlrm', 'Bnmrm')
   end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%