Weird behaviour in the nLMS code.
CASTalk.com Forum Index CASTalk.com
Discussion of DSP, FPGA, storage and embedded system.
 
 FAQFAQ   MemberlistMemberlist     RegisterRegister 
 ProfileProfile   Log in to check your private messagesLog in to check your private messages   Log inLog in 
 
Google
 
Web castalk.com
Weird behaviour in the nLMS code.
Goto page 1, 2  Next
 
Post new topic   Reply to topic    CASTalk.com Forum Index -> DSP
Author Message
Jack L.
Guest
Posted: Sat Dec 18, 2004 6:43 pm    Post subject: Weird behaviour in the nLMS code. Reply with quote
Dear group,
I have a problem with my nLMS code written in C. When the Euclidean norm
isn't applied to the code, i.e., it's an ordinary LMS algorithm, the system
destabilises. Once applied, however, the weight vectors (seem to) converge.
I can't believe the system stabilises only when normalisation is applied;
this would suggest the LMS never worked. o_O

Also, if you look at the code below you'll see two different versions of the
Euclidean norm. I got the first version from a former group who had also
been involved in this project, but I've never understood how it evaluates
the Euclidean norm, where's the sqrt? Is it an approximation? I wrote the
second version and I believe this is the correct way of evaluating the
Euclidean norm though it makes the code slow. However, the second version
doesn't seem to make system more stable.

The algorithm is used in removing environmental noise and noise from the
body.

The algorithm has been tested in MatLab's Simulink by implementing it as a
S-function; it seems to be performing worse than the 'lmsdemo' in MatLab.
It's then implemented in an Analog Devices SHARC ADSP-21161N (also in C)
with few modifications.

I admit this is a school project (final) which isn't welcome in the group,
but I've fought with the stability problem for a long time without luck and
you're my last resort. I've found few LMS algorithms here and there on the
Internet and my code seems to be similar to theirs.

All in all, some inputs on the code would be appreciated.

Yours,
Jack

Code:

Global initialisations:
#define TAPS // filter length
float ref_0[TAPS] // delay line
float w_0[TAPS] // weight vectors
float euclidean;

void Process_Sample()
{
int n;
double y, error, norm_my, d;

/* Version 1 of the Euclidean norm BEGIN */
euclidean -= ref_0[TAPS-1]*ref_0[TAPS-1];
euclidean += ref_0[0]*ref_0[0];
/* Version 1 of the Euclidean norm END */

/* Version 2 of the Euclidean norm BEGIN */
euclidean = 0;
for (n = 0; n < TAPS; n++) {
euclidean += ref_0[n] * ref_0[n];
}
euclidean = sqrt(euclidean);
/* Version 2 of the Euclidean norm END */

// Delay-line
for (n=TAPS-1; n>=1; n--) {
ref_0[n] = ref_0[n-1];
}

// Inserting new sample into the delay line
ref_0[0] = u0_[0];

// FIR-filter
y = w_0[0]*ref_0[0];
for (n=1; n<TAPS; n++) {
y += w_0[n]*ref_0[n];
}

error = d - y;

norm_my = stepsize_0 * error / (b + euclidean);
for (n=0; n<TAPS; n++) {
w_0[n] += ref_0[n] * norm_my;
}
}
Back to top
Peter K.
Guest
Posted: Sat Dec 18, 2004 7:03 pm    Post subject: Re: Weird behaviour in the nLMS code. Reply with quote
Jack L. wrote:

Quote:
I have a problem with my nLMS code written in C. When the
Euclidean norm isn't applied to the code, i.e., it's an
ordinary LMS algorithm, the system destabilises. Once
applied, however, the weight vectors (seem to) converge.
I can't believe the system stabilises only when
normalisation is applied; this would suggest the
LMS never worked. o_O

That's part of the point of the NLMS algorithm: it converges over a
much better range of step-sizes.

What value are you using for stepsize_0? Make it smaller and the LMS
should converge quite happily... provided the code does implement the
LMS :-)

