% simulate_cma.m
% This file generates received T/-spaced, noisy signals from
% empirically-measured T/2-spaced AppSigTec channels
% and simulates a blind, adaptive equlaizer updated with CMA.
% The source sequence is QAM.
% All user selectable variables are in the Input Parameters section,
% except for the channel which is loaded in the Load and Normalize
% Channel section.
% tje 6/12/96

clear
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% INPUT PARAMETERS %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Lg=16; % equalizer length
n=5000; % number of source symbols generated
snrdb=35; % SNR (in dB) at Eq. input
SNR=10.^(snrdb/10);
d=4; % d^2-level QAM

channel_selector=16; % 0==> Use full, T/2 channel, no frequency decimation
% xx==> Use frequency-decimated, length xx channel
% xx should be a power of 2

mu=1e-3; % stepsize
c=zeros(Lg,1); % Equalizer INITIALIZATION
c(floor(Lg/2))=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% LOAD and NORMALIZE CHANNEL %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load chan2
if channel_selector==0
chan=C;
clear C rt
else
freq_full=fft(C,256);
freq_dec=freq_full(1:256/channel_selector:256);
chan=ifft(freq_dec);
clear C rt freq_full freq_dec
end

Lc=length(chan);rand('seed',sum(100*clock))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% QAM Source Sequence %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rand('seed',sum(100*clock))
alphabet=-(d-1):2:(d-1);
alphabet=(alphabet/sqrt(mean(alphabet.^2)));

x1=d*rand(n/2,1);
x1=floor(x1)+1;
x2=d*rand(n/2,1);
x2=floor(x2)+1;
s=zeros(n,1);
s(1:2:n)=alphabet(x1)+j*alphabet(x2);
clear x1 x2
s=sqrt(0.5)*s; % T-spaced sequence is unit variance
% T/2-spaced sequence has variance 0.5

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Received Sequence %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rcvd=filter(chan,1,s); % received sequence
clear s
noise=randn(length(rcvd),2)*([1 j].'); % Complex Noise
noise=noise*sqrt(1/SNR)*(1/sqrt(2));
rcvd=rcvd+noise;
clear noise

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Dispersion Constant: QAM Source %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alphabet=-(d-1):2:(d-1); % values for PAM letters
beta=alphabet*j;
constellation=100*ones(d^2,1);
ccccc=0;
for pp=1:d
for qq=1:d
ccccc=ccccc+1;
constellation(ccccc)=alphabet(pp)+beta(qq);
end
end
constellation=constellation/sqrt(mean(abs(constellation).^2));
rd=mean(abs(constellation).^4)/mean(abs(constellation).^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% INITIALIZE STORAGE VECTORS %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% assume 500 data points saved across eq. iteration
isie=100*ones(500,1);
isio=100*ones(500,1);
nge=100*ones(500,1);
ngo=100*ones(500,1);
residuale=100*ones(500,1);
residualo=100*ones(500,1);
zz=n/1000; % keep 500 data points across iteration run
ckeep=c*ones(1,500);
ykeep=zeros(500,1);
iter=zeros(500,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% CMA BLIND ADAPTIVE EQUALIZER %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

count=0;
for i=Lg+1:2:n-Lg+1 %update on baud
kkk=(i-1)/2;
regg=rcvd(i:-1:i-Lg+1); %regressor of inputs to equalizer
y=regg.'*c; % equalizer output vector
cmaerror=y*(rd-(abs(y))^2);
c=c+mu*conj(regg)*cmaerror;

% save zzth data points
if rem(kkk,zz)==0
count=count+1;
ykeep(count)=y; % equalizer output vector
ckeep(:,count+1)=c;
iter(count)=i;

end % save data loop
end % received data loop

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% EQUALIZER PERFORMANCE %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% (stupid way to do this, but is out of the iteration loop)
for lll=1:500
jjj=conv(chan,ckeep(:,lll)); % Combined T/2-spaced Chan-Eq. Response
jjje=jjj(1:2:length(jjj)); % Downsampled (even) Chan-Eq response
jjjo=jjj(2:2:length(jjj)); % Downsampled (odd) Chan-Eq response
isie(lll)=isi(jjje);
isio(lll)=isi(jjjo);
nge(lll)=sum(abs(ckeep(:,lll).^2));
ngo(lll)=sum(abs(ckeep(:,lll).^2));
residuale(lll)=residual(jjje);
residualo(lll)=residual(jjjo);
end

mse_e=residuale + nge*(1/SNR);
mse_o=residualo + ngo*(1/SNR);

semilogy(iter(1:490),mse_e(1:490)) % Neglect last samples/edge effects

%%%%%%%%%%%%%%%%%%%%%%%
%%% SAVE DATA %%%
%%%%%%%%%%%%%%%%%%%%%%%
clear r
complexity=flops;
%save cma_output_data