%--------------------------- Begin Lab 5 --------------------
% Colin Joye
% ECE 3770-01
% Lab 5, Decision Directed LMS
% due 4/26/01R
clear all;

% 1) Setup.
n  = 1:2000; % Sample vector [samples]
Ts = 10^-3;  % Sampling period [s]
f1 = 1;      % Message frequency [Hz] 
in = 2*(wgn(1,length(n),0)>0)-1; % in-phase
quad = 2*(wgn(1,length(n),0)>0)-1; % quadrature
x = in+j*quad; % clean QPSK message transmitted.
snr = 10;  % Signal to Noise Ratio

% 2) Plotting the message free from noise or channel.
figure(1)
subplot(211);
scatter(real(x), imag(x),'*'); grid;
axis([-2 2 -2 2]);
ylabel('Quadrature');
s=['QPSK Constellation, No Noise'];
title(s);

% 3) Putting Signal through channel with H(z)=1-0.3z^-1
B=[1 -0.3]; % Numerator polynomial coefficients.
A=[1];      % Denominator polynomial coefficients.
[h,t]=impz(B,A,1/Ts);% Generating impulse response.
h=conv(h,x); % sending info through filter.
y=h(2:2000); % Taking only relevant samples.
y=[y(1) y];
figure(1)
subplot(212);
scatter(real(y), imag(y),'*'); grid;
axis([-2 2 -2 2]);
xlabel('In-Phase');
ylabel('Quadrature');
s=['QPSK Constellation through channel, no noise'];
title(s);

% 4) Adding noise to message after it goes through channel.
wi=wgn(1,length(y),-snr); % in-phase WG Noise at SNR=20dB.
wq=wgn(1,length(y),-snr); % quadrature WG Noise at SNR=20dB.
z=y+wi+j*wq;
figure(2)
scatter(real(z), imag(z),'*'); grid;
axis([-2 2 -2 2]);
xlabel('In-Phase');
ylabel('Quadrature');
s=['QPSK Constellation through channel, SNR='];
s=[s num2str(snr) 'dB'];
title(s);

% 5) Implementing adaptive filter.
u=0.005; % step size.
L=16;    % length of estimator.
wh=zeros(L,length(z)+L); % w^, weights.

% Buffered by 8 zeros so that index=9 corresponds to first data.
I=[zeros(1,L/2) x];    % Training sequence is the clean bits.
Ih=zeros(1,length(I)); % I-hat, the estimate.
It=Ih;    % I-tilda, the decision directed output.

e=zeros(1,length(z)+L); % Error vector is 2016 long.
z=[zeros(1,L) z]; % buffering input by 16 zeros.

S=L/2+500; % Sample to switch from training to decision.
a=[1:L]; % Loop limiter vector avoids negative indices.
a=[a L*ones(1,length(z)-length(a))];

for N=9:length(z)-1,

   for k=1:a(N)
     Ih(N-8)=Ih(N-8)+z(N-k+1)*wh(k,N);
   end
   
   %Decision directed.
   It(N-8)=2*(real(Ih(N-8))>0)-1;
   It(N-8)=It(N-8)+j*(2*(imag(Ih(N-8))>0)-1);
   
   if N>S, % then used decision directed.
     e(N)=It(N-8)-Ih(N-8);
   else
     e(N)=I(N-8)-Ih(N-8);
   end
   
   for k=1:a(N),
       wh(k,N+1)=wh(k,N)+u*conj(z(N-k+1))*e(N);
   end  
   
end    

% Removing buffering zeros.
z=z(find(real(z)~=0 | imag(z)~=0));
It=It(find(real(It)~=0 | imag(It)~=0));
Ih=Ih(find(real(Ih)~=0 | imag(Ih)~=0));
e=e(L:length(e));
wh=wh(:,L:2000);

% plotting MSE error and weights.
figure(3)
plot(1:length(e),10*log10(abs(e).^2),'b',1:length(wh),10*abs(wh),'g');
axis([0 2000 -70 10]);
grid;
s=['DDLMS: L='];
s=[s num2str(L) ];
s=[s ', u='];
s=[s num2str(u) ];
s=[s ', Switch value='];
s=[s num2str(S-L/2) ];
title(s);
s=['Weights*10 (green), MSE (blue, dB), '];
ylabel(s);
xlabel('Sample');

% plotting input and output constellations.
figure(4);
clf(4);
hold on;
scatter(real(z),imag(z),'b*');
scatter(real(Ih(S:length(Ih))),imag(Ih(S:length(Ih))),'r*');
h=scatter(real(It(S:length(It))),imag(It(S:length(It))),'go');
set(h,'LineWidth',2);
grid;
hold off;
axis([-2 2 -2 2]);
s=['DDLMS: L='];
s=[s num2str(L) ];
s=[s ', u='];
s=[s num2str(u) ];
s=[s ', Switch value='];
s=[s num2str(S-L/2) ];
s=[s ': Input (blue), Output (red), '];
title(s);
xlabel('In-Phase');
ylabel('Quadrature');

%--------------------------- End Lab 5 --------------------