clear all
close all
show=3;
%lms :1
%nlms :2
%3rd algo :3
p = 4;
sat = 10;
PC=2;
expand = 1;
f = 1;
dt = 0.1;
t = 0 : dt : 500;
amp = 1;
% Initializations
INPUTVEC = 0;
mu = 1e-3;
% Initialize a global variablea
% Convergence factor
len = 8;%filter order
% lin = 200;% input length
lin = length(t);
% xin = randn(lin,1) - 0.5;%input variable
% chaotic input generator
a(1 : p) = randn(1, p);
m(1: p) = randn(1,p);
kk = p + 1;
for kk = p + 1 : length(t)
m(kk) = expand * a * m(kk-1 : -1 : kk - p)';
if m(kk) > sat
reminder = rem(m(kk), sat);
m(kk) = -sat + reminder;
% m(kk) = sat;
elseif m(kk) < -sat
reminder = rem(abs(m(kk)), sat);
m(kk) = sat - reminder;
% m(kk) = -sat;
end
% end
end
% plot(m)
% break
xin= m;
%xin = amp * cos(2 * pi * f * t);
xin = xin';
%xin=amp*sin(2*pi*f*t);
mse = 0;
% unsys
hu = zeros(len, 1);%unknown s/m response
hu(end) = 4;
hu(1) = -4;
x = xin(1:8)';%taking 1 to 8
w = zeros(len, 1);%weight of adaptive filter
% Calculate first iteration
d = x * hu;% desired out put is the conv of input and impulse response of unknwm s/m
y = w' * x';% out put of adaptive filter
dc(1) = d;
yc(1) = y;
e = d - y;
mse(1) = e^2;
w = w + mu * e * x';%next iteration value
if show==1
for j = 2 : lin - len %lms algorithm
x = xin(j : j + len - 1)';
d(j) = x * hu;
y = w' * x';
e = d(j) - y;
mse(j) = e^2;
w = w + mu * e * x';
dc(j) = d(j);
yc(j) = y;
% mu adaptation
% if mse(j) > mse(j-1)
% mu = 10 * mu * rand * 2;
% else
% mu = (mu / 20) * rand * 2;
% end
mu_p(j) = mu;
end
f = figure('Position',[20 20 1240 700]);
set(f,'name','LMS Algorithm')
subplot(3,1,1)
plot(t(1:length(dc)), mse)
% hold on
grid on
title ('Plot of Mean Square Error vs. Iteration Number')
ylabel('MSE')
xlabel('time')
subplot(3,1,2)
plot(t(1:length(dc)), dc);
title('desired output')
ylabel('amplitude')
xlabel('time')
hold on
subplot(3,1,3)
plot(t(1:length(yc)), yc, 'r');
hold on
title('predicted output')
ylabel('amplitude')
xlabel('time')
hold on;
% figure(3)
% title('mu var plot')
% plot(mu_p)
% measure for convergence
%conv_measure = sum( mse .* t(1:length(mse))) / ( length(mse) * t(length(mse)) * mean(mse))
mse1(1)=mse(1);
zh = figure('Position',[20 20 1240 700]);
set(zh,'name','Adapting Phase of LMS Algorithm')
dck=zeros(1,20);
dck=dc(1:20);
yck=zeros(1,20);
yck=yc(1:20);
subplot(2,1,1)
plot(1:20, dck,'g');
ylabel('amplitude')
xlabel('time')
title('Desired Signal')
hold on
subplot(2,1,2)
plot(1:20, yck, 'r');
ylabel('amplitude')
xlabel('time')
title('Actual Signal')
hold on;
end
if show==2
for j = 2 : lin - len
x = xin(j : j + len - 1)';
d(j) = x * hu;
y = w' * x';
e = d(j) - y;
mse1(j) = e^2;
w = w + mu * e * x'/norm(x);
dc(j) = d(j);
yc(j) = y;
mu_p1(j) = mu;
end
g = figure('Position',[20 20 1240 700]);
set(g,'name','NLMS Algorithm')
subplot(3,1,1)
plot(t(1:length(dc)), mse1)
% hold on
grid on
title ('Plot of Mean Square Error vs. Iteration Number')
ylabel('MSE')
xlabel('time')
subplot(3,1,2)
plot(t(1:length(dc)), dc);
title('desired output')
ylabel('amplitude')
xlabel('time')
hold on
subplot(3,1,3)
plot(t(1:length(yc)), yc, 'r');
hold on
title('predicted output')
ylabel('amplitude')
xlabel('time')
hold on
zh = figure('Position',[20 20 1240 700]);
set(zh,'name','Adapting Phase of LMS Algorithm')
dck=zeros(1,20);
dck=dc(1:20);
yck=zeros(1,20);
yck=yc(1:20);
subplot(2,1,1)
plot(1:20, dck,'g');
ylabel('amplitude')
xlabel('time')
title('Desired Signal')
hold on
subplot(2,1,2)
plot(1:20, yck, 'r');
ylabel('amplitude')
xlabel('time')
title('Actual Signal')
hold on;
end
mse2(1)=mse(1);
if show==3
N=2000;
%inp = randn(N,1);
n = randn(lin,1); % generate a random matrix
[b,a] = butter(2,0.25); % design a low pass iir butterworth filter with cut off frequency .25 and order 2
Gz = tf(b,a,-1); % calculate transfer function with b num a denom -1 as sapling freq
h=[b -a(2:length(a))];
inporder=3;
outorder=2;
sysorder = inporder + outorder ;
y = lsim(Gz,xin); % generate a LTI system output, where the system is defined by the tranfer function Gz and input xin.
%add some noise
n = n * std(y)/(15*std(n));
d = y + n;
totallength=size(d,1);
%Take only 50 points for training ( N - inporder 47 = 50 - 3 )
N=50 ;
%begin of the algorithm
%forgetting factor
lamda = 0.999 ;
%initial P matrix
delta = 1e2 ;
P = delta * eye (sysorder ) ;
w = zeros ( sysorder , 1 ) ;
for n = inporder : N
%u(n),u(n-1),u(n-2)
u = xin(n:-1:n-inporder+1) ;
%d(n-1),d(n-2)
outp= d(n-1:-1:n-outorder) ;
u=[u ; outp];
phi = u' * P ;
k = phi'/(lamda + phi * u );
y(n)=w' * u;
e(n) = d(n) - y(n) ;
w = w + k * e(n) ;
P = ( P - k * phi ) / lamda ;
% Just for plotting
Recordedw(1:sysorder,n)=w;
end
%check of results
for n = N+1 : totallength
%u(n),u(n-1),u(n-2)
u = xin(n:-1:n-inporder+1) ;
%d(n-1),d(n-2)
outp= d(n-1:-1:n-outorder) ;
u=[u ; outp];
y(n) = w' * u ;
e(n) = d(n) - y(n) ;
end
figure()
hold on
subplot(3,1,1)
plot(1:length(d),d)
subplot(3,1,2)
plot(1:length(y),y,'r');
subplot(3,1,3)
plot(1:length(e),e,'r');
title('System output') ;
xlabel('Samples')
ylabel('True and estimated output')
zh = figure('Position',[20 20 1240 700]);
set(zh,'name','Adapting Phase of RLS Algorithm')
dck=zeros(1,20);
dck=d(1:20);
yck=zeros(1,20);
yck=y(1:20);
subplot(2,1,1)
plot(1:20, dck,'g');
ylabel('amplitude')
xlabel('time')
title('Desired Signal')
hold on
subplot(2,1,2)
plot(1:20, yck, 'r');
ylabel('amplitude')
xlabel('time')
title('Actual Signal')
hold on;
end