I've had a brief look, and it seems to, but the application seems
different from stuff I've done before.

Quote:
Also, if you look at the code below you'll see two
different versions of the Euclidean norm. I got the
first version from a former group who had also been
involved in this project, but I've never understood
how it evaluates the Euclidean norm, where's the
sqrt? Is it an approximation? I wrote the second
version and I believe this is the correct way of
evaluating the Euclidean norm though it makes the
code slow. However, the second version
doesn't seem to make system more stable.

Version 1 is just assuming euclidean is initialized according to your
version 2, then all you have to do is delete the oldest old sample and
include the newest sample).

Ciao,

Peter K.
Back to top
Randy Yates
Guest
Posted: Sat Dec 18, 2004 8:34 pm    Post subject: Re: Weird behaviour in the nLMS code. Reply with quote
"Jack L." <jack_nospam@nospam.dk> writes:
Quote:
[...]
/* Version 1 of the Euclidean norm BEGIN */
euclidean -= ref_0[TAPS-1]*ref_0[TAPS-1];

This should be

euclidean -= ref_0[TAPS]*ref_0[TAPS];
--
% Randy Yates % "The dreamer, the unwoken fool -
%% Fuquay-Varina, NC % in dreams, no pain will kiss the brow..."
%%% 919-577-9882 %
%%%% <yates@ieee.org> % 'Eldorado Overture', *Eldorado*, ELO
http://home.earthlink.net/~yatescr
Back to top
Randy Yates
Guest
Posted: Sat Dec 18, 2004 8:41 pm    Post subject: Re: Weird behaviour in the nLMS code. Reply with quote
The correction in my other post was incomplete. I
think this will do it:

"Jack L." <jack_nospam@nospam.dk> writes:

Quote:
Global initialisations:
#define TAPS // filter length
float ref_0[TAPS+1] // delay line
float w_0[TAPS] // weight vectors
float euclidean;

void Process_Sample()
{
int n;
double y, error, norm_my, d;

/* Version 1 of the Euclidean norm BEGIN */
euclidean -= ref_0[TAPS]*ref_0[TAPS];
euclidean += ref_0[0]*ref_0[0];
/* Version 1 of the Euclidean norm END */

/* Version 2 of the Euclidean norm BEGIN */
euclidean = 0;
for (n = 0; n < TAPS; n++) {
euclidean += ref_0[n] * ref_0[n];
}
euclidean = sqrt(euclidean);
/* Version 2 of the Euclidean norm END */

// Delay-line
for (n=TAPS; n>=1; n--) {
ref_0[n] = ref_0[n-1];
}

// Inserting new sample into the delay line
ref_0[0] = u0_[0];

// FIR-filter
y = w_0[0]*ref_0[0];
for (n=1; n<TAPS; n++) {
y += w_0[n]*ref_0[n];
}

error = d - y;

norm_my = stepsize_0 * error / (b + euclidean);
for (n=0; n<TAPS; n++) {
w_0[n] += ref_0[n] * norm_my;
}
}



--
% Randy Yates % "Watching all the days go by...
%% Fuquay-Varina, NC % Who are you and who am I?"
%%% 919-577-9882 % 'Mission (A World Record)',
%%%% <yates@ieee.org> % *A New World Record*, ELO
http://home.earthlink.net/~yatescr
Back to top
Jack L.
Guest
Posted: Sun Dec 19, 2004 2:57 am    Post subject: Re: Weird behaviour in the nLMS code. Reply with quote
Randy Yates wrote:
Quote:
The correction in my other post was incomplete. I
think this will do it:


Hi Randy, thanks for your correction. However, I don't see what changes they
do to the code. As I see it, if the desired tap length is 600 then I'd
define TAPS to 599 with your correction, while defining TAPS to 600 with my
version?

--
Mvh / Best regards,
Jack, Copenhagen

