cmasurf_eq.m

有关信道估计和信道均衡的仿真程序

% This plots the CMA surface for the BERGulator and returns a handle to it
 
% Copyright 1997-1998 Phil Schniter

function h_cont = cmasurf_eq(C,F,diag_E,kappa,sigma2_n,delta_opt,org,ranges,real_s,real_n,N_v)
 
 % channel and equalizer lengths 
 [N_f,N_h] = size(F);		% BS/FS equalizer length

 % define a grid in f space 
 num_pts = 80;			% usually use 80 (or 150)
 org = org(:).';
 f1_grid = linspace(-ranges(1),ranges(1),num_pts)+org(1);
 f2_grid = linspace(-ranges(2),ranges(2),num_pts)+org(2);
 L1 = length(f1_grid); L2 = length(f2_grid);

 % calculate CMA cost values on f grid  (sigma2_s=1)
 CMA = zeros(L1,L2);
 if real_s&real_n		% PAM source with real noise
   for i=1:L1,
     for k=1:L2,
       f = [f1_grid(i);f2_grid(k)];
       h = C*f;
       hTh = h'*h; 		fTf = f'*f;
       CMA(i,k) = kappa^2 - 2*kappa*(hTh + sigma2_n*fTf)...
	 + 3*hTh^2 + (kappa-3)*(h.^2)'*(h.^2)...
	 + 6*sigma2_n*hTh*fTf + 3*sigma2_n^2*fTf^2;
     end
   end
 elseif real_s,			% PAM source with rot.inv. noise
   for i=1:L1,
     for k=1:L2,
       f = [f1_grid(i);f2_grid(k)];
       h = C*f;
       hTh = h'*h; 		fTf = f'*f;
       CMA(i,k) = (kappa-3)*sum(abs(h).^4) + 2*hTh^2 + abs(h.'*h)^2 ...
	+ 2*sigma2_n^2*fTf^2 + 4*sigma2_n*hTh*fTf...
	- 2*kappa*(hTh + sigma2_n*fTf) + kappa^2;
     end
   end
 else				% QAM source with rot.inv. noise
   for i=1:L1,
     for k=1:L2,
       f = [f1_grid(i);f2_grid(k)];
       h = C*f;
       hTh = h'*h; 		fTf = f'*f;
       CMA(i,k) = (kappa-2)*sum(abs(h).^4) + 2*hTh^2 ...
	+ 2*sigma2_n^2*fTf^2 + 4*sigma2_n*hTh*fTf...
	- 2*kappa*(hTh + sigma2_n*fTf) + kappa^2;
     end
   end
 end

 % automatically select contour lines
 org_cost = CMA(num_pts/2,num_pts/2);
 [dum,delta_wrst] = max(diag_E);
 wnr_cost = CMA( round(num_pts/2*(1+(F(1,delta_wrst)-org(1))/ranges(1))),...
		    round(num_pts/2*(1+(F(2,delta_wrst)-org(2))/ranges(2))) );
 min_cost = min(min(CMA));
 v = linspace( min_cost,max(0.9*org_cost,1.5*wnr_cost),N_v );
%v = min(min(CMA))+log10([[1:0.15:1.8],[2:0.4:5]]); 
%v = min(min(CMA))+log10([[1:0.005:1.02],[1.03:0.15:1.8],[2:0.4:10]]);
%v = [linspace(0.02*min(min(CMA)),0.25*CMA(num_pts/2,num_pts/2),6),...
%     logspace(log10(0.3*CMA(num_pts/2,num_pts/2)),...
%     log10(0.7*CMA(num_pts/2,num_pts/2)),6)];

 % plot CMA contours in f space
 [dum,h_cont] = contour(f1_grid,f2_grid,CMA',v,'c'); h_cont = h_cont.';
%[dum,h_cont] = contour(f1_grid,f2_grid,CMA',v); h_cont = h_cont.';
 grid off;
 hold on;
 for delta = 1:N_h, 	% plot optimal equalizers for delta delay
   h_cont = [h_cont, plot(F(1,delta),F(2,delta),'x'),... 
   	     plot(-F(1,delta),-F(2,delta),'x')]; 
 end;
 for delta = 1:N_h, 	% plot optimal equalizers for delta delay
   if abs(diag_E(delta) - diag_E(delta_opt)) < 100*eps,
     h_cont = [h_cont, plot(F(1,delta),F(2,delta),'*'),... 
               plot(-F(1,delta),-F(2,delta),'*')]; 
   end;
 end;
 hold off;
 xlabel('f_0'); ylabel('f_1'); 
 title('CMA Error Surface');

 % plot MSE ellipse axes
 axe_org = [-(3/4)*ranges(1),(3/4)*ranges(2)] + org;
 [V,Lam] = eig(C'*C+sigma2_n*eye(N_f));
 vec1 = V(:,1)*sqrt(1./max(Lam(1,1),eps))/6; 	% protection for singular chans
 vec2 = V(:,2)*sqrt(1./max(Lam(2,2),eps))/6; 
 hold on;
 h_cont = [h_cont,...
   plot(axe_org(1)+[-vec1(1),vec1(1)],axe_org(2)+[-vec1(2),vec1(2)]),...
   plot(axe_org(1)+[-vec2(1),vec2(1)],axe_org(2)+[-vec2(2),vec2(2)]),...
   text((1.21)*(axe_org(1)-org(1))+org(1),(1.20)*(axe_org(2)-org(2))+org(2),...
	'MSE ellipse axes')];
 hold off;