% costsurf.m script  -- cost surface generation via Eq (7)
%
% This script displays the cost contours for SE-CMA with 2-tap 
% fractionally spaced equalizers given the following input parameters:
%
%     c       a vector of channel coefficients
%     SNR     the SNR of the system
%     M       the alphabet size
%     src_var the source variance
%
%
% The following parameters affect the display and accuracy of the graph:
%
%     f0             a vector of the values of f0 over which the program should iterate
%     f1             a vector of the values of f1 over which the program should iterate
%                      (Note: f0 and f1 must have the same length)
%     contour_lines  the number of contour lines to display
%
%
% The following flags affect WHAT is displayed on the graph:
% 
%     USE_NOISELESS_GAMMA  determines which way to calculate gamma, i.e. via the 
%                                Schniter/Brown paper (1) or via the Ghosh paper (0)
%     SHOW_GRANULARITY     flag to display yellow "granularity" dots that show granularity
%                                of possible grpah locations for SE-CMA and CMA 2-2 minima
%     SHOW_CMA22_MIN       flag to display CMA 2-2 minima (green diamonds)
%     SHOW_SECMA_MIN       flag to display SE-CMA minima  (blue squares)
%     SHOW_MMSE_MIN        flag to display MMSE minima    (red stars, with optimal minima circled)
%
% All of the input parameters are defined in the beginning of the script file.

% A. Klein  2/28/98

% -------------------------------- user defined variables and flags
c=[0.2 0.5 1 -0.1]';
%c=[0.1 0.3 1 -0.1 0.5 0.2]';
SNR=15;
M=8;
src_var=1;

num_pts=30;
contour_lines=30;
f0=linspace(-4,4,num_pts);
f1=linspace(-2,2,num_pts);

USE_NOISELESS_GAMMA=1;
SHOW_GRANULARITY=1;
SHOW_CMA22_MIN=1;
SHOW_SECMA_MIN=1;
SHOW_MMSE_MIN=1;


% -------------------------------- generate source, noise, channel matrix
if USE_NOISELESS_GAMMA
  beta=sqrt(SEgamma(M,src_var));
else
  beta=sqrt(SEgamma_noisy(M,src_var,SNR));
end;
gamma=(9*M^2-21)/(5*M^2-5);  %closed form expression for kurtosis of PAM sequence -- for CMA 2-2
w_var=src_var*10^(-SNR/10);
Ns=length(c)/2;
s=allpam(M,Ns,src_var);
C=convmtx(c,2);
C=C(2:2:length(C),:);
m=s'*C;

% -------------------------------- create constants
sig=sqrt(src_var*10^(-SNR/10));
Sq2=sqrt(2);   % constant
SqPi=sqrt(pi); % constant

% -------------------------------- iterate through range of f and create cost surface
% -------------------------------- matrices for CMA 2-2 and CMA 1-1
SE=zeros(num_pts);
CMA=zeros(num_pts);
for i0=1:num_pts
   for i1=1:num_pts
      f=[f0(i0); f1(i1)];
      x=m*f;
      h=conv(c,f);
      h=h(2:2:length(h));
      sig=sqrt(w_var)*norm(f);
      SE(i1,i0)=sum(-beta+(beta+x).*erf((beta+x)/sig/Sq2)+(beta-x).*erf((beta-x)/sig/Sq2)...
                   -x.*erf(x/sig/Sq2)+Sq2/SqPi*sig*(exp(-(beta+x).^2/2/sig^2)...
                   +exp(-(beta-x).^2/2/sig^2)-exp(-x.^2/2/sig^2)))*M^-Ns;
      CMA(i1,i0)=src_var^2*(gamma-3)*sum(h.^4)+3*src_var^2*norm(h)^4 ...
                   +3*w_var^2*norm(f)^4+6*src_var*w_var*norm(f)^2*norm(h)^2 ...
                   -2*src_var*gamma*(src_var*norm(h)^2+w_var*norm(f)^2)+src_var^2*gamma^2;
   end
end

% --------------------------------- plot resulting contour
lambda=10^(-SNR/10);
temp=inv(C'*C+lambda*eye(2))*C';

f_mmse=zeros(2,Ns); mse=zeros(1,Ns); h_mtx=eye(Ns);
for k=1:Ns
   f_mmse(:,k)=temp*h_mtx(:,k);
   mse(k)=src_var*norm(h_mtx(k)-C*f_mmse(:,k))^2+(src_var*10^(-SNR/10))*norm(f_mmse(:,k))^2;
end

contour(f0,f1,SE,linspace(min(min(SE)),0.9*SE(num_pts/2,num_pts/2),contour_lines))
title('SE-CMA cost surface'); xlabel('f_0'); ylabel('f_1');
hold on

if SHOW_MMSE_MIN
  for k=1:Ns
    plot(f_mmse(1,k),f_mmse(2,k),'rp');
    plot(-f_mmse(1,k),-f_mmse(2,k),'rp');
    if mse(k)==min(mse)
        plot(f_mmse(1,k),f_mmse(2,k),'ro');
        plot(-f_mmse(1,k),-f_mmse(2,k),'ro');
    end
  end
end

if SHOW_SECMA_MIN
  minSE=surfmin(SE);
  for k=1:size(minSE,1)
    plot(f0(minSE(k,2)),f1(minSE(k,1)),'bs')
  end
end

if SHOW_CMA22_MIN
  minCMA=surfmin(CMA);
  for k=1:size(minCMA,1)
    plot(f0(minCMA(k,2)),f1(minCMA(k,1)),'gd')
  end
end

if SHOW_GRANULARITY
  for i=1:length(f1)
    plot(f0,f1(i)*ones(length(f0),1),'y.')
  end
end

hold off