The email address is for real. :)
Back to top
Randy Yates
Guest
Posted: Sun Dec 19, 2004 7:35 am    Post subject: Re: Weird behaviour in the nLMS code. Reply with quote
"Jack L." <jack_nospam@nospam.dk> writes:

Quote:
Randy Yates wrote:
The correction in my other post was incomplete. I
think this will do it:


Hi Randy, thanks for your correction. However, I don't see what changes they
do to the code. As I see it, if the desired tap length is 600 then I'd
define TAPS to 599 with your correction, while defining TAPS to 600 with my
version?

Hi Jack,

No. Define TAPS to be the number of taps of the filter.

Other than the absence of the square root operation, the only thing
version 2 is doing that is different than version 1 is using a
different structure for computing the averaging filter.

A moving average filter of length N is an FIR filter with impulse
response b = 1/N * ones(1, N). The straight-forward way to compute
it is non-recursively (as in version 1):

y[n] = (1/N) * sum_{j=0}^{N-1} x[n-j].

However, there is a recursive way to compute the filter as well:

y[n] = y[n-1] + (1/N) * (x[n] - x[n-N]).

This can be easily derived by reexpressing the non-recursive form as

y[n-1] = (1/N) * (x[n - 1] + x[n - 2] + ... + x[n - N])

and

y[n] = (1/N) * (x[n - 0] + x[n - 1] + x[n - 2] + ... + x[n-N+1])

Then

y[n] - y[n - 1] = (1/N) * (x[n - 0] - x[n - N])

and we solve for y[n]:

y[n] = y[n - 1] + (1/N) * (x[n - 0] - x[n - N])

Thus in order to do a length N filter recursively we need to store N+1
input values. Slightly more memory storage but a heckuva lot less
computation.

So if your desired filter tap length is 600 you will need an array of
length 601 to do the filter recursively. The filter you have
implemented originally in version 2 will probably be very close in
performance to the one in version 1 since version 2 implements a 599-tap
filter instead of a 600-tap filter, but the results won't be identical.

I have no idea why one method uses a square root and the other does not.
That seems to be an error.
--
% Randy Yates % "Maybe one day I'll feel her cold embrace,
%% Fuquay-Varina, NC % and kiss her interface,
%%% 919-577-9882 % til then, I'll leave her alone."
%%%% <yates@ieee.org> % 'Yours Truly, 2095', *Time*, ELO
http://home.earthlink.net/~yatescr
Back to top
Jack L.
Guest
Posted: Sun Dec 19, 2004 7:57 am    Post subject: Re: Weird behaviour in the nLMS code. Reply with quote
Randy Yates wrote:
Quote:
"Jack L." <jack_nospam@nospam.dk> writes:

Hi Randy,

thanks a bunch for your time on replying. And it's surely gotten my brain
cells working! :) Been re-reading the post for about half an hour now. >_<

I get the math there but I'm a bit confused about the text itself. Are you
showing how to evaluate the Euclidean norm or the FIR-filter faster? I
/think/ you're explaining why the FIR-filter should be 1 tap longer than
TAPS. However, at the same time you mention about the square root operation
which is found in the Euclidean norm, and you refer between 'version 1' and
'version 2' which I used to differentiate the two versions of the Euclidean
norm (they're not used at the same time though).

So, I'm not so sure whether you're explaining the reason for TAPS+1 long
FIR-filter so as to accomodate the faster Euclidean norm evaluation, or are
you explaining about a more efficient FIR-filter evaluation? If it's the
latter, the recursive method seems interesting, it seems to suggest that I
only need to multiply two sample inputs with two of their corresponding
impulse response samples. I've been looking through my school books, can't
seem to find anything about 'recursive FIR-filter'.

I know my question is confusing, it's 5 AM here, don't know why I'm still
awake. It's just the post a bit confusing, as in, few things look mixed up
together to me. :D

--
Mvh / Best regards,
Jack, Copenhagen

The email address is for real. :)
Back to top
Randy Yates
Guest
Posted: Mon Dec 20, 2004 3:30 am    Post subject: Re: Weird behaviour in the nLMS code. Reply with quote
"Jack L." <jack_nospam@nospam.dk> writes:

