Adaptive filters

function LMSADF

%Program to illustrate adaptive filtering using the LMS algorithms

% X delayed input data vector

% Y measured signal

% W coefficient vector

% E enhanced signal

N=30; % filter length

M=0; % delay

w0=1; % initial value for adaptive filter coefficients

SF=2048; % factor for reducing the data samples - 11 bit ADC assumed

mu=0.04;

X = zeros(N,1);

delay = zeros(1,M+1);

W = w0*ones(N,1);

in = fopen('ADF.dat','r'); %read input data from specified data file

Y = fscanf(in,'%g',inf)/SF;

fclose(in);

if w0==0

sf = SF; % scaling factor for display

else

sf = SF/N/w0;

end

for i=1:length(Y)

if M>0

delay(2:M+1) = delay(1:M); % shift data for delay

end

delay(1) = Y(i);

X(2:N) = X(1:N-1); % update buffer

X(1) = delay(M+1);

E(i) = Y(i)-W'*X; % the enhanced signal

W = W + 2*mu*E(i)*X; % update the weights

end

subplot(2,1,1),plot(1:length(Y),Y*SF); title('Input Signal');

subplot(2,1,2),plot(1:length(E),E*sf); title('Enhanced Signal');

====================================================

function UDUADF

% program to illustrate adaptive filtering using

% the RLS algorithm via the UDU factorization

% X delayed input data vector

% Y measured signal

% W coefficient vector

% E enhanced signal

clear all;

N = 30; % filter length

M = 1; % delay

npt = N*(N+1)/2;

SF = 2048; % 12-bit ADC scaling

p0 = 0.05;

w0 = 1;

gamma = 0.98;

RemoveMean = 0; % 1 - remove the mean from the data, 0 - otherwise

delay = zeros(1,M);

U=zeros(1,npt);

U(1)=p0;

W = w0*ones(N,1);

X = zeros(N,1);

for i=1:N-1

ik=(i*(i+1)-2)/2+1;

U(ik)=p0;

end

if w0==0

sf = SF; % scaling factor for display

else

sf = SF/N/w0;

end

in = fopen('ADF.dat','r'); %read input data from specified data file

Y = fscanf(in,'%g',inf)/SF;

fclose(in);

if RemoveMean % remove the mean from the data if required

Y = Y - sum(Y)/length(Y);

end

for i=1:length(Y)

if M>0

delay(2:M+1) = delay(1:M); % shift input data in delay registers

end

delay(1) = Y(i);

X(2:N) = X(1:N-1); % update buffer

X(1) = delay(M+1);

E(i) = Y(i) - X'*W; % the enhanced signal

W = uduflt(W,X,U,E(i),gamma ,N);

end

subplot(2,1,1),plot(1:length(Y),Y*SF); title('Input Signal');

subplot(2,1,2),plot(1:length(E),E*sf); title('Enhanced Signal');

==========================================

function w=uduflt(w,x,u,ek,gamma,N)

% udu algorithm - a numerically stable form of

% the recursive least squares algorithm

%

% inputs:

% x() input vector

% dn latest input data value

% w() coefficient vector

% u() vector containing elements of U and D

%

% outputs:

% en error signal

% yn digital filter output

% w() updated coefficient vector

% u() updated elements of U and D

%

sf = 1/gamma;

m=1; % update the UD elements

v=zeros(1,N);

v(1)=x(1);

for j=2:N

v(j)=x(j);

for k=1:j-1

m=m+1;

v(j)=v(j)+u(m)*x(k);

end

m=m+1;

b(j)=u(m)*v(j);

end

b(1)=u(1)*x(1);

alpha=gamma+b(1)*v(1);

delta=1/alpha;

u(1)=u(1)*delta;

m=1;

for j=2:N

beta1=alpha;

alpha=alpha+b(j)*v(j);

p=-v(j)*delta;

delta=1/alpha;

for k=1:j-1

m=m+1;

beta=u(m);

u(m)=beta+b(k)*p;

b(k)=b(k)+b(j)*beta;

end

m=m+1;

u(m)=u(m)*beta1*delta*sf;

end

perr=ek/alpha;

for j=1:N % update the weights

