Jump to content


Toggle shoutbox Shoutbox Open the Shoutbox in a popup

* Answers to registration questions are: Alternating Current (AC) inventor is "Nicolas Tesla", The unlincensed frequency ranges are "2.4 and 5.8 GHz", and internal bleeding procedure is "Lavage".

* Requests are not permitted in the Shoutbox. Also, requests will be ignored if you have less than 10 useful posts.
@  dhampu : (06 May 2017 - 04:31 PM) :)  :)  :)  :)  :)
@  aliahmad : (27 May 2016 - 07:04 AM) hi guys.i need MDI jade 6
@  Oldschool : (24 November 2015 - 08:54 AM) Hi guys :)
@  dhampu : (15 November 2015 - 03:13 PM) :ph34r:  :ph34r:  :ph34r:  :ph34r:  :ph34r:  B)
@  user_128 : (28 April 2015 - 03:10 AM) Halo . Anyone got Minix operting system by Andrew S. Tanenbaum solution manual ? Any version of Minix solution manual would do But prefer minix 3 as it support x86 n arm . Present studyingt n trying to port minix to raspberry pi 2 , hope u guys could share minix solution manual . Thx in advance n regards.
@  Oldschool : (13 March 2015 - 11:43 AM) RHG :) !!!!
@  nsdfxela : (22 February 2015 - 02:01 PM) hi everyone!
@  kazugawa : (21 February 2015 - 10:46 AM) hi...good afternoon all
@  RHGLazyAzz : (17 February 2015 - 10:29 PM) Hi Tom!!!!!
@  Oldschool : (30 January 2015 - 08:10 PM) Let's get this place active, guys :)
@  moha : (10 December 2014 - 09:40 AM) http://eknowz.com/in...ral/&c=LZTREMHX
@  moha : (10 December 2014 - 09:36 AM) Hi :)  :D
@  Oldschool : (01 December 2014 - 07:28 AM) Hi :)
@  moha : (11 October 2014 - 11:41 AM) http://ezproxy.catsboard.com/
@  ducthanh206 : (24 September 2014 - 01:44 PM) can you give the solution of microelectronics circuit analysis and design 4th edition
@  barr : (23 September 2014 - 09:52 AM) hi there
@  moha : (07 August 2014 - 07:16 AM) Hello :D
@  Oldschool : (04 August 2014 - 12:05 PM) Hi guys :)
@  dhampu : (29 June 2014 - 02:07 AM) B)  B)  B)  B)  B)  B)  :ph34r:  :ph34r:  :ph34r:  :ph34r:  :ph34r:  :ph34r:  :ph34r:
@  moha : (25 May 2014 - 05:40 PM) :)

Photo

Recursive Least Squares (RLS)


  • Please log in to reply
4 replies to this topic

#1 atsil

atsil

    Instructor

  • Researcher
  • 111 posts

Posted 16 October 2009 - 09:30 AM

http://en.wikipedia...._squares_filter

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

Attached Files



#2 atsil

atsil

    Instructor

  • Researcher
  • 111 posts

Posted 16 October 2009 - 10:05 AM

Next step is to implement fast recursive algorithm with linear complexity.
Original FORTRAN implementation was already published in reffered book of S.L.Marple -
Digital spectral analysis: with applications, Prentice-Hall, Inc. Upper Saddle River, NJ, USA -
but the problem for many years was not only in coding it in C ))

The problem was that algorithm was rather unstable with respect to rounding errors... so we need follow to stability analisys issues
Fortunately this problem was solved resently, let we follow this article

Stability Analisys of Fast RLS

So We have after all