Quote:
Randy Yates wrote:
"Jack L." <jack_nospam@nospam.dk> writes:

Hi Randy,
thanks a bunch for your time on replying. And it's surely gotten my brain
cells working! :) Been re-reading the post for about half an hour now. >_

I get the math there but I'm a bit confused about the text itself. Are you
showing how to evaluate the Euclidean norm or the FIR-filter faster? I
/think/ you're explaining why the FIR-filter should be 1 tap longer than
TAPS.

Jack,

I think I see why you're confused. The "FIR filter" I was referring to
is in relation to the Euclidean norm and had nothing to do with the
FIR filter further down in your code.

The technique I described is indeed a technique that can be applied to
any FIR filter of the form h[n] = a, nBegin <= n <= nEnd. Part of the
Euclidean norm computation can be seen as an FIR filter of this form
where the input data to the filter is ref_0^2[n] and the factor "a"
above is 1.

In any case, as I had stated in my previous post, TAPS should be
defined to be the number of taps of the filter, which is also the
number of "taps" in the Euclidean norm computation.

Quote:
However, at the same time you mention about the square root operation
which is found in the Euclidean norm, and you refer between 'version 1' and
'version 2' which I used to differentiate the two versions of the Euclidean
norm (they're not used at the same time though).

To do a Euclidean norm, you must first compute the sum of the squares, then
take the square root. Both Version 1 and Version 2 of the code in your original
post compute the sum of the squares. However, Version 2 fails to subsequently
take the square root. As I stated, this failure appears to be an error.

Quote:
So, I'm not so sure whether you're explaining the reason for TAPS+1 long
FIR-filter so as to accomodate the faster Euclidean norm evaluation, or are
you explaining about a more efficient FIR-filter evaluation?

In no case is any filter, whether it be the adaptive filter or the Euclidean
norm "sum of the squares filter," TAPS+1 taps long. As I've stated and restated,
TAPS is the number of taps in the filter.

Quote:
If it's the
latter, the recursive method seems interesting, it seems to suggest that I
only need to multiply two sample inputs with two of their corresponding
impulse response samples. I've been looking through my school books, can't
seem to find anything about 'recursive FIR-filter'.

It only works when an FIR's impulse response is a constant value, thus it's more
like a "trick" than a theoretical result so you may not see it in a textbook.

Quote:
I know my question is confusing, it's 5 AM here, don't know why I'm still
awake. It's just the post a bit confusing, as in, few things look mixed up
together to me. :D

The only thing that I neglected to say is that the sum of the squares in the
Euclidean norm is a type of FIR, and it is that FIR I was referring to.

I'm sorry I wasn't more clear.
--
% Randy Yates % "Ticket to the moon, flight leaves here today
%% Fuquay-Varina, NC % from Satellite 2"
%%% 919-577-9882 % 'Ticket To The Moon'
%%%% <yates@ieee.org> % *Time*, Electric Light Orchestra
http://home.earthlink.net/~yatescr
Back to top
Randy Yates
Guest
Posted: Mon Dec 20, 2004 3:35 am    Post subject: Re: Weird behaviour in the nLMS code. Reply with quote
Randy Yates <yates@ieee.org> writes:
Quote:
[...]
"Jack L." <jack_nospam@nospam.dk> writes:
However, at the same time you mention about the square root operation
which is found in the Euclidean norm, and you refer between 'version 1' and
'version 2' which I used to differentiate the two versions of the Euclidean
norm (they're not used at the same time though).

To do a Euclidean norm, you must first compute the sum of the squares, then
take the square root. Both Version 1 and Version 2 of the code in your original
post compute the sum of the squares. However, Version 2 fails to subsequently
take the square root. As I stated, this failure appears to be an error.

Swap "Version 1" and "Version 2" - I mixed them up here. Sorry.
--
% Randy Yates % "Midnight, on the water...
%% Fuquay-Varina, NC % I saw... the ocean's daughter."
%%% 919-577-9882 % 'Can't Get It Out Of My Head'
%%%% <yates@ieee.org> % *El Dorado*, Electric Light Orchestra
http://home.earthlink.net/~yatescr
Back to top
Jack L.
Guest
Posted: Mon Dec 20, 2004 5:56 am    Post subject: Re: Weird behaviour in the nLMS code. Reply with quote
Peter K. wrote:
Quote:
Jack L. wrote:

That's part of the point of the NLMS algorithm: it converges over a
much better range of step-sizes.

What value are you using for stepsize_0? Make it smaller and the LMS
should converge quite happily... provided the code does implement the
LMS :-)

Hi Peter, so far, I don't seem to be able to exceed 0.2 for nLMS and about

0.08 for LMS (otherwise, the system diverges/goes unstable). As the former
group managed to get up to 2.0 I'm wondering if my code is incorrect. o_O

--
Mvh / Best regards,
Jack, Copenhagen

The email address is for real. :)
Back to top
Jack L.
Guest
Posted: Mon Dec 20, 2004 6:15 am    Post subject: Re: Weird behaviour in the nLMS code. Reply with quote
Randy Yates wrote:
Quote:
"Jack L." <jack_nospam@nospam.dk> writes:


