tf_multiwin.m

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

function [Psi, Mask, v2]= tf_multiwin(N, M, L, K, cut, proto)
% function [Psi, Mask, v2]= tf_multiwin(N, M, L[, K, cut, proto])
%   This file is part of the TFPM toolbox v1.0 (c)
%   [e-mail] and underlies the GPL.
% 
% Generates a positive semidefinite (nonseparable) TF window starting
% from a prototype:
%   0: Hanning/Hanning
%   1: Gaussian
%   2: Choi-Williams?
% The variable cut means:
%   0: nothing happens
%   1: the prototype is multiplied by a rect(M, L)
%   2: "   "         "  "               rolloff(M, L).
%
% The largest positive K eigenvalues are used to form the
% posdef. function, alpha= 1/2! If K==0 the K= N/M/L. Kmax=
% sqrt(N), Mmax, Lmax= N/2-1!
%

% Stefan: Psi= tf_multiwin(N, M, L, 0, 2, 1)!!!!!!!!!!!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    = 128;
M    =   3;
L    =   3;
K    =   0;
cut  =   0;
proto=   0;


Psi= tf_multiwin(N, M, L);
tf_show(Psi)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha= 1/2;

if( ~M && ~L)
   Psi= ones(N);
   Mask= Psi;
   v2= [];
   return;
end;

if(M>N/2-1)
   M= N/2-1;
end;

if(L>N/2-1)
   L= N/2-1;
end;

if(nargin<=5)
   proto= 1;
end;

if(nargin<=4)
   cut= 2;
end;

if(nargin<=3)
   K= 0;
end;

if(~K)
   K= max(1, floor(2*N/(M+1)/(L+1)));
end;
if(K>sqrt(N))
   K= floor(sqrt(N));
end;
%if(K<3)
%   K= 3;
%end;

