[转帖]RLS,Demonstration of Wiener filter,LMS filter,Steep...

happy

1. % Demonstration of Wiener filter,LMS filter,Steep-descent algorithm
   
2. clear;
   
3. clc;
   
4. 
   
5. N = 10000; %----- the length of the observation sequence
   
6. M = 2; %----- the filter length
   
7. 
   
8. v = randn(1,N); %----- white process as the AR excitation
   
9. a = poly(sign(randn(1,M)).*rand(1,M)); %----- coefficient of AR process
   
10. 
    
11. u = filter(1,a,v); %-----the input sequence
    
12. 
    
13. d = v; %----- the desired response
    
14. 
    
15. rf = xcorr(u,M,'biased');
    
16. rv = rf(M+1:2*M+1);
    
17. 
    
18. R = toeplitz(rv); %----- the correlation matrix of the input
    
19. 
    
20. pf = xcorr(d,u,M,'biased');
    
21. 
    
22. pv = pf(M+1:2*M+1).'; %----- the cross-correlation vector between the input and the desired response
    
23. 
    
24. %----- the optimal tap weight vector for Wiener filter-----
    
25. wopt = inv(R) * pv;
    
26. [V,D] = eig(R); %-----selection of a stable step size mu
    
27. lambda_max = max(diag(D));
    
28. mu = 0.9 * 2/lambda_max;
    
29. 
    
30. %----- the steepest descent learning-----
    
31. wsd = randn(M+1,1); %-----initial weight vector for steepest descent
    
32. total_iteration_number = 100; %-----total iteration number
    
33. 
    
34. for i=1:total_iteration_number
    
35. wsd = wsd + mu * (pv - R*wsd);
    
36. end
    
37. 
    
38. %----- the LMS learning-----
    
39. wlms = randn(M+1,1); %-----initial weight vector for LMS
    
40. uv = zeros(M+1,1); %-----initial input vector
    
41. 
    
42. mu = 0.1*2/lambda_max %-----step size mu for LMS
    
43. 
    
44. for n=1:N;
    
45. uv(2:M+1) = uv(1:M);
    
46. uv(1) = u(n);
    
47. 
    
48. y = wlms' * uv;
    
49. e = d(n) - y;
    
50. wlms = wlms + mu * uv * conj(e);
    
51. end;

复制代码

happy

1. % RLS Adaptive Noise cancellation
   
2. 
   
3. %-----Filter Parameters-----;
   
4. M = 20;
   
5. delta = 1; %%%---------------diffenrence
   
6. lamda = 0.99; %%%---------------diffenrence
   
7. mu = 0.05;
   
8. e_max = 400; %-----maximum of epochs
   
9. 
   
10. %-----Contants-----
    
11. pi = 3.14;
    
12. Fs = 0.01; %-----signal frequency
    
13. Fn = 0.05; %-----noise frequency
    
14. 
    
15. %-----Initialize-----
    
16. w = zeros(M,1); %%%---------------diffenrence
    
17. d = zeros(M,1); %%%---------------diffenrence
    
18. u = zeros(M,1); %%%---------------diffenrence
    
19. P = eye(M)/delta; %%%---------------diffenrence
    
20. 
    
21. %-----Generate desired signal and input(signal+noise)-----
    
22. for t=1:M-1
    
23. d(t) = sin(2*pi*Fs*t);
    
24. u(t) = d(t) + 0.5*sin(2*pi*Fn*t) + 0.09*randn;
    
25. end
    
26. t = M;
    
27. epoch = 0;
    
28. 
    
29. while epoch<e_max %-----generate new input
    
30. input = sin(2*pi*Fs*t);
    
31. for i=2:M %-----shift new input into array
    
32. d(M-i+2) = d(M-i+1);
    
33. u(M-i+2) = u(M-i+1);
    
34. end
    
35. d(1) = input;
    
36. u(1) = input +0.5*sin(2*pi*Fn*t)+0.09*randn; %-----add undesired freq & random noise
    
37. 
    
38. 
    
39. output = w'*u; %-----compute filter output
    
40. 
    
41. 
    
42. %-----RLS algorithm-----
    
43. k = (P*u)/(lamda + u'*P*u); %-----compute error-----
    
44. 
    
45. E = d(1) - w'*u;
    
46. w = w + k*E;
    
47. 
    
48. P = (P/lamda) - (k*u'*P/lamda);
    
49. 
    
50. int(t-M+1) = u(1);
    
51. out(t-M+1) = output;
    
52. err(t-M+1) = E;
    
53. 
    
54. t = t+1;
    
55. epoch = epoch + 1;
    
56. 
    
57. %-----plot noise and filtered signal-----
    
58. figure(2);
    
59. subplot(211),plot(t,u(1)),axis([0 e_max -2.5 2.5]),title('RLS Filter Input(Siganal + Noise)'),drawnow,hold on
    
60. subplot(212),plot(t,output),axis([0 e_max -2.5 2.5]),title('RLS Filtered Signal'),drawnow,hold on
    
61. end

复制代码

yyfei11

学习下，谢谢:@) :victory: