# Weird behaviour in the nLMS code.

Started by December 18, 2004
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. :)