- % 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
|
