% Generic integrated-error-fxn based surface plotting using the BERGulator.% Valid only in the absence of noise. % Copyright 1997-1998 Phil Schniterfunction 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