Mask= zeros(N);
switch(proto)
 case 0% Hanning/Hanning
  L1= ceil(0.8*N/L);
  L1= ceil(L1/3*2);M= ceil(1.5*M);
  winm= [zeros(1, N/2-M) hanning(2*M+1).' zeros(1, N/2-M-1)];
  winl= [zeros(N/2-L1, 1); hanning(2*L1+1); zeros(N/2-L1-1, 1)];
  Winl= fft([winl(N/2+1:N); winl(1:N/2)]);
  Winl= real([Winl(N/2+1:N); Winl(1:N/2)]);
  Mask= Winl*winm;
 case 1% Gaussian
  if(M)
     Maskm= exp(-(-N/2:N/2-1).^2/M^2);
  else
     Maskm= ones(1, N);
  end;
  if(L)
     Maskl= exp(-(-N/2:N/2-1).^2/L^2).';
  else
     Maskl= ones(N, 1);
  end;
 Mask= Maskl*Maskm;
 case 2% Choi-Williams?NOT WORKING NOW!!!
  if(~M)
     M= 1;
  end;
  if(~L)
     L= 1;
  end;
  for m= -N/2:N/2-1
     for l= -N/2:N/2-1
        Mask(l+1+N/2, m+1+N/2)=exp(-m^2*l^2/M^2/L^2);
     end
  end
end;

switch(cut)
 case 0
  Mask1= ones(N);
 case 1
  Mask1= zeros(N);
  Mask1(N/2+1-L:N/2+1+L, N/2+1-M:N/2+1+M)= ones(2*L+1, 2*M+1);
 case 2
  Mask1= zeros(N);
  Mask1(N/2+1-L:N/2+1+L, N/2+1-M:N/2+1+M)= ones(2*L+1, 2*M+1);
  Mask1(N/2+1-L:N/2+1+L, N/2+1-M)= ones(2*L+1, 1)/2;
  Mask1(N/2+1-L:N/2+1+L, N/2+1+M)= ones(2*L+1, 1)/2;
  Mask1(N/2+1-L, N/2+1-M:N/2+1+M)= ones(1, 2*M+1)/2;
  Mask1(N/2+1+L, N/2+1-M:N/2+1+M)= ones(1, 2*M+1)/2;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(94);clf;cesh(Mask)
figure(95);clf;imagesc(abs(Mask))
figure(96);clf;cesh(Mask1)
figure(97);clf;imagesc(Mask1)
figure(98);clf;cesh(Mask.*Mask1)
figure(99);clf;imagesc(Mask.*Mask1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Mask= Mask.*Mask1;
Mask= (Mask+sym_ambi(Mask))/2;%:-)

psi= ml_to_nm(Mask);
psi= lag_to_ker(psi, 1/2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(97);clf;tf_show(psi);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[U, V]= eig(psi);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[max(max(abs(imag(U)))) max(max(abs(imag(V))))]
figure(98);clf;tf_show(U)
figure(99);clf;clot(diag(V))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

V2 = real(diag(pospart(V)));
v2= flipud(sort(V2));
v2= v2(1:K);
V2= V2.*(V2>=min(v2));
psi2 = U*diag(V2)*U';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(97);clf;tf_show(psi2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Psinm= ker_to_lag(psi2, -1, 1/2);
Psiml= nm_to_ml(Psinm);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(97);clf;tf_show(Psinm);
figure(98);clf;tf_show(Psiml);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Psi= Psiml/Psiml(N/2+1, N/2+1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Psi(N/2+1, N/2+1)
figure(97);clf;cesh(Psi)
figure(98);clf;cesh(nm_to_nk(ml_to_nm(Psi)))
figure(99);clf;cesh(Psi-sym_ambi(Psi))
norm(Psi-sym_ambi(Psi))


N=128;
M=  4;
L=  4;
cut= 0;
[Psi, Mask, v2]= tf_multiwin(N, M, L, 0, cut, 2);

v2
figure(1);clf;
imagesc(abs(Psi))
set(gca, 'XTick', [1 N/2+1-M N/2+1 N/2+1+M N])
set(gca, 'XTicklabel', [-N/2 -M 0 M N/2-1])
set(gca, 'YTick', [1 N/2+1-L N/2+1 N/2+1+L N])
set(gca, 'YTicklabel', [-N/2  -L 0 L N/2-1])
colormap(flipud(gray.^2))
%print('-depsc2', '-tiff', '-r300', sprintf('/users/mjachan/tex/prj/TF/figures/kernel%d%d%d1', M, L, cut))

figure(2);clf;
subplot(2, 1, 1);
plot(abs(Psi.'), 'Color', 'k')
axis([1 N 0 1.1])
set(gca, 'XTick', [1 N/2+1-M N/2+1 N/2+1+M N])
set(gca, 'XTicklabel', [-N/2 -M 0 M N/2-1])

subplot(2, 1, 2);
plot(abs(Psi), 'Color', 'k')
axis([1 N 0 1.1])
set(gca, 'XTick', [1 N/2+1-L N/2+1 N/2+1+L N])
set(gca, 'XTicklabel', [-N/2  -L 0 L N/2-1])
%print('-depsc2', '-tiff', '-r300', sprintf('/users/mjachan/tex/prj/TF/figures/kernel%d%d%d2', M, L, cut))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N= 128
figure(94);clf;imagesc(abs(tf_multiwin(N,  0,  0, 10, 2, 1)))
figure(95);clf;imagesc(abs(tf_multiwin(N,  0,  1, 10, 2, 1)))
figure(96);clf;imagesc(abs(tf_multiwin(N,  0,  2, 10, 2, 1)))
figure(97);clf;imagesc(abs(tf_multiwin(N,  1,  0, 10, 2, 1)))
figure(98);clf;imagesc(abs(tf_multiwin(N,  1,  1, 10, 2, 1)))
figure(99);clf;imagesc(abs(tf_multiwin(N,  1,  2,  6, 2, 1)))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N= 128
figure(94);clf;cesh(tf_multiwin(N,  0,  0, 10, 2, 1))
figure(95);clf;cesh(tf_multiwin(N,  0,  1,  9, 2, 1))
figure(96);clf;cesh(tf_multiwin(N,  0,  2,  8, 2, 1))
figure(97);clf;cesh(tf_multiwin(N,  1,  0,  7, 2, 1))
figure(98);clf;cesh(tf_multiwin(N,  1,  1,  6, 2, 1))
figure(99);clf;cesh(tf_multiwin(N,  1,  2,  5, 2, 1))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N= 128
for M= 0:10
   for L= 0:10
      [M, L]
      [Psi, Mask, v2]= tf_multiwin(N, M, L, 0, 2, 1);
      figure(94);clf;cesh(Psi)
      figure(95);clf;imagesc(abs(Psi))
      figure(96);clf;plot(v2);
      pause
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
MAR  =   3;
LAR  =   3;
MMA  =   0;
LMA  =   0;
N    =  64;
re_im= 'r';
mo_no= 'm';
tfpm_file_gen;
%-------------
figure(94);clf;mesh(abs(tfarma_ambi(Aml, Bml, N, 1/2)));
figure(95);clf;mesh(tfarma_wvsp(Aml, Bml, N, 1/2));
Psi= tf_multiwin(N, 2*MAR, 2*LAR, 0, 2, 0);
figure(96);clf;mesh(abs(Psi));

y= tfarma_gen(randn(N, 1), Aml, Bml, 1/2);
Qalpha= corr_est(y, y, -1, 1/2);
Aalpha= fft(Qalpha);
Aalpha= [Aalpha(N/2+1:N, :); Aalpha(1:N/2, :)];
%figure(97);clf;mesh(abs(Aalpha))
Awlpha= Aalpha.*conj(Psi);
figure(98);clf;mesh(abs(Awlpha))
%Palpha= nm_to_nk(ml_to_nm(Aalpha));
%figure(98);clf;cesh(Palpha)
Pwlpha= real(nm_to_nk(ml_to_nm(Awlpha)));
figure(99);clf;mesh(Pwlpha)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
N= 64;
M= 8;
L= 8;

Psi0= tf_multiwin(N, M, L, 0, 0, 0);
Psi1= tf_multiwin(N, M, L, 0, 0, 1);
figure(1);tf_show(Psi0)
figure(2);tf_show(Psi1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;% TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%