[转帖]RLS,Demonstration of Wiener filter,LMS filter,Steep...
% Demonstration of Wiener filter,LMS filter,Steep-descent algorithmclear;
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;
= 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(),title('RLS Filter Input(Siganal + Noise)'),drawnow,hold on
subplot(212),plot(t,output),axis(),title('RLS Filtered Signal'),drawnow,hold on
end 学习下，谢谢:@) :victory:
页:
[1]