%--------------------------- Begin Lab 4 --------------------
% Colin Joye
% ECE 3770-01
% Lab 4, Linear Adaptive Prediction
% due 4/9/01
clear;

% 1) Setup.
n=1:2000; % Sample vector [samples]
Ts=10^-3; % Sampling period [s]
f1=1;     % Message frequency [Hz] 
x=cos(2*pi*f1*n*Ts); % Message 
snr=20;  % SNR ratio

% 2) Plotting the message.
figure(1)
subplot(211);
plot(n,x,'b',[1 length(x)],[0 0],'k'); grid;
ylabel('Message (blue)');
title('2) Message and Noise');

% 3) Adding noise to message.
w=wgn(1,length(x),-snr); % WG Noise at SNR=-20dB.
y=x+w;

% 4) Plotting noisy message.
figure(1)
subplot(212);
plot(n,y,'b',[1 length(x)],[0 0],'k'); grid;
ylabel('Message + Noise');
xlabel('Sample');

% 5) Implementing adaptive filter.
u=0.1; % step size.
p=8;   % tap number.
wh=zeros(p,length(x)); % w^, the approximated weights.
yh=zeros(1,length(x)); % y^, the filtered output.
e=yh;                  % the error vector.
a=[1 1:8];
a=[a 8*ones(1,length(x)-length(a))];
for N=2:length(x),
    for k=1:a(N),
       wh(k,N)=wh(k,N-1)+u*y(N-k)*e(N-1);
    end
    yh(N)=y(1+N-a(N):N)*wh(1:a(N),N);
    e(N)=x(N)-yh(N); % the error vector;
end

% 6)
figure(2)
plot(n,e,'b',n,wh,'g',n,y,'y',n,yh,'r'); 
grid;
s=['6) Adaptive filter: p='];
s=[s num2str(p) ];
s=[s ', u='];
s=[s num2str(u) ];
title(s);
ylabel('Weights (green), Output (red), Input (yellow), Error (blue)');
xlabel('Sample');

% 7)
A=1;
figure(3);
N=[1 50 2000];
[h1,f,s] = freqz(wh(:,N(1)),A,512,1/Ts);
[h2,f,s] = freqz(wh(:,N(2)),A,512,1/Ts);
[h3,f,s] = freqz(wh(:,N(3)),A,512,1/Ts);
s.plot   = 'both';    % Plot the magnitude and phase.
s.xunits = 'hz';      % Label the frequency units correctly.
s.yunits = 'linear';  % Plot the magnitude linearly.
h=[h1 h2 h3];
freqzplot(h,f,s);
s=['7) Plotting Filter Response at n=1 (blue), n=50 (green)'];
s=[s ', and n=2000 (red)'];
title(s);

%--------------------------8-----------------
% 8)
x=cos(2*pi*f1*n*Ts)+cos(2*pi*2*f1*n*Ts); %New Message 
snr=20;  % SNR ratio

% 8a) Plotting the message.
figure(4)
subplot(211);
plot(n,x,'b',[1 length(x)],[0 0],'k'); grid;
ylabel('Message (blue)');
title('8a) Message and Noise');

% 8b) Adding noise to message.
w=wgn(1,length(x),-snr); % WG Noise at SNR=-20dB.
y=x+w;

% 8c) Plotting noisy message.
figure(4)
subplot(212);
plot(n,y,'b',[1 length(x)],[0 0],'k'); grid;
ylabel('Message + Noise');
xlabel('Sample');

% 8d) Implementing adaptive filter.
u=0.1; % step size.
p=8;   % tap number.
wh=zeros(p,length(x)); % w^, the approximated weights.
yh=zeros(1,length(x)); % y^, the filtered output.
e=yh;                  % the error vector.
a=[1 1:8];
a=[a 8*ones(1,length(x)-length(a))];
for N=2:length(x),
    for k=1:a(N),
       wh(k,N)=wh(k,N-1)+u*y(N-k)*e(N-1);
    end
    yh(N)=y(1+N-a(N):N)*wh(1:a(N),N);
    e(N)=x(N)-yh(N); % the error vector;
end

% 8e)
figure(5)
plot(n,e,'b',n,wh,'g',n,y,'y',n,yh,'r'); 
grid;
s=['8e) Adaptive filter: p='];
s=[s num2str(p) ];
s=[s ', u='];
s=[s num2str(u) ];
title(s);
ylabel('Weights (green), Output (red), Input (yellow), Error (blue)');
xlabel('Sample');

