in article 39etu6F60u5gkU1@individual.net, Jon Harris at goldentully@hotmail.com wrote on 03/11/2005 20:12:> "robert bristow-johnson" <rbj@audioimagination.com> wrote in message > news:BE57780C.52D9%rbj@audioimagination.com......>> (i've also done >> 1st order IIR filtering to smooth out audio parameters, such as pitch, but >> hadn't liked it because the tail *never* goes away.) > > A side note, sometimes I've followed the 1st order IIR smoothing filter with a > little routine that checks the result of the filter and if it is "close > enough" and if it is, makes it exactly so. This "post-check" routine is run > infrequently so it doesn't take much processing. It's not elegant, but it is > effective (and quite important in the case of say smoothing the coefficients > of a very narrow low-frequency notch filter).that's one way to deal with it. dunno why (perhaps it's just a matter of taste), but i tend to stay away from inelegant solutions. the moving average filter is nice because it always gets to its final value in N samples. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Question about computing a "moving variance"
Started by ●February 27, 2005
Reply by ●March 11, 20052005-03-11
Reply by ●March 12, 20052005-03-12
in article 1110583041.947845.231670@l41g2000cwc.googlegroups.com, JS at jshima@timing.com wrote on 03/11/2005 18:17:> Hi Robert, > > Yeah, I tried to imply that the filter is not an LTI,if i had bothered to read it through, before typing half a book (which i didn't wanna throw away after doing it), i would have seen it. the thing that threw me off is you used a different "n" and "N", but they're really the same in your case.> For sample mean using sample input x[n]: > > m[n+1] = m[n] + 1/(N+1)*(x[n+1] - m[n]) > = N/(N+1)*m[n] + 1/(N+1)*x[n+1] = a nice first-order IIR filterreally, what you're saying is> m[N+1] = m[N] + 1/(N+1)*(x[N+1] - m[N]) > = N/(N+1)*m[N] + 1/(N+1)*x[N+1] not= a nice first-order IIR filter> it is time > variant since you dont have linear constant coefficients. Basically I > was trying to give Rick another idea of a recursive mean and variance.right. there's a lot of LPFs that give you some meaning of the word "mean".> This is not a "moving variance" as defined,yup.> but it can mock one if you > want to reset the filter at a given number of points "M". In this case > you have a variance that will be equal to the block variance at each > point. > > Definition wise for the mean computation, start with: > > xbar[n] = 1/n* sum(i=1 to n) {x[i]} // standard block mean for n > points > > let: > > xbar[n+1] = 1/(n+1)* sum(i=1 to n+1){x[i]} > > so that: > > xbar[n+1] = 1/(n+1)* sum(i=1 to n){x[i]} + 1/(n+1)*x[n+1] > xbar[n+1] = n/(n+1)*xbar[n] + 1/(n+1) * x[n+1] > > And you get the equation I listed. It is a time-variant filter that > "looks" like an alpha filter with a = n/(n+1) and (1-a) = 1/(n+1), but > certainly cannot be analyzed as such.nope> At any rate, I didnt mean to confuse the issue at hand... :) > > I will look at your coeffs above for the LTI alpha filter, it looks > interesting. Thanks...all's it is, is two things. 1. it's a replacement of the CIC filter that is a moving weighted average of the past samples (where the weights are 1/N for the most current N samples and 0 for all else) with a 1 pole recursive LPF (where the weights are (1-p)*p^(-n-1)). both of these averaging filters have DC gain of 1 (as they should) no matter what N and p are. so how do you relate the two so to make a comparable substitution? i believe the way to make the substitution the most similar is, keeping the DC gain the same (which is automatically done) is to equate the effective delay of the two filters by setting p = (N-1)/(N+1) to substitute for an N tap CIC filter. that way, in both cases, there is a delay (at least at DC) of (N-1)/2 samples and, in both cases, there is some way we can say "we're averaging, at least in some sense, the most current N samples and the squares of the most current N samples". 2. then, if in both cases you are are averaging (in some sense of the word) the most current N samples, in the CIC case, we know we have to multiply the "biased" variance by N/(N-1) to get unbiased variance, and i would argue that you to the same for the 1 pole IIR averager as long as p = (N-1)/(N+1). otherwise, i would not have any idea on how to unbias the variance coming out of the IIR averager. actually i *do* have an idea on how to do it!! in both cases, we define the expectations: mean: E{ x[n] } = integral{ x * p(x) dx} = M variance: E{ (x[n]-m)^2 } = integral{ (x-M)^2 * p(x) dx} = sigma^2 and we know that mean square: E{ (x[n])^2 } = integral{ x^2 * p(x) dx} = sigma^2 + M^2 where the integral is over all x and p(x) is the p.d.f. of x[n]. now with the CIC estimations we have: N-1 est. mean: m[n] = (1/N) * SUM{ x[n-i] } i=0 N-1 est. variance: v[n] = 1/(N-1) * SUM{ (x[n-i] - m[n])^2 } i=0 N-1 = N/(N-1))*[1/N * SUM{(x[n-i])^2} - ( m[n] )^2] i=0 = N/(N-1))*[ s[n] - ( m[n] )^2] where N-1 s[n] = (1/N) * SUM{ (x[n-i])^2 } i=0 here 1/N for 0 <= i < N is just the impulse response, h[i], of the CIC averager. but with the 1 pole LPF we have: +inf est. mean: m[n] = SUM{ (1-p)/p*p^(-i) * x[n-i] } i=0 +inf = (1-p)/p * SUM{ p^(-i) * x[n-i] } i=0 (looks a lot like a Z transform) est. variance: +inf v[n] = (1+p)/(2*p) * SUM{ (1-p)/p * p^(-i)*(x[n-i] - m[n])^2 } i=0 +inf = (1+p)/(2*p) * (1-p)/p * SUM{ p^(-i)*(x[n-i] - m[n])^2 } i=0 the leading (1+p)/(2*p) factor is what i am saying we need to unbias the estimation of the variance so that E[ v[n] ] = sigma^2. (1-p)/p*p^(-i) is just the impulse response, h[i], of the IIR averager. now i am not sure in this case, but does +inf SUM{ (1-p)/p*p^(-i)*(x[n-i] - m[n])^2 } = s[n] - ( m[n] )^2 ? i=0 where +inf s[n] = SUM{ (1-p)/p*p^(-i) * (x[n-i])^2 } i=0 i think it does. so i predict that the unbiased estimated variance is: v[n] = (1+p)/(2*p) * [ s[n] - ( m[n] )^2 ] just like with the CIC filters but instead of multiplying by N/(N-1) where multiplying by (1+p)/(2*p) because we are equating p/(1-p) and (N-1)/2 . this comes from functioning 1 pole IIR filters: estimated mean: m[n] = p*m[n-1] + (1-p)*x[n] estimated mean square: s[n] = p*s[n-1] + (1-p)*(x[n])^2 where p is the pole value. we've already worked out that the group delay (and phase delay) at DC is tau(0) = p/(1-p) and we want to set that to tau(w) = (N-1)/2 for the N delay CIC filter. we know that for the N delay CIC averager, we unbias the estimate for variance by multiplying by N/(N-1) so if tau(0) = p/(1-p) = (N-1)/2 and we solve for N, it's N = (1+p)/(1-p) and N/(N-1) = (1+p)/(2*p) . oh geez, this is starting to look laborious. hey can you, Jim, or Clay or Randy or Rick or anybody work this out and figure out what the expectations of m[n] and v[n] are? as before we know that when n <> i, then E{ x[n]*x[i] } = E{ x[n] } * E{ x[i] } = M^2 because x[n] and x[i] are controlled by the same statistics. and when n = i, then E{ x[n]*x[i] } = E{ x[n]^2 } = sigma^2 + M^2 who wants to blast through that? (i am weary, but still curious.) doing so earns you, well nothing tangible, but maybe some comp.dsp status. :-) who wants to slug that out? i wonder if, and will sorta bet that, the result comes out to be E{ s[n] - ( m[n] )^2 } = 2*p/(1+p) * sigma^2 that's what i predict based on my "equate the two filter delays" strategy. that is tau(0) = p/(1-p) = (N-1)/2 (this is an N tap CIC, not an N+1 tap CIC as was your case, Jim.) if this is true then in one case you multiply by N/(N-1) to unbias the estimated variance and in the other you multiply by (1+p)/(2*p) which makes, in both cases E{ v[n] } = sigma^2 which is right on the money. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●March 13, 20052005-03-13
in article BE57E7F5.52F4%rbj@audioimagination.com, robert bristow-johnson at rbj@audioimagination.com wrote on 03/12/2005 00:27:> oh geez, this is starting to look laborious. hey can you, Jim, or Clay or > Randy or Rick or anybody work this out and figure out what the expectations > of m[n] and v[n] are?i worked it out and my prediction was confirmed. if you use a 1 pole IIR filter (with DC gain = 1) to do averaging for the moving estimated mean and moving estimated mean square of some data, the scaler you need to unbias the estimate of variance derived from that is (1+p)/(2*p) where p is the pole value. this corresponds to a CIC averager of N samples and the scaler to unbias the estimate of variance is N/(N-1). the neat coincidence (if it really is a conincidnce) is that if p is set so that the unbiasing scalers are the same, that is (1+p)/(2*p) = N/(N-1) then the delay of the IIR averager, p/(1-p), is the same as the delay of the N sample CIC, (N-1)/2. that implies, in both cases, that they are averaging, in some sense of the word, the same number of samples which is N in the CIC case. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●March 14, 20052005-03-14
"robert bristow-johnson" <rbj@audioimagination.com> wrote in message news:BE57B4E4.52EC%rbj@audioimagination.com...> in article 39etu6F60u5gkU1@individual.net, Jon Harris at > goldentully@hotmail.com wrote on 03/11/2005 20:12: > > > "robert bristow-johnson" <rbj@audioimagination.com> wrote in message > > news:BE57780C.52D9%rbj@audioimagination.com... > ... > >> (i've also done > >> 1st order IIR filtering to smooth out audio parameters, such as pitch, but > >> hadn't liked it because the tail *never* goes away.) > > > > A side note, sometimes I've followed the 1st order IIR smoothing filter witha> > little routine that checks the result of the filter and if it is "close > > enough" and if it is, makes it exactly so. This "post-check" routine is run > > infrequently so it doesn't take much processing. It's not elegant, but itis> > effective (and quite important in the case of say smoothing the coefficients > > of a very narrow low-frequency notch filter). > > that's one way to deal with it. dunno why (perhaps it's just a matter of > taste), but i tend to stay away from inelegant solutions. the moving > average filter is nice because it always gets to its final value in N > samples.Agreed. OTOH, often times I don't have the memory/computation bandwidth to do an FIR filter for smoothing for the fairly-long smoothing times I've needed in my applications. So that is where the first order IIR + post-check routine is useful. But given the option, I would prefer the moving average as well.
Reply by ●March 15, 20052005-03-15
robert bristow-johnson wrote:> i worked it out and my prediction was confirmed. if you use a 1 pole IIR > filter (with DC gain = 1) to do averaging for the moving estimated mean and > moving estimated mean square of some data, the scaler you need to unbias the > estimate of variance derived from that is (1+p)/(2*p) where p is the pole > value. this corresponds to a CIC averager of N samples and the scaler to > unbias the estimate of variance is N/(N-1). the neat coincidence (if it > really is a conincidnce) is that if p is set so that the unbiasing scalers > are the same, that is > > (1+p)/(2*p) = N/(N-1) > > then the delay of the IIR averager, p/(1-p), is the same as the delay of the > N sample CIC, (N-1)/2. that implies, in both cases, that they are > averaging, in some sense of the word, the same number of samples which is N > in the CIC case. >I need some clarification. If you set p such that (1+p)/(2*p) = N/(N-1) then p becomes p = (N+1)/(N-1) Assuming N is positive (I think that's safe), doesn't this make p>1? Isn't that a Bad Thing for an IIR filter, or am I totally confused? Seems likely since I'm highly prone to confusion over simple matters. To wit, my head musta been on backwards when I suggested using a Goertzel to get the DC value. This was probably obvious to a lot of you, but it took me a while to come to the realization. Although it would work, you'd hafta recalculate it for every new sample. In which case, you may as well go the traditional route. -- Jim Thomas Principal Applications Engineer Bittware, Inc jthomas@bittware.com http://www.bittware.com (603) 226-0404 x536 Having a smoking section in a restaurant is like having a peeing section in a swimming pool.
Reply by ●March 15, 20052005-03-15
in article 113dq4c29g9k066@corp.supernews.com, Jim Thomas at jthomas@bittware.com wrote on 03/15/2005 08:52:> I need some clarification. If you set p such that > > (1+p)/(2*p) = N/(N-1) > > then p becomes > > p = (N+1)/(N-1) >are you sure? i think it's the other way around? p = (N-1)/(N+1) lesseee... (1+p)/(2*p) = N/(N-1) (1+p) = N/(N-1)*(2*p) 1 = N/(N-1)*(2*p) - p 1 = (2*N/(N-1) - 1)*p p = 1/(2*N/(N-1) - 1) = (N-1)/(2*N - (N-1)) = (N-1)/(N+1) -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●March 15, 20052005-03-15
robert bristow-johnson wrote:> in article 113dq4c29g9k066@corp.supernews.com, Jim Thomas at > jthomas@bittware.com wrote on 03/15/2005 08:52: >>p = (N+1)/(N-1) >> > are you sure? i think it's the other way around? p = (N-1)/(N+1)Of course you're right. But I haven't had much luck getting this to compute a moving average (and thus, a moving variance). Is this correct? y[n] = p*y[n-1] + (1-p)x[n] -- Jim Thomas Principal Applications Engineer Bittware, Inc jthomas@bittware.com http://www.bittware.com (603) 226-0404 x536 Having a smoking section in a restaurant is like having a peeing section in a swimming pool.
Reply by ●March 15, 20052005-03-15
in article 113ec4tmt8h1h89@corp.supernews.com, Jim Thomas at jthomas@bittware.com wrote on 03/15/2005 13:59:> robert bristow-johnson wrote: >> in article 113dq4c29g9k066@corp.supernews.com, Jim Thomas at >> jthomas@bittware.com wrote on 03/15/2005 08:52: >>> p = (N+1)/(N-1) >>> >> are you sure? i think it's the other way around? p = (N-1)/(N+1) > > Of course you're right. But I haven't had much luck getting this to > compute a moving average (and thus, a moving variance). Is this correct? > > y[n] = p*y[n-1] + (1-p)x[n]yes. that is an LPF with DC gain of 1, so it has some legitimacy in being called an "averager". but it isn't the same moving average as the N tap CIC filter. the impulse response is different, so the weighting factors are different. in the N tap CIC, you have weighting factors of 1/N for the most current N samples and 0 for every other sample. for the 1 pole IIR (above) your weighting factors are (1-p)/p * p^(-n+i) for all i<=n and zero for i>n. they are both averagers, they can both be used for a sort of weighted moving average, if p = (N-1)/(N+1), then they both average, in a sense, the most current N samples, but the IIR is actually averaging the past infinity samples, it's just that some of them have negligible weight. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●March 20, 20052005-03-20
On Tue, 15 Mar 2005 13:59:38 -0500, Jim Thomas <jthomas@bittware.com> wrote:>robert bristow-johnson wrote: >> in article 113dq4c29g9k066@corp.supernews.com, Jim Thomas at >> jthomas@bittware.com wrote on 03/15/2005 08:52: >>>p = (N+1)/(N-1) >>> >> are you sure? i think it's the other way around? p = (N-1)/(N+1) > >Of course you're right. But I haven't had much luck getting this to >compute a moving average (and thus, a moving variance). Is this correct? > >y[n] = p*y[n-1] + (1-p)x[n] >Hi Jim, just saw you post. To add 2 cents to Robert B-J's words: your y[n] = p*y[n-1] + (1-p)x[n] difference equation is, as you said, and exponential averager [also called a "leaky integrator"],when p is some real number less that one. I'll bet you already knew that. It's the IIR filter (used for noise reduction) that I described in Section 8.5 of the 1st edition of my book [Section 11.5 of the 2nd edition]. Please note: my alpha is equal to your above (1-p). Ya' know Jim, I was thinking. When I see block diagrams of that exponential averager in the literature, I always see the input sequence multiplied by (1-p) [my alpha] before the feedback loop as shown below. x(n)--->(*)------>(+)-----------+----> y(n) ^ ^ | | | v [1-p] | [z^-1] FIG 1. | | +----(*)-----+ ^ | p It seems to me that multiplying x(n) by a number less than one [(1-p)] reduces the input signal-to-quantization noise ratio if we use a fixed-point number system. If that's true, maybe it's better to use the following FIG 2 block diagram: x(n)------->(+)-----------+------->(*)---> y(n) ^ | ^ | v | | [z^-1] [1-p] | | +----(*)-----+ ^ FIG 2. | p I wonder if anyone has modeled (in fixed-pt math) those two versions of an exponential averager, and knows which version should be used. (I don't have an easy way to model "fixed-point" math.) See Ya, [-Rick-]
Reply by ●March 20, 20052005-03-20
Rick Lyons wrote: ...> x(n)--->(*)------>(+)-----------+----> y(n) > ^ ^ | > | | v > [1-p] | [z^-1] FIG 1. > | | > +----(*)-----+ > ^ > | > p...> x(n)------->(+)-----------+------->(*)---> y(n) > ^ | ^ > | v | > | [z^-1] [1-p] > | | > +----(*)-----+ > ^ FIG 2. > | > pRick, What about overflow? Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