Randy, thanks a lot, that cleared all confusions! :-)

--
Mvh / Best regards,
Jack, Copenhagen

The email address is for real. :)
Back to top
Peter J. Kootsookos
Guest
Posted: Mon Dec 20, 2004 7:58 am    Post subject: Re: Weird behaviour in the nLMS code. Reply with quote
"Jack L." wrote

Quote:
Hi Peter, so far, I don't seem to be able to exceed 0.2 for nLMS and about
0.08 for LMS (otherwise, the system diverges/goes unstable). As the former
group managed to get up to 2.0 I'm wondering if my code is incorrect. o_O

I've seen cases where the LMS blows up for mu greater than 0.02, so that's
not an issue.

What was the value of "b" you used?

Ciao,

Peter K.
Back to top
Shawn Steenhagen
Guest
Posted: Mon Dec 20, 2004 7:58 am    Post subject: Re: Weird behaviour in the nLMS code. Reply with quote
"Jack L." <jack_nospam@nospam.dk> wrote in message
news:%tWwd.76207$Vf.3615633@news000.worldonline.dk...
Quote:
[SNIP]
Code:

Global initialisations:
#define TAPS // filter length
float ref_0[TAPS] // delay line
float w_0[TAPS] // weight vectors
float euclidean=eps;

Process_Sample()
Quote:
{
int n;
double y, error, norm_my, d;

/* Version 1 of the Euclidean norm BEGIN */
euclidean -= ref_0[TAPS-1]*ref_0[TAPS-1];
euclidean += ref_0[0]*ref_0[0];
/* Version 1 of the Euclidean norm END */

/* Version 2 of the Euclidean norm BEGIN */
euclidean = 0;
for (n = 0; n < TAPS; n++) {
euclidean += ref_0[n] * ref_0[n];
}
euclidean = sqrt(euclidean);
/* Version 2 of the Euclidean norm END */

// Delay-line
for (n=TAPS-1; n>=1; n--) {
ref_0[n] = ref_0[n-1];
}

// Inserting new sample into the delay line
ref_0[0] = u0_[0];

// FIR-filter
y = w_0[0]*ref_0[0];
for (n=1; n<TAPS; n++) {
y += w_0[n]*ref_0[n];
}

error = d - y;

norm_my = stepsize_0 * error / (b + euclidean);
for (n=0; n<TAPS; n++) {
w_0[n] += ref_0[n] * norm_my;
}
}


Jack,




A couple of comments:



