Matlab - Algoritmo RLS

   
Vista:

Algoritmo RLS

Publicado por Aitor (9 intervenciones) el 29/04/2015 12:48:57
Buenos días/tardes a todos,

Necesito programar el algoritmo RLS para mi proyecto fin de carrera, pero tras mucho estrellarme no he logrado conseguirlo. Mi código es éste:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
clear all
 clf
 N=6; % Number of antennas
 SNR=20; % SIgnal to noise ratio
 mu=.01; % LMS step
 niter=50000; % Number of iterations
 Npro=10; % Number of promedies
 
 lambda = 0.99; % Forgetting factor RLS
 delta = 0.001; % Initialization of P(0)
 
 % Scenario
 a=25/180*pi; % Desired signal
 b1=125/180*pi; % Interferring signal 1
 b2=-90/180*pi; % Interferring signal 2
 g=1; % Gain of our desired signal
 s2n=10.^(-SNR/10);
 
 % Definition of signals
 sa=1/sqrt(N)*exp(j*pi*(0:N-1)*sin(a))';
 sb1=1/sqrt(N)*exp(j*pi*(0:N-1)*sin(b1))';
 sb2=1/sqrt(N)*exp(j*pi*(0:N-1)*sin(b2))';

 % Correlation matrices
 Rs=sa*sa';
 Rint=sb1*sb1'+sb2*sb2'+s2n*eye(N);
 Rx=Rs+Rint;
 
 % Projection and blocking matrices
 A=eye(N);
 A(:,1)=sa;
 
 % Gram-Schmidt algorithm
 Aort=zeros(size(A));
 Aort(:,1)=A(:,1)/sqrt(real((A(:,1)'*A(:,1))));
 for k1=2:N
    for k2=1:(k1-1)
     Aort(:,k1)=Aort(:,k1)+(Aort(:,k2)'*A(:,k1))*Aort(:,k2);
    end
    Aort(:,k1)=A(:,k1)-Aort(:,k1);
    Aort(:,k1)=Aort(:,k1)/sqrt(real((Aort(:,k1)'*Aort(:,k1))));
 end
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 C=Aort(:,1); % Restrictions subspace
 B=Aort(:,2:end); % Restrictions' orthogonal subspace
 P_C=C*inv(C'*C)*C';
 P_C_ort=eye(N)-P_C;
 
 % Optimal solution
 wopt=inv(Rx)*C*inv(C'*inv(Rx)*C)*g;
 wq=sa;

 e_gslc_LMS=zeros(Npro,niter); % Error sequence LMS
 e_gslc_RLS=zeros(Npro,niter); % Error sequence RLS

 % GSLC algorithm
 for kk=1:Npro
    Xa = kron(sa,ones(1,niter)).*kron(ones(N,1),(sign(randn(1,niter))+1i*sign((randn(1,niter)))/sqrt(2)));
    Xb1 = kron(sb1,ones(1,niter)).*kron(ones(N,1),(sign(randn(1,niter))+1i*sign((randn(1,niter)))/sqrt(2)));
    Xb2 = kron(sb2,ones(1,niter)).*kron(ones(N,1),(sign(randn(1,niter))+1i*sign((randn(1,niter)))/sqrt(2)));
    X = Xa + Xb1 + Xb2 + (randn(size(Xa)) + 1i*randn(size(Xa)))/sqrt(2)*sqrt(s2n);

    w_gslc_LMS = zeros(N,niter);
    w_gslc_RLS = zeros(N,niter);
    xa = B'*X;
 
    wa_LMS = ones(N-1,niter);
    wa_RLS = ones(N-1,niter);
 
    P = inv(delta)*eye(N-1);
 
    for k = 2:niter
       y_LMS = X(:,k)'*w_gslc_LMS(:,k-1);
       y_RLS = X(:,k)'*w_gslc_RLS(:,k-1);
 
       % LMS
       wa_LMS(:,k) = wa_LMS(:,k-1) + mu*xa(:,k)*(X(:,k)'*wq - xa(:,k)'*wa_LMS(:,k-1));
 
       % RLS
       z = P*xa(:,k);
       alfa = X(:,k)'*wq - xa(:,k)'*wa_RLS(:,k-1);
       ge = z/(lambda + xa(:,k)'*z);
       P = (P - ge*z')/lambda;
       wa_RLS(:,k) = wa_RLS(:,k-1) + alfa*ge; 
       w_gslc_LMS(:,k) = wq - B*wa_LMS(:,k);
       w_gslc_RLS(:,k) = wq - B*wa_RLS(:,k);
    end
    e_gslc_LMS(kk,:) = sum(abs(w_gslc_LMS - kron(wopt,ones(1,niter))).^2,1)/N;
    e_gslc_RLS(kk,:) = sum(abs(w_gslc_RLS - kron(wopt,ones(1,niter))).^2,1)/N;
 end
 
 figure(1)
 plot(10*log10(mean(e_gslc_LMS)),'b'), hold on
 plot(10*log10(mean(e_gslc_RLS)),'r')
 grid
 xlabel('Iterations')
 ylabel('dB')
 title('GSLC algorithm - Learning curve of LMS/RLS algorithms')

Para intentarlo me he basado en la información que aparece en su entrada de Wikipedia:
http://en.wikipedia.org/wiki/RLS_algorithm#Lattice_recursive_least_squares_filter_.28LRLS.29

Por favor, ¡necesito vuestra ayuda!


Un saludo.
Valora esta pregunta
Me gusta: Está pregunta es útil y esta claraNo me gusta: Está pregunta no esta clara o no es útil
0
Responder

Algoritmo RLS

Publicado por Aitor (9 intervenciones) el 03/05/2015 22:21:05
¿Alguien me puede ayudar por favor?

Sigo sin descifrar la respuesta, y verdaderamente lo necesito.

Quedo eternamente agradecido con aquél/lla que lo consiga.

Un saludo.
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar

Algoritmo RLS

Publicado por Aitor (9 intervenciones) el 06/05/2015 13:31:28
¿Hola? :( ¿Hay alguien?
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar

Algoritmo RLS

Publicado por AItor (9 intervenciones) el 11/05/2015 20:52:07
Por favor, ayuda...
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar