% This plots the CMA surface for the BERGulator and returns a handle to it % Copyright 1997-1998 Phil Schniterfunction 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;