新建 文本文档 (4).txt

RLS,Demonstration of Wiener filter,LMS filter,Steep-descent algorithm

RLS,Demonstration of Wiener filter,LMS filter,Steep-descent algorithm
% Demonstration of Wiener filter,LMS filter,Steep-descent algorithm
clear;
clc;

N = 10000; %----- the length of the observation sequence
M = 2; %----- the filter length

v = randn(1,N); %----- white process as the AR excitation
a = poly(sign(randn(1,M)).*rand(1,M)); %----- coefficient of AR process

u = filter(1,a,v); %-----the input sequence

d = v; %----- the desired response

rf = xcorr(u,M,'biased');
rv = rf(M 1:2*M 1);

R = toeplitz(rv); %----- the correlation matrix of the input

pf = xcorr(d,u,M,'biased');

pv = pf(M 1:2*M 1).'; %----- the cross-correlation vector between the input and the desired response

%----- the optimal tap weight vector for Wiener filter-----
wopt = inv(R) * pv;
[V,D] = eig(R); %-----selection of a stable step size mu
lambda_max = max(diag(D));
mu = 0.9 * 2/lambda_max;

%----- the steepest descent learning-----
wsd = randn(M 1,1); %-----initial weight vector for steepest descent
total_iteration_number = 100; %-----total iteration number

for i=1:total_iteration_number
wsd = wsd mu * (pv - R*wsd);
end

%----- the LMS learning-----
wlms = randn(M 1,1); %-----initial weight vector for LMS
uv = zeros(M 1,1); %-----initial input vector

mu = 0.1*2/lambda_max %-----step size mu for LMS

for n=1:N;
uv(2:M 1) = uv(1:M);
uv(1) = u(n);

y = wlms' * uv;
e = d(n) - y;
wlms = wlms mu * uv * conj(e);
end;

 

% RLS Adaptive Noise cancellation

%-----Filter Parameters-----;
M = 20;
delta = 1; %%%---------------diffenrence
lamda = 0.99; %%%---------------diffenrence
mu = 0.05;
e_max = 400; %-----maximum of epochs

%-----Contants-----
pi = 3.14;
Fs = 0.01; %-----signal frequency
Fn = 0.05; %-----noise frequency

%-----Initialize-----
w = zeros(M,1); %%%---------------diffenrence
d = zeros(M,1); %%%---------------diffenrence
u = zeros(M,1); %%%---------------diffenrence
P = eye(M)/delta; %%%---------------diffenrence

%-----Generate desired signal and input(signal noise)-----
for t=1:M-1
d(t) = sin(2*pi*Fs*t);
u(t) = d(t) 0.5*sin(2*pi*Fn*t) 0.09*randn;
end
t = M;
epoch = 0;

while epoch<e_max %-----generate new input
input = sin(2*pi*Fs*t);
for i=2:M %-----shift new input into array
d(M-i 2) = d(M-i 1);
u(M-i 2) = u(M-i 1);
end
d(1) = input;
u(1) = input 0.5*sin(2*pi*Fn*t) 0.09*randn; %-----add undesired freq & random noise


output = w'*u; %-----compute filter output


%-----RLS algorithm-----
k = (P*u)/(lamda u'*P*u); %-----compute error-----

E = d(1) - w'*u;
w = w k*E;

P = (P/lamda) - (k*u'*P/lamda);

int(t-M 1) = u(1);
out(t-M 1) = output;
err(t-M 1) = E;

t = t 1; 
epoch = epoch 1;

%-----plot noise and filtered signal-----
figure(2);
subplot(211),plot(t,u(1)),axis([0 e_max -2.5 2.5]),title('RLS Filter Input(Siganal Noise)'),drawnow,hold on
subplot(212),plot(t,output),axis([0 e_max -2.5 2.5]),title('RLS Filtered Signal'),drawnow,hold on
end