Variable step-size NLMS adaptive algorithm

%=================================================
clear all; close all; clc;

k           = 5000;      % number of samples
run         = 100;
tap         = 64;        % number of filter weights
sigma_n2    = 0.01;  	% noise power
sigma_x2    = 1;        % signal power

unknown_system = fircband(tap-1,[0 0.4 0.5 1], [1 1 0 0],[1 0.2], {'w' 'c'});
x = sqrt(sigma_x2)*randn(k,1); % Creating the input signal
n = sqrt(sigma_n2)*randn(1,k);       % noise

d = zeros(k,1);     % desired signal = unknown_system' * x + n
X = zeros(tap,1);
for p = 1:k
	X       =   [x(p,1); X(1:(tap-1),1)];     % input signal (tapped delay line)
	d(p)    =   (unknown_system*X(:,1)) + n(p);       % desired signal
end

xk = [zeros(1,tap-1) x'];
w0 = zeros(tap,1);
weight(:,1) = reshape(w0,length(w0),1);

mu = 1;
alpha = 0.998;
beta = 15;
epsilon = 0.01;
sigma2_e = 0;
sigma2_x = 0;
RexHat = 0;

for i = 1:k;
    Xk(1:tap,1) = xk(i+(tap-1):-1:i);
    Yk(i) = weight(:,i)' * Xk;
    e(i) = d(i) - Yk(i);

      sigma2_e = alpha*sigma2_e + (1-alpha)*e(i)^2;		% Estimate MSE
      sigma2_x = alpha*sigma2_x + (1-alpha)*mean(Xk)^2;	% Estimate input power
      RexHat = alpha*RexHat + (1-alpha)*Xk*e(i);
      sigma2_v = sigma2_e - (1/sigma2_x)*(RexHat'*RexHat);

      mu = alpha*mu + (1-alpha)*sigma2_e/(beta*sigma2_v);

      weight(:,i+1) = weight(:,i) + (mu/(epsilon+Xk'*Xk))*conj(e(i))*Xk;
end

plot(1:length(e),20*log10(e.^2));