RLS algorithm is one of two algorithms usually used for adaptive linear system identification.

There are a lot of scientific articles reffered to RLS. But You hardly ever find C/C++ implementation

in The Net, especially while speaking about most general complex implementation.

I can fill the gap, as I have already did it.

First of all we need implementation of class for complex datatype if we haven't yet. Well, I place it in attachment.

It's not something special. Then RLS algorithm will be

#include "complex.h"

#include "rnd.h"

#include <math.h>

/*

complex conjg(complex x)

{

complex y(x.real(), -x.imag());

return y;

}

*/

#define SYMMETRIC

#define Ord (2)

complex H[Ord];

complex G[Ord];

#ifndef SYMMETRIC

complex Ri[Ord*Ord]; // (Ri)tr = (Ri)*

#else

complex Ri[Ord*(Ord+1)/2];

#endif

complex X[Ord];

complex Wrk[Ord];

void crls_init(int n, complex *X, complex *G, complex *H, complex *Ri, double sigma)

{

double diag=100.0/sigma;

int ind=0;

int i, j;

for (i=0; i<n; i++) {

H[i] = 0.0;

G[i] = 0.0;

X[i] = 0.0;

#ifndef SYMMETRIC

for (j=0; j<n; j++) {

if (i==j) Ri[ind++] = diag;

else Ri[ind++] = 0.0;

}

#else

Ri[ind++] = diag;

for (j=i+1; j<n; j++) {

Ri[ind++] = 0;

}

#endif

}

}

// complixity = 2*N*N + 3*N, or 3*N*N/2 + 3*N for SIMMETRIC

int crls(int n, double omega, complex x, complex y, complex *X, complex *G, complex *Ri0, complex *h)

{

// n - filter length

// omega - forgetting factor

// x - filter input

// y - filter output

// X - input delay line

// G - Kalman gain, vector of length n

// Ri - inverse of correlation matrix n x n, might be simmetrical packed

// h - impulse response, vector of order n

int i, j, ind=0;

double ssum;

complex csum;

// input delay line shift

for (i=n-1; i>0; i--) {

X[i] = X[i-1];

}

X[0] = x;

// adaptation

// calculate Kalman Gain

// ind=0

for (i=0; i<n; i++) {

Wrk[i] =complex(0.,0.);

#ifndef SYMMETRIC

for (j=0; j<n; j++)

Wrk[i] = Wrk[i]+ Ri[ind++]*conjg(X[j]);

#else

ind = i;

for (j=0; j<i; j++) {

Wrk[i] = Wrk[i]+ conjg(Ri[ind]*X[j]);

ind += (n-1-j);

}

for (j=i; j<n; j++)

Wrk[i] = Wrk[i]+ Ri[ind++]*conjg(X[j]);

#endif

}

csum = complex(omega, 0.);

for (i=0; i<n; i++) {

csum = csum + X[i]*Wrk[i];

}

csum = complex(1.0,0.0)/csum;

for (i=0; i<n; i++) {

G[i] = Wrk[i]*csum.re;

}

ssum = 1./omega;

ind =0;

// calculate inverse correlation

for (i=0; i<n; i++) {

#ifndef SYMMETRIC

for (j=0; j<n; j++) {

#else

Ri[ind] = (Ri[ind] - G[i]* conjg(Wrk[i])).re*ssum;

ind++;

for (j=i+1; j<n; j++) {

#endif

Ri[ind] = (Ri[ind] - G[i]* conjg(Wrk[j]))*ssum;

ind++;

}

}

// filtering

// error scalar

csum = y;

for (i=0; i<n; i++) {

csum = csum - h[i]*X[i];

}

// impulse response update

for (i=0; i<n; i++) {

h[i] = h[i] + G[i]*csum;

}

return 0;

}

SYMMETRIC implementation is more efficient and more stable as it use property of conjugate transpose ( Hermitian transpose, or adjoint matrix ) of Ri matrix. But algorithm is not fast yet...

Formulas used are from this exelent book:

http://portal.acm.or...FTOKEN=66044442