// comp.load = 5*ip
// delay line more depth (>=0) for reinitialization
#define DEPTH 100
double omegaN;
int fastrls_in(int &init, double omega, int ip, complex *x, double &pf, complex *af, double &pb, complex *ab, double γ, complex *c)
// delay line x[ip+1] (ip+1+DEPTH if reinitialization used)
// Kalman Gain c[ip+1]
// Forward prediction coeff. af[ip]
// Backward prediction coeff ab[ip]
// flag init=0
{
complex sample, temp, hold, eb, epb, ef, epf;
double save, gamma1, tmp, pf1;
complex eb1;
int k, mk;

if (init<=ip) {
if (init==0) {
omegaN = pow(omega,(ip-1));
gamma=1.0;
pf = x[0].real()*x[0].real() + x[0].imag()*x[0].imag();
init++;
return init;
}
ef=x[0];
if (init!=1) {
for (k=0; k<(init-1); k++)
ef = ef + af[k]*x[k+1];
}
epf = ef/gamma;
pf = omega*pf;
hold = conjg(ef)/pf;
gamma = gamma + (ef.real()*hold.real()-ef.imag()*hold.imag());

temp = x[init];
af[init-1] = -ef/temp;
if (init!=1) {
for (k=init-1; k>0; k--)
c[k]= c[k-1] + hold*af[k-1];
}
c[0] = hold;
hold =-temp/(gamma);
if (init>=ip) { // ==ip
pb = (temp.real()*temp.real() + temp.imag()*temp.imag())/gamma;
for (k=0; k<init; k++)
ab[k] = hold*c[init-1-k];
}
init++;
return init;
}

save = omega*pf;
ef = x[0];
for (k=0; k<ip; k++) {
ef= ef + af[k]*x[k+1];
}
epf = ef/gamma;
hold = conjg(ef)/save;
// pf1 = pf; // pf prev
pf1 = save + epf.real()*ef.real() + epf.imag()*ef.imag();
// pf>0

for (k=ip-1; k>=0; k--) {
temp=c[k];
c[k+1] = temp + hold*af[k];
af[k] = af[k] - epf*temp;
}
c[0] = hold;
save=omega*pb;
// eb = complex(save*c[ip].real(), -save*c[ip].imag());
eb = conjg(c[ip]);

#define nMARPLE_ORIGINAL
#ifdef MARPLE_ORIGINAL
eb = eb*save;
gamma = gamma + ef.real()*hold.real() - ef.imag()*hold.imag();
hold = c[ip];
gamma = gamma - eb.real()*hold.real() + eb.imag()*hold.imag();

#else
tmp = pf/(gamma*omegaN);
// 0.25 0.75(?)
#define mu (0.25)
// eb0 = (eb*tmp); // rf1
// eb = (eb * save); // rf0
temp = eb*(save+(tmp-save)*mu); // (1-mu)*rf0+mu*rf1
eb1 = x[ip]; // rc
for (k=0; k<ip; k++) {
eb1 = eb1 + x[ip-k-1]*ab[k];
}
temp = eb1 - temp; // ksi
eb = eb1 + temp; // propagate ksi

gamma1 = gamma + ef.real()*hold.real() - ef.imag()*hold.imag();
hold = c[ip];


#define OLD
#ifdef OLD
// old
// gamma = gamma1 - (eb*hold).real();
gamma = gamma1 - eb.real()*hold.real() + eb.imag()*hold.imag();
#else
// new
// gamma = gamma1 - (eb*hold).real()*gamma*(omegaN*omega)*pb/pf;
gamma = gamma1 - (eb*hold).real()*omega*pb/tmp;
#endif

#endif
pf = pf1;

// reinitialize
if (gamma<0.99) {
init=0;
for (k=0; k<ip+DEPTH; k++) {
fastrls_in(init, omega, ip, x+ip-k+DEPTH, pf, af, pb, ab, gamma, c);
}
return -1;
}

epb = eb/gamma;
pb = save + epb.real()*eb.real() + epb.imag()*eb.imag();
// pb >0;
for (k=0; k<ip; k++) {
mk = ip-1-k;
temp = c[k] - hold*ab[mk];
c[k] = temp;
ab[mk] = ab[mk] - epb*temp;
}
return init;
}

// ip - order
// omega - forgetting factor
// X - input
// Y - output
// x - input delay line
// W - impulse response Y=W convolved X, vector of order ip
// pf, pb - forward & backward prediction error power
// af, ab - forward & backward prediction coeff.
// c - Kalman Gain

int fastrls(int &init, double omega, int ip, complex X, complex Y, complex *x, double &pf, complex *af, double &pb, complex *ab, double γ, complex *c, complex *W)
{
int k;
complex sample = x[0];
complex er;

for (k=ip+DEPTH; k>0; k--)
x[k] = x[k-1];

x[0] = X;
// prediction part
fastrls_in(init, omega, ip, x, pf, af, pb, ab, gamma, c);
// filtering part
er = Y; // filter
for (k=0; k<ip; k++) {
er = er + W[k]*x[k];
}
// adaptation
for (k=0; k<ip; k++) {
W[k] = W[k] + c[k]*er*gamma;
}
return 0;
}

The complexity of algorithm is only 5*ip for prediction part and 2*ip - for filtering.
Voila! it is almost always stable! So it's quite good for modem equalizer or echocanceller design

This fast RLS implementation with test module also placed in attachment

Attached Files



#3 live_11

live_11

    Assistant

  • Student
  • 1 posts

Posted 23 June 2010 - 12:49 PM

Hi atsil,


I am just wondering if you the RLS program you posted above is in C or C++? I am still new to programming in these two languages and still can distinguish between them

Thank u

#4 atsil

atsil

    Instructor

  • Researcher
  • 111 posts

Posted 11 July 2010 - 10:20 AM

Hi atsil,
I am just wondering if you the RLS program you posted above is in C or C++? I am still new to programming in these two languages and still can distinguish between them

Thank u


This is C++ but for only reason - complex numbers arithmetics
You need to define structure complex and write arithmetic operations (*,+,-,/, conjugate) for that manualy - in order to convert program to C

#5 skywolker

skywolker

    Assistant

  • Student
  • 1 posts

Posted 19 October 2013 - 08:53 PM

Thank you for the RLS code, really saved me from trouble.