% Zero forcing equalizer with additive noise in the received signals.
% 
% Rich Kozick, ELEC 470, Spring 1998

T = 1;      % Bit period
tau = 3;    % Time constant of channel
SNR = 100;   % Ratio of signal power to noise power (NOT in dB)

dt = 0.01;  % Sampling time in simulation
N = 100;    % Number of bits to generate

clear t1 t2 c x y

% Create output pulse: rectangular pulse convolved with first-order
% low-pass filter impulse response.

t1 = (dt:dt:T)';
c(1:100,1) = 1 - exp(-t1/tau);
t2 = (T+dt:dt:T+5*tau)';
c(101:100+length(t2),1) = c(100)*exp(-(t2-T)/tau);

figure(1)
plot([t1; t2], c)
xlabel('Time (sec)')
ylabel('c(t)')
title('Smeared pulse c(t)')

% Generate bit stream

b = rand(N,1);
z0 = find(b < 0.5);
z1 = find(b >= 0.5);
b(z0) = -1*ones(size(z0));
b(z1) = +1*ones(size(z1));

% Create received signal with ISI

nT = T/dt;
nc = length(c);
nx = N*nT;
x = zeros(nx, 1);
for n=1:N
  i1 = (n-1)*nT;
  y = [zeros(i1,1); b(n)*c; zeros(N*nT-i1-nc,1)];
  x = x + y(1:nx);
end

% Add noise of the specified level

sp = sum(x.*x)/length(x);  % Signal power
np = sp/SNR;               % Noise power
noise = sqrt(np) * randn(length(x),1);
x = x + noise;

% Plot eye diagram

figure(2)
t3 = dt:dt:2;
plot(t3, x(1:200));
hold on
for n=3:2:N
  plot(t3, x((n-1)*nT+1:(n+1)*nT));
end
hold off
xlabel('Time (sec)')
title('EYE DIAGRAM WITH NO EQUALIZATION')

% Compute number of bit errors

xT = x(nT:nT:nx);
dz0 = find(xT < 0);
dz1 = find(xT >= 0);

db = b;
db(dz0) = -1*ones(size(dz0));
db(dz1) = +1*ones(size(dz1));

err = find(db ~= b);

fprintf('No equalizer: %d bits out of %d in error\n', length(err), N);

% Define length of zero-forcing equalizer

Ne = 5*tau/T;

% Get samples of c(t) to solve for ZF equalizer weights

cT = c(nT:nT:nc);
csamp = [zeros(2*Ne,1); cT; zeros(2*Ne+1-length(cT),1)];

% Construct the matrix on the left side of ZF equalizer equation

C = zeros(2*Ne+1,2*Ne+1);
for ne = 1:2*Ne+1
  C(ne,:) = csamp(2*Ne+ne:-1:ne)';
end

% Right side of ZF equalizer weight equation

r = [zeros(Ne,1); 1; zeros(Ne,1)];

% Solve for ZF equalizer weights

w = C \ r;

% Process the received signal with the ZF equalizer

nw = length(w);
z99 = [1; zeros(nT-1,1)];
hzf = kron(w, z99);        % Impulse response of equalizer
yall = conv(x, hzf);       % Do the equalization filtering
y = yall((Ne*nT+1):(length(yall)-(Ne+1)*nT)+1);


% Eye diagram of equalized signal

figure(3)
plot(t3, y(1:200));
hold on
for n=3:2:N
  plot(t3, y((n-1)*nT+1:(n+1)*nT));
end
hold off
xlabel('Time (sec)')
title('EYE DIAGRAM AFTER ZF EQUALIZER')

% Compute number of bit errors

yT = y(nT:nT:nx);
dz0 = find(yT < 0);
dz1 = find(yT >= 0);

db = b;
db(dz0) = -1*ones(size(dz0));
db(dz1) = +1*ones(size(dz1));

err = find(db ~= b);

fprintf('ZF equalizer: %d bits out of %d in error\n', length(err), N);