w(j)=w(j)+b(j)*perr;

end

============================================

function SQRTADF

% program to illustrate adaptive filtering using

% the square root RLS algorithm

% X delayed input data vector

% Y measured signal

% W coefficient vector

% E enhanced signal

N = 30; % filter length

M = 1; % delay

npt = N*(N+1)/2;

SF = 2048; % 12-bit ADC scaling

p0 = 0.05;

w0 = 1;

gamma = 0.98;

RemoveMean = 0; % 1 - remove the mean from the data, 0 - otherwise

delay = zeros(1,M);

W = w0*ones(N,1);

X = zeros(N,1);

S = zeros(1,npt);

S(1)=p0;

for i=1:N-1

ik=(i*(i+1)-2)/2+1;

S(ik)=p0;

end

if w0==0

sf = SF; % scaling factor for display

else

sf = SF/N/w0;

end

in = fopen('ADF.dat','r'); %read input data from specified data file

Y = fscanf(in,'%g',inf)/SF;

fclose(in);

if RemoveMean % remove the mean from the data if required

Y = Y - sum(Y)/length(Y);

end

for i=1:length(Y)

if M>0

delay(2:M+1) = delay(1:M); % shift input data in delay registers

end

delay(1) = Y(i);

X(2:N) = X(1:N-1); % update buffer

X(1) = delay(M+1);

E(i) = Y(i) - X'*W; % the enhanced signal

W = sqrtflt(W,X,E(i),S,gamma,N);

end

subplot(2,1,1),plot(1:length(Y),Y*SF); title('Input Signal');

subplot(2,1,2),plot(1:length(E),E*sf); title('Enhanced Signal');

==================================================

function w=sqrtflt(w,x,perr,s,gamma,N)

% A simple square root RLS adaptive filter

% For details, see:

% Digital Signal Processing: A Practical Approach

% E C Ifeachor and B W Jervis, Pearson, 2002

forgt=sqrt(gamma);

sig=forgt;

sigsq=forgt*forgt;

ij=1; ji=1;

for j=2:N

fj=0.0;

for i=1:j-1

ji=ji+1;

fj=fj+s(ji)*x(i);

end

a=sig/forgt;

b=fj/sigsq;

sigsq=sigsq+fj*fj;

sig=sqrt(sigsq);

a=a/sig;

g(j)=s(ji)*fj;

s(ji)=a*s(ji);

for i=1:j-1

ij=ij+1;

sqp=s(ij);

s(ij)=a*(sqp-b*g(i));

g(i)=g(i)+sqp*fj;

end

ij=ij+1;

end

w = w + g'*perr/sigsq;

=============================

function RLSadf

% program to illustrate adaptive filtering using

% the RLS algorithm

% X delayed input signal

% Y measured signal

% W coefficient vector

% E enhanced signal

N = 30; % filter length

M = 1; % stages of delay

SF = 2048; % 12-bit ADC scaling

p0 = 0.05;

w0 = 100;

gamma = 0.98;

RemoveMean = 0; % 1 to remove the mean, 0 otherwise

W = w0*ones(N,1); % adaptive filter weights

X = zeros(N,1);

delay = zeros(1,M+1);

P = p0*diag(ones(1,N),0);

if w0==0

sf = SF; % scaling factor for display

else

sf = SF/N/w0;

end

in = fopen('ADF.dat','r'); %read input data from specified data file

Y = fscanf(in,'%g',inf)/SF;

fclose(in);

if RemoveMean % remove the mean from the data if required

Y = Y - sum(Y)/length(Y);

end

for i=1:length(Y)

if M>0

delay(2:M+1) = delay(1:M); % shift input data in delay registers

end

delay(1) = Y(i);

X(2:N) = X(1:N-1); % update buffer

X(1) = delay(M+1);

E(i) = Y(i) - X'*W; % the enhanced signal

G = P*X/(gamma + X'*P*X);

P = (P - G*X'*P)/gamma;

W = W + G*E(i); % update the weights

end

subplot(2,1,1),plot(1:length(Y),Y*SF); title('Input Signal');

subplot(2,1,2),plot(1:length(E),E*sf); title('Enhanced Signal');