channel_estold.m

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
alpha= 1/2;
load updown_h_rusk.mat;
N= 2048;
M0= 28;
figure(1);mesh(abs(h_I2))

MCH= 5;
%I= 1;
Mmax= 1;
Lmax= 1;

hnm= h_I2(113:end, M0:M0+MCH);
hnm= hnm/max(max(abs(hnm)));

figure(2);mesh(abs(hnm))
AMBI= {};
for tau= 0:MCH
   for taup= 0:MCH
%      if(abs(tau-taup)<=I)
	 [tau taup]
         Rttp= corr_est(hnm(:, tau+1), hnm(:, taup+1), Mmax, alpha);
         Attp= fft(Rttp);
         Attp= [Attp(N/2+1:N, :); Attp(1:N/2, :)];
         mesh(abs(Attp(N/2+1-3*Lmax:N/2+1+3*Lmax, :)))
         drawnow
         AMBI{tau+1, taup+1}= Attp(N/2+1-3*Lmax:N/2+1+3*Lmax, :);
%      end;
   end;
end;
filename= sprintf('channelAMBI%d%d', MCH, MCH)
save(filename, 'AMBI', 'M0', 'hnm', 'N');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm
alpha= 1/2;
load updown_h_rusk.mat;

MCH= 5;
I= 5;
filename= sprintf('channelAMBI%d%d', MCH, I)
load(filename);
[L, M]= size(AMBI{1, 1});
M= (M-1)/2;
L= (L-1)/6;

A= zeros(M*(2*L+1), (MCH+1)*M*(2*L+1));
for tau= 1:MCH
   for l= -L:L
      for lp= -L:L
	 for m= 1:M
	    for mp= 1:M
%	       [tau m mp l lp]
	       [(m-1)*(2*L+1)+L+1+l tau*(M*(2*L+1))+(mp-1)*(2*L+1)+L+1+lp]
	       A((m-1)*(2*L+1)+L+1+l, tau*(M*(2*L+1))+(mp-1)*(2*L+1)+L+1+lp)=...
		   AMBI{tau+1, tau}(3*L+1+l-lp, M+1+m-mp);
	    end;
	 end;
      end;
   end;
end;

B= zeros(M*(2*L+1), (MCH+1)*M*(2*L+1));
for tau= 0:MCH
   for l= -L:L
      for lp= -L:L
	 for m= 1:M
	    for mp= 1:M
%	       [tau m mp l lp]
	       [(m-1)*(2*L+1)+L+1+l tau*(M*(2*L+1))+(mp-1)*(2*L+1)+L+1+lp]
	       B((m-1)*(2*L+1)+L+1+l, tau*(M*(2*L+1))+(mp-1)*(2*L+1)+L+1+lp)=...
		   AMBI{tau+1, tau+1}(3*L+1+l-lp, M+1+m-mp);
	    end;
	 end;
      end;
   end;
end;

C= zeros(M*(2*L+1), (MCH+1)*M*(2*L+1));
for tau= 0:MCH-1
   for l= -L:L
      for lp= -L:L
	 for m= 1:M
	    for mp= 1:M
%	       [tau m mp l lp]
	       [(m-1)*(2*L+1)+L+1+l tau*(M*(2*L+1))+(mp-1)*(2*L+1)+L+1+lp]
	       C((m-1)*(2*L+1)+L+1+l, tau*(M*(2*L+1))+(mp-1)*(2*L+1)+L+1+lp)=...
		   AMBI{tau+1, tau+2}(3*L+1+l-lp, M+1+m-mp);
	    end;
	 end;
      end;
   end;
end;

a= zeros((MCH+1)*M*(2*L+1), (MCH+1));
for tau= 0:MCH
   for taup= 0:MCH
      if(abs(tau-taup)<=I)
         for m= 1:M
            for l= -L:L
               a(tau*(M*(2*L+1))+(m-1)*(2*L+1)+L+1+l, taup+1)= AMBI{tau+1, taup+1}(3*L+1+l, M+1+m);
	    end;
	 end;
      end;
   end;
end;

Ainv= invert_blocktri(A, B, C);

theta= -Ainv*a;

CML= {};
for tau= 0:MCH
   for taup= 0:MCH
      [tau, taup]
      Aml= [[zeros(L, 1); 1; zeros(L, 1)] ...
	    param_destack_ml(theta(tau*M*(2*L+1)+1:(tau+1)*M*(2*L+1), taup+1), M-1, L)];
      Aml= param_stabilize(Aml, N, .98, 1);
      CML{tau+1, taup+1}= {Aml 1};
   end;
end;

filename= sprintf('channelPARA%d%d', MCH, I)
save(filename, 'CML', 'M0', 'hnm', 'N');