gensurf_eq.m

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

% Generic integrated-error-fxn based surface plotting using the BERGulator.
% Valid only in the absence of noise.
 
% Copyright 1997-1998 Phil Schniter

function h_cont = gensurf_eq(C,F,diag_E,kappa,sigma2_n,delta_opt,org,ranges,alphabet,se_gamma,dse_gamma,ncma_alpha,alg,N_v)

 %global L0 L1 dg0 dg1 i k L a J g0_grid g1_grid grad_g0 grad_g1 m0 m1 f0_grid f1_grid

 % parse inputs
 alphabet = alphabet(:).';
 org = org(:).';
 if sigma2_n > 0,
   fprintf('\nWarning: calculation of cost contours ignores noise power;\n');
   fprintf('         comparison to Wiener receivers may be misleading.\n');
 end;

 % channel and equalizer lengths 
 [N_f,N_h] = size(F);           % BS/FS equalizer length

 % coordinate transformation 
 norm_f = zeros(1,N_f); 
 for i=1:N_f,
   norm_f(i) = norm(F(:,i));
 end;
 [dum,indx] = max(norm_f);
 g0 = F(:,indx)/norm(F(:,indx)); 
 g1 = [-g0(2);g0(1)];
 G = [g0,g1];
 range = abs([(sign(g0).*ranges)'*g0,(sign(g1).*ranges)'*g1]);

 % define a grid in transformed equalizer space 
 half_pts = 50;                  % usually use 40 
 g0_grid = linspace(-range(1),range(1),2*half_pts+1)+org(1);
 g1_grid = linspace(-range(2),range(2),2*half_pts+1)+org(2);
 L0 = length(g0_grid); L1 = length(g1_grid);
 dg0 = g0_grid(2) - g0_grid(1);
 dg1 = g1_grid(2) - g1_grid(1);

 % use finite-alphabet source... 
 M_s = length(alphabet);
 S = zeros(M_s^N_h,N_h);      	% every possible input vector
 for i=1:N_h,
   S(:,i) = alphabet(rem(floor([0:M_s^N_h-1]/M_s^(N_h-i)),M_s)+1).';
 end
 SC = S*C*G;			% every possible (transformed) regressor

 % modified regressors
 SCmod = zeros(M_s^N_h,N_f);
 switch alg,
   case {2,10},
     for i=1:M_s^N_h,
       SCmod(i,:) = SC(i,:)/(norm(SC(i,:))^2+ncma_alpha);  	% normalized 
     end
   case {4,5},
     SCmod = sign(real(SC))+j*sign(imag(SC));  			% signed 
 end;

 % compute gradient for one of various algorithms
 grad_g0 = zeros(L0,L1);
 grad_g1 = zeros(L0,L1);
 switch alg
 case {1,8},
   for i=1:L0,
     for k=1:L1,
       y_cur = SC*[g0_grid(i);g1_grid(k)];
       grad_cur = SC'*((abs(y_cur).^2 - kappa).*y_cur);		% CMA-22 
       grad_g0(i,k) = grad_cur(1);
       grad_g1(i,k) = grad_cur(2);
     end;
   end;
   title_str = 'CMA ';
 case 2,
   for i=1:L0,
     for k=1:L1,
       y_cur = SC*[g0_grid(i);g1_grid(k)];
       grad_cur = SCmod'*((abs(y_cur).^2- kappa).*y_cur);   	% NCMA-22 
       grad_g0(i,k) = grad_cur(1);
       grad_g1(i,k) = grad_cur(2);
     end;
   end;
   title_str = 'N-CMA ';
 case 3,
   for i=1:L0,
     for k=1:L1,
       y_cur = SC*[g0_grid(i);g1_grid(k)];
       grad_cur = SC'*sign((abs(y_cur).^2-se_gamma).*y_cur);	% SE-CMA 
       grad_g0(i,k) = grad_cur(1);
       grad_g1(i,k) = grad_cur(2);
     end;
   end;
   title_str = 'SE-CMA ';
 case 4,
   fprintf('\aWarning: Integrated SR-CMA surfaces usually inaccurate!\n\a');
   for i=1:L0,
     for k=1:L1,
       y_cur = SC*[g0_grid(i);g1_grid(k)];
       grad_cur = SCmod'*((abs(y_cur).^2 - kappa).*y_cur);	% SR-CMA
       grad_g0(i,k) = grad_cur(1);
       grad_g1(i,k) = grad_cur(2);
     end;
   end;
   title_str = 'SR-CMA ';
 case 5,
   for i=1:L0,
     for k=1:L1,
       y_cur = SC*[g0_grid(i);g1_grid(k)]; 			% SS-CMA 
       grad_cur = SCmod'*sign((abs(y_cur).^2-se_gamma).*y_cur);	
       grad_g0(i,k) = grad_cur(1);
       grad_g1(i,k) = grad_cur(2);
     end;
   end;
   title_str = 'SS-CMA ';
 case 6,
   for i=1:L0,
     for k=1:L1,
       y_cur = SC*[g0_grid(i);g1_grid(k)]; 			% DSE-CMA 
       grad_cur = SC'*min(max((abs(y_cur).^2-dse_gamma).*y_cur,-1),1);
       grad_g0(i,k) = grad_cur(1);
       grad_g1(i,k) = grad_cur(2);
     end;
   end;
   title_str = 'DSE-CMA ';
 case 7,
   fprintf('\aError: No surface available for WSGA.\n\a'); h_cont = pi;
   return;
 case 9,	% from baykal & tanrikulu
   for i=1:L0,
     for k=1:L1,
       y_cur = SC*[g0_grid(i);g1_grid(k)];
       grad_cur = SC'*((abs(y_cur) - gamGS).*y_cur);		% "GS" 
       grad_g0(i,k) = grad_cur(1);
       grad_g1(i,k) = grad_cur(2);
     end;
   end;
   title_str = 'GS ';
 case 10,	% from baykal & tanrikulu
   gamGS = sum(abs(alphabet).^3)/sum(abs(alphabet).^2);
   for i=1:L0,
     for k=1:L1,
       y_cur = SC*[g0_grid(i);g1_grid(k)];
       grad_cur = SCmod'*((abs(y_cur) - gamGS).*y_cur);		% SCS-1 
       grad_g0(i,k) = grad_cur(1);
       grad_g1(i,k) = grad_cur(2);
     end;
   end;
   title_str = 'SCS-1 ';
 end;% switch 
 grad_g0 = real(grad_g0);	% clean up numerical leftovers...
 grad_g1 = real(grad_g1);

 % integrate cost-gradient to obtain cost using coefficients "a"
 m0 = (L0+1)/2; 	% midpoint
 m1 = (L1+1)/2; 	% midpoint
 J = zeros(L0,L1);
 a = [0,55/24,-59/24,37/24,-9/24];	% third-order extrapolative
 a = [3/8,19/24,-5/24,1/24];		% third-order
 a = [1/2,1/2];				% first-order (trapezoid rule)
 L = length(a);
 for k=m1,			% integrate along "g0" direction
   for i=1:m0-(L-1),
     J(m0+i,k) = J(m0+i-1,k) + a*(grad_g0(m0+i+[-1:L-2],k))*dg0; 
     J(m0-i,k) = J(m0-i+1,k) - a*(grad_g0(m0-i-[-1:L-2],k))*dg0; 
   end
 end
 for i=(L-1):L0-(L-2),		% integrate along "g1" direction
   for k=1:m1-(L-1),
     J(i,m1+k) = J(i,m1+k-1) + (grad_g1(i,m1+k+[-1:L-2]))*a'*dg1; 
     J(i,m1-k) = J(i,m1-k+1) - (grad_g1(i,m1-k-[-1:L-2]))*a'*dg1; 
   end;
 end
 J = J-min(min(J));

 % derive 2D grid in f space
 f0_grid = zeros(L0,L1);
 for i=1:L1,
   f0_grid(:,i) = G(1,1)*g0_grid' + G(1,2)*g1_grid(i);  
 end;
 f1_grid = zeros(L0,L1);
 for i=1:L0,
   f1_grid(i,:) = G(2,1)*g0_grid(i) + G(2,2)*g1_grid;  
 end;

 % automatically select contour lines
 org_cost = J(m0,m1);
 wnr_cost = J( round(half_pts*(1+(F(1,3-delta_opt)-org(1))/ranges(1))),...
               round(half_pts*(1+(F(2,3-delta_opt)-org(2))/ranges(2))) );
 min_cost = 0;
 v = linspace( min_cost,max(0.9*org_cost,1.5*wnr_cost),N_v );
 %v = log10([[1:0.15:1.8],[2:0.4:5]]);
 %v = [linspace(0.02,0.25*J(m0,m1),6),...
 %     logspace(log10(0.3*J(m0,m1)),log10(0.9*J(m0,m1)),8)];
 %v = linspace(0,3.6*J(m0,m1),N_v);
  v = linspace(0,0.9*J(m0,m1),N_v);

 % plot cost contours in equalizer space
 [dum, h_cont] = contour(f0_grid,f1_grid,J,v,'c'); h_cont = h_cont.'; 
 axis([-ranges(1),ranges(1),-ranges(2),ranges(2)]);
 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(['Integrated ',title_str,'Cost Surface']);

 % plot MSE ellipse axes
 [V,Lam] = eig(C'*C+sigma2_n*eye(N_f));
 vec0 = V(:,1)*sqrt(1./max(Lam(1,1),eps))/6;    % protection for singular chans
 vec1 = V(:,2)*sqrt(1./max(Lam(2,2),eps))/6; 
 axe_org = [-(3/4)*ranges(1),(3/4)*ranges(2)] + org;
 hold on;
 h_cont = [h_cont,...
   plot(axe_org(1)+[-vec0(1),vec0(1)],axe_org(2)+[-vec0(2),vec0(2)]),...
   plot(axe_org(1)+[-vec1(1),vec1(1)],axe_org(2)+[-vec1(2),vec1(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;
 zoom on;

return

 grid on;
 figure(2); clf;
 %quiver(g0_grid,-g1_grid,-grad_g1,grad_g0,10) 
 contour(g0_grid,g1_grid,((grad_g0.^2+grad_g1.^2).^0.5)',linspace(0,4.0,20) ) 
 xlabel('g0'); ylabel('g1');
 zoom on
 axis('equal')
 grid on