On 8/12/11 11:33 AM, maury wrote:
> On Aug 11, 7:35 am, JohnPower<asdfghjkl...@googlemail.com> wrote:
>> On Aug 11, 3:16 am, robert bristow-johnson<r...@audioimagination.com>
>> wrote:
>>
>>
>>
>>
>>
>>> On 8/10/11 1:56 PM, JohnPower wrote:
>>
>>>> On 3 Aug., 15:36, robert bristow-johnson<r...@audioimagination.com>
>>>> wrote:
>>>>> On 8/3/11 8:04 AM, JohnPower wrote:
>>
>>>>>> For my master's thesis I'm currently trying to implement an active
>>>>>> noise cancellation (ANC) system on a ADSP-21369 EZ-KIT Lite. I've
>>>>>> chosen the Normalized Leaky LMS algorithm to update (451) coefficients
>>>>>> of a FIR filter.
>>
>>>>> John, would you mind plopping down the relevant equations for the NLLMS?
>>
>>>>> i know how NLMS works and i have an idea how to make the coefficient
>>>>> updating to be leaky, but i dunno how to get a handle on your problem
>>>>> without seeing the equations.
>>
>>>>> --
>>
>>>>> r b-j r...@audioimagination.com
>>
>>>>> "Imagination is more important than knowledge."
>>
>>>> this being my first thread in a usenet group, I've somehow expected to
>>>> get a notification on new messages... well, sorry for the delay (and
>>>> also the double post...).
>>
>>>> my equation for the NLLMS is (vectors uppercase, scalars lowercase):
>>
>>>> W[n+1] = W[n] * l + a*e[n]*X[n] / (s + power(X[n])
>>
>>> i guessed pretty good:
>>
>>> g[n] = (2*mu) / mean{ x[n]^2 }
>>
>>> h[k]<-- p*h[k] - g[n] * e[n]*x[n-k] at time n.
>>
>>> different symbols and i missed the "s" term.
>>
>>>> where
>>>> W - coefficient vector
>>>> X - vector which contains last N samples from my template which was
>>>> treated with the impulse response P I've measured via MLS
>>
>>> since you're not using the shorthand i was using, then let's explicitly
>>> expand the notation so we know exactly what is being done. using C-like
>>> notation for 2-dim arrays:
>>
>>> coef vector: W[n] = { w[0][n], w[1][n], ... w[k][n], ... w[K-1][n] }
>>
>>> FIR input vector: X[n] = { x[n], x[n-1], ... x[n-k], ... x[n-K+1] }
>>
>>> we need to do this, because i see an ambiguity in your equation below.
>>> can you confirm this?
>>
>> my FIR input vector Xo[n] has the same dimensions as the coefficients
>> vector W[n]:
>>
>> Xo[n] = { x[0], x[1], ... , x[N] }
this can't be right. i think you mean { x[n], x[n-1], ... x[n-N+1] }.
and why isn't this vector named "X"? if this is not X, then what is in
the X vector?
>> W[n] = { w[0], w[1], ... , w[N] }
this should have two indices on each vector element.
W[n] = { w[0][n], w[1][n], ... , w[N-1][n] }
>>
>> where N is the filter order. When writing W[n+1] I mean the
>> coefficient vector which is used
>> to compute the next output. I've used the nomenclature from Diniz'
>> book on 'Adaptive Filtering'
>>
>>
>>
>>
>>
>>>> l - leakage 0< l< 1, in my case it is 0.999
>>
>>> which is like a first-order filter pole. just like my "p".
>>
>>>> a - factor to influence speed of convergence, I get rather good
>>>> results with setting this to -0.017
>>
>>> this, or something proportional to it, is called the "adaptation gain".
>>
>>>> e - current error value from error microphone
>>
>>> so is the differencing done in the air? this does not make sense for an
>>> LMS. e[n] is, as far as i can tell, *computed* from the output of the
>>> FIR (what is defined by W[n] and what i thought is y[n]) and from
>>> subtracting the "desired" signal (what i would think comes from this
>>> microphone). you need to clear this up.
>>
>> yes, it's done in the air.
>> The microphone serves two purposes:
>> 1) record the 'desired' signal
>> 2) measure the error betwen output signal and current noise
>>
>> It's something like a hybrid between a feedforward and feedback
>> system, I guess this is
>> what makes it so hard to understand - and explain :-/
>>
>>
>>
>>>> power(X[n]) - the power of X[n] estimated by doing a scalar
>>>> multiplication of X[n] with itself,
>>
>>> i don't get this. "scalar multiplication" of a vector with a vector?
>>> do you mean dot product?:
>>
>> sorry, my bad. Yes I mean dot product (in German, Skalarprodukt is the
>> common term
>> which I mixed up with scalar multiplication...)
>>
>>
>>
>>> N-1
>>> power(X[n]) = SUM{ x[n-k]^2 }
>>> k=0
>>
>>> this is what i called "mean{ x[n]^2 }".
>>
>>> > this is the "normalize" part
>>
>>> dividing by that is the "normalize" part, i think you mean.
>>
>> yes
>>
>>
>>
>>>> s - small value to prevent coefficients getting too big when
>>>> power(X[n]) is small
>>
>>> well, we don't ever want to divide by 0 but when the power is small,
>>> e[n] and X[n] should also be small. so when the power is small, the
>>> adaptation gets slower (which may be just fine, what you want). let's
>>> assume that power(X[n]) is large enough that s is negligible, okay?
>>
>>> how are you computing the division? explicitly? sometimes, for NLMS,
>>> when division is expensive. it can be computed (tracked) with a
>>> multiply, compare, and negative feedback.
>>
>> my actual C code looks like this:
>>
>> tempUpdateFactor= a * e / templPower;
>>
>> As one can see, I've left out the s as it doesn't seem to be the
>> problem here...
>> So yes, I do explicit division. I suppose the "mulitply, compare,
>> negative feedback"
>> you've suggested would be to have the algorithm run faster? I haven't
>> looked into
>> this exact line for optimization yet because the vector calculations
>> are much
>> more time expensive.
>>
>>
>>
>>
>>
>>
>>
>>>> after calculating W[n+1], I've got an FIR Filter doing
>>
>>>> K-1
>>>> y[n] = SUM{ W[k] * Xo[n-k] }
>>>> k=0
>>
here it's K-1
>>>> where
>>>> Xo - template which did not go through the impulse response P
>>
>>> right here, i see a problem regarding the difference between the vector
>>> X[n] and, apparently, the scalar values Xo[n-k]. is the vector X[n]
>>> made up of the samples Xo[n-k]? if yes, your notation is bad (the
>>> elements of X[n] should be small-case, just like y[n] is). if no, this
>>> NLMS algorithm is messed up.
>>
>> I don't get it, if it was messed up this bad, it wouldn't work at all,
>> would it?
>> But my darned system works excellent for some time and then fades out
>> gracefully...
>>
>> corrected notation:
>>
>> N-1
>> y[n] = SUM{ w[k] * xo[N-k] }
>> k=0
>>
first it was K, then it's N, now it N-1.
need to get the symbols down, John. i don't think you have it clear
about Xo and X
>> I actually do this by using a library function to
>> compute the dot product of the two vectors W[n] an Xo[n].
>>
>> Xo and X are two different vectors. I've uploaded a new signal graph.
>> It was
>> revised to be coherent with the symbols in my equations:http://www.abload.de/img/signal_flow_coherentdusl.png
>>
>>
>>
>>> is P the impulse response of the room or whatever other path (comparable
>>> to W[n])?
>>
>> Yes, P is the impulse response of the room.
>> To be more specific, it is the plant model of the system
>> DAC->amplifier->transducer->microphone->amplifier->ADC
so it's a system identification thing (which adaptive filters are
supposed to do).
>
> O.K., here is what may be happening. I took a quick look at the cited
> paper. Read the paper again. At the output of the adaptive filter you
> have an amplifier. You also state that P(z) is the room impulse
> response. I think you have made a mistake here. In the paper M(z) is
> the impulse response of the room, the speaker, etc. There is a H(z)
> where you have the amplifier at the output of the adaptive filter.
> H(z) represents the impulse response between the adaptive filter
> output to the speaker to the microphone (page 286 - *H(z) the transfer
> function between the control loudspeaker and the error microphone*).
> Then, P(z) is the same as H(z), not the room impulse response. Take
> another look at the paper.
>
and i'm not convinced that the schematic diagram you refer us to is
right (that is, if you're trying to do noise cancellation). i don't get
any of this "trigger" stuff. doesn't belong there.
i believe what you're trying to accomplish is pretty close to the
standard LMS or NLMS (or LNLMS), but instead of trying to match the
"plant", which is the path from the ambient noise microphone to the
earpiece, and computing the "error" or difference signal e[n] in the
DSP, you're computing the e[n] in the air. in both cases, the adaptive
filter is trying to minimize e[n]. you can start with a clean LMS
filter (or NLMS or LNLMS) and change the topology a little (it really
isn't changing it, just drawing the border around e[n] differently).
i'm gonna force changing the notation to something more conventional in
the lit.
the regular LMS problem is this:
.-----------. y[n]
.---->| h[k][n] |------------.
| '-----------' |
x[n] --->---| (+)----> e[n] = y[n] - d[n]
| .--------. d[n] |
'---->| p[k] |-------->[-1]--'
'--------'
vector H[n] = { h[0][n], h[1][n], ... h[K-1][n] } is the FIR
vector P = { p[0], p[1], ... p[K-1] } is the plant.
the whole idea of LMS adaptive filtering is to get your FIR impulse
response (vector H[n]) to match the plant (vector P) and when the match
is good, e[n] should start to get small in magnitude. sometimes "d[n]"
is called the "desired" signal and you're trying to get the FIR output,
y[n], to match the plant output d[n]. the error signal e[n] feeds back
somehow to change the h[k] so that h[k] matches p[k] by some measure.
your FIR computes y[n] as
K-1
y[n] = SUM{ h[k][n] * x[n-k] }
k=0
and your "plant" hypothetically creates d[n] as
K-1
d[n] = SUM{ p[k] * x[n-k] }
k=0
but no one is doing that second summation. the physics of the plant is
doing it. now your LNLMS (leaking normalized least mean square) filter
is updating the FIR coefficients
h[k][n+1] = h[k][n]*(1-leak) - mu*e[n]*x[n-k]/mean{ x[n]^2 }
where mu is the adaptation gain (small and positive) and
0 < leak << 1
K-1
mean{ x[n]^2 } = 1/K * SUM{ (x[n-k])^2 }
k=0
the 1/K doesn't have to be computed if it gets incorporated into the mu
coefficient. and this is a simple moving summation, after squaring, add
one term in and subtract the term that is falling off the edge. there
are other ways (like a first-order IIR) to calculate mean{x[n]^2}. you
can add a small positive "s" term to this before dividing, if you want.
you can show (using partial derivatives) that this term with the mu in
it, "nudges" the coefficients in such a way as to drive down the
opposite direction of the gradient of e[n]^2 as a function of h[k][n].
conventionally, you have a microphone for d[n] and you do the
subtraction to of d[n] out of y[n] to get e[n].
now, John, first you need to get this mechanism understood. once it's
understood, we can change it a little. first, since the subtraction
will be done in the air, we recognize that it's really an addition that
is physically done and fiddle with the signs a little:
.-----------. y[n]
.---->| h[k][n] |------------.
| '-----------' |
x[n] --->---| (+)----> e[n] = y[n] + d[n]
| .--------. d[n] |
'---->| p[k] |-------->------'
'--------'
now both y[n] and e[n] are the negatives of what they used to be. just
a sign change, but we want e[n]^2 to be small just the same. the result
is that h[k][n] should converge to the *negative* of p[k]. the FIR is
the negative of the plant, so that when they add, it gets closer to zero.
you do the same song-and-dance to get the gradient of e[n]^2 as a
function of all of the h[k] and it comes out the same:
(d/dh[k])( e[n] )^2 = 2*e[n] * x[n-k]
so, to march down the gradient (with little steps), it's the same:
h[k][n+1] = h[k][n]*(1-leak) - mu*e[n]*x[n-k]/mean{ x[n]^2 }
and there is nothing wrong with drawing the boundary differently, let
the addition be done in the air, and assume the microphone (at the
earpiece) is measuring e[n] instead of d[n]. makes no difference.
i'm wondering if you get the polarity on the mic wrong, what that will
do is drive the h[k] coefficients into the opposite direction. but,
independent of the mic polarity, the h[k] impulse response should come
out to be the negative of the hypothetical p[k]. i wonder if that's a
problem. Tim or Eric or someone, can you comment? can the wrong
polarity in the e[n] mic make this diverge instead of converge? it
seems to me that it should still be adjusting h[k] to drive down e[n]^2,
whether e[n] became negated or not, but i might be wrong about that. if
so, that might answer the OP's problem.
hmmm, as i ponder this, i think that the mic polarity, relative to the
headphone speaker, must be consistent or the sign in front of mu has to
change in the filter adaptation equation.
--
r b-j rbj@audioimagination.com
"Imagination is more important than knowledge."