By definition, Normalize LMS divides the update term by the total power in
the regressor i.e. sum of w_0[i]*w_0[i]. not the the sqrt of the sum of
squares. Version 1 is using a trick to do this sum of squares in one step
by subtracting the oldest sq term first, and then adding the contribution of
the newest. So you must initialize it to zero, (or initialize it to
epsilon). I have seen problems with version 1 on floating point machines
when the input gets small or goes away, the subtract-then-add-approach runs
into precision problems and goes negative! (which will cause herrendeous
problems). Also you will note that the first time through your loop the
power is zero, a situation which should be avoided. This means you should
calculate the regressor power after doing the delay line update. Version 2
is doing the sum of squares directly, but then taking the sqrt, so the
normalization term is not right (may often be too small which means divding
by too small of a number, thus boosting your update term too much, thus
requiring smaller step sizes.)



You have to be careful of the order of doing the estimate, the LMS update
and the buffer update: Here is an excerpt from a matlab m-file that I have
used countless times over the years:



for icnt=1:nump

x(1)=xin(icnt); % add newest input to buffer.

yh=x'*b; % calculate estimate

e=ydes(icnt)-yh; % calculate error

b = b + mu*e*x/(eps+x'*x); % update forward weights

x(2:nb)=x(1:nb-1); % update buffer for estimate

end

Hope this helps

-Shawn

www.appliedsignalprocessing.com
Back to top
Jack L.
Guest
Posted: Tue Dec 21, 2004 2:32 pm    Post subject: Re: Weird behaviour in the nLMS code. Reply with quote
Peter J. Kootsookos wrote:
Quote:
"Jack L." wrote


What was the value of "b" you used?

Ciao,

Peter K.

Thanks for your response, Peter.

Interesting, mu = 0.02 for the LMS? Could that have been the result of a
very long filter? (given that the 0 < mu < 2 / M*Smax formula was used,
where M = filter length, Smax = power spectral density of the input vector)

I had set "b" to 0.001 or 0.0001. I've realised this number can be a bit too
small because a small 'euclidean' (the variable) + a small 'b' can give a
huge norm_my. I may have put a bit too much emphasis on Haykin's words:
"delta is a small value". I'm going to try 0.5 and 1.

The DSP board's voltage level is between +/- 1.4V, and generally, the signal
exhibiting from the human ears we are interested in, is about 100 mV.

--
Mvh / Best regards,
Jack, Copenhagen

The email address is for real :)
Back to top
Jack L.
Guest
Posted: Sun Dec 26, 2004 5:02 am    Post subject: Re: Weird behaviour in the nLMS code. Reply with quote
Shawn Steenhagen wrote:
Quote:
"Jack L." <jack_nospam@nospam.dk> wrote in message
news:%tWwd.76207$Vf.3615633@news000.worldonline.dk...

By definition, Normalize LMS divides the update term by the total
power in the regressor i.e. sum of w_0[i]*w_0[i]. not the the sqrt of
the sum of squares.

Shawn, thanks for your response. It's sort of confirming what I've seen,
that is, the algorithm seems to perform much better by not taking the sqrt
of the sum of sqares. But I get confused because in Simon Haykin's book
"Adaptive Filter Theory", 4th edition, he says that mu is normalised with
respect to the Euclidean norm. If I'm not wrong, the Euclidean norm is given
by: sqrt(sum(w_0[i]*w_0[i])) <-- half MatLab, half C-notation. I don't know
what to believe, the book or what I observe (+ what you say)?

Quote:
Also you will note that the first time through your loop the power is
zero, a situation which should be avoided. This means you should
calculate the regressor power after doing the delay line update.


Yes, you're right. Thanks again for the correction and for the m-file - it
shows how simple the nLMS algorithm can be. :)

--
Mvh / Best regards,
Jack, Copenhagen

The email address is for real :)
Back to top
 
Post new topic   Reply to topic    CASTalk.com Forum Index -> DSP All times are GMT
Goto page 1, 2  Next
Page 1 of 2

 
 You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum




VoIP Electronics Powered by phpBB