Weird behaviour in the nLMS code.

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;
 }
}

Jack L. wrote:

> 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.

> 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.

"Jack L." <jack_nospam@nospam.dk> writes:
> [...]
>  /* 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

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

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

> 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

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?

--
Mvh / Best regards,
Jack, Copenhagen

The email address is for real. :)

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

> 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

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. 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. :)

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

> 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. 

> 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. 

> 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. 

> 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.

> 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

Randy Yates <yates@ieee.org> writes:
> [...]
> "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

Peter K. wrote:
> 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. :)