% 8f)
A=1;
figure(6);
N=[1 50 2000];
[h1,f,s] = freqz(wh(:,N(1)),A,512,1/Ts);
[h2,f,s] = freqz(wh(:,N(2)),A,512,1/Ts);
[h3,f,s] = freqz(wh(:,N(3)),A,512,1/Ts);
s.plot   = 'both';    % Plot the magnitude and phase.
s.xunits = 'hz';      % Label the frequency units correctly.
s.yunits = 'linear';  % Plot the magnitude linearly.
h=[h1 h2 h3];
freqzplot(h,f,s);
s=['8f) Plotting Filter Response at n=1 (blue), n=50 (green)'];
s=[s ', and n=2000 (red), u='];
s=[s num2str(u)];
title(s);

%-------------------- 9a ----------------
% 9a)
x=cos(2*pi*f1*n*Ts)+cos(2*pi*2*f1*n*Ts); %New Message 
snr=20;  % SNR ratio

% 9a2) Adding noise to message.
w=wgn(1,length(x),-snr); % WG Noise at SNR=-20dB.
y=x+w;

% 9a4) Implementing adaptive filter.
u=0.01; % step size.
p=8;   % tap number.
wh=zeros(p,length(x)); % w^, the approximated weights.
yh=zeros(1,length(x)); % y^, the filtered output.
e=yh;                  % the error vector.
a=[1 1:8];
a=[a 8*ones(1,length(x)-length(a))];
for N=2:length(x),
    for k=1:a(N),
       wh(k,N)=wh(k,N-1)+u*y(N-k)*e(N-1);
    end
    yh(N)=y(1+N-a(N):N)*wh(1:a(N),N);
    e(N)=x(N)-yh(N); % the error vector;
end

% 9a5)
figure(7)
plot(n,e,'b',n,wh,'g',n,y,'y',n,yh,'r'); 
grid;
s=['9a5) Adaptive filter: p='];
s=[s num2str(p) ];
s=[s ', u='];
s=[s num2str(u) ];
title(s);
ylabel('Weights (green), Output (red), Input (yellow), Error (blue)');
xlabel('Sample');

% 9a6)
A=1;
figure(8);
N=[1 50 2000];
[h1,f,s] = freqz(wh(:,N(1)),A,512,1/Ts);
[h2,f,s] = freqz(wh(:,N(2)),A,512,1/Ts);
[h3,f,s] = freqz(wh(:,N(3)),A,512,1/Ts);
s.plot   = 'both';    % Plot the magnitude and phase.
s.xunits = 'hz';      % Label the frequency units correctly.
s.yunits = 'linear';  % Plot the magnitude linearly.
h=[h1 h2 h3];
freqzplot(h,f,s);
s=['9a6) Plotting Filter Response at n=1 (blue), n=50 (green)'];
s=[s ', and n=2000 (red), u='];
s=[s num2str(u)];
title(s);

%-------------------- 9b ----------------
% 9b)
x=cos(2*pi*f1*n*Ts)+cos(2*pi*2*f1*n*Ts); %New Message 
snr=20;  % SNR ratio

% 9b2) Adding noise to message.
w=wgn(1,length(x),-snr); % WG Noise at SNR=-20dB.
y=x+w;

% 9b4) Implementing adaptive filter.
u=0.001; % step size.
p=8;   % tap number.
wh=zeros(p,length(x)); % w^, the approximated weights.
yh=zeros(1,length(x)); % y^, the filtered output.
e=yh;                  % the error vector.
a=[1 1:8];
a=[a 8*ones(1,length(x)-length(a))];
for N=2:length(x),
    for k=1:a(N),
       wh(k,N)=wh(k,N-1)+u*y(N-k)*e(N-1);
    end
    yh(N)=y(1+N-a(N):N)*wh(1:a(N),N);
    e(N)=x(N)-yh(N); % the error vector;
end

% 9b5)
figure(9)
plot(n,e,'b',n,wh,'g',n,y,'y',n,yh,'r'); 
grid;
s=['9b5) Adaptive filter: p='];
s=[s num2str(p) ];
s=[s ', u='];
s=[s num2str(u) ];
title(s);
ylabel('Weights (green), Output (red), Input (yellow), Error (blue)');
xlabel('Sample');

% 9b6)
A=1;
figure(10);
N=[1 50 2000];
[h1,f,s] = freqz(wh(:,N(1)),A,512,1/Ts);
[h2,f,s] = freqz(wh(:,N(2)),A,512,1/Ts);
[h3,f,s] = freqz(wh(:,N(3)),A,512,1/Ts);
s.plot   = 'both';    % Plot the magnitude and phase.
s.xunits = 'hz';      % Label the frequency units correctly.
s.yunits = 'linear';  % Plot the magnitude linearly.
h=[h1 h2 h3];
freqzplot(h,f,s);
s=['9b6) Plotting Filter Response at n=1 (blue), n=50 (green)'];
s=[s ', and n=2000 (red), u='];
s=[s num2str(u)];
title(s);
%--------------------------- End Lab 4 --------------------