On 6/16/12 1:43 PM, niarn wrote:>> it *shouldn't* be flat. the roundoff error goes through a simple >> first-order "differentiator" with a zero right on z=1. the noise floor >> should be zero at DC and twice as big at Nyquist (compared to if no >> noise shaping was done). that's what the "noise shaping" is supposed to be. > > True, right at DC and an excruciatingly narrow band around DC the noise > floor goes to -inf. Apart from that it is flat as a pancake. > > I'm not sure I follow you. I believe that the error (neglecting the error > from the multiplication) 'sees' the same transfer function as the input > does...you are correct. i had to draw it out and look at it again. i misremembered it. in *comparison* with the noise spectrum that would result of noise shaping wasn't used, the noise was steered away from DC (where it goes from some level to zero) and is steered toward Nyquist (where it goes from some level to twice the level). but, either way, it is inside the filter loop and, in fact, the noise transfer function is the same as the DC blocking transfer function. if noise shaping wasn't used, you would see the noise spectrum having an LPF shape. and putting it through a differentiator makes it flat nearly everywhere (except at DC). -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
DC Blocking Filter with Noise Shaping
Started by ●June 12, 2012
Reply by ●June 16, 20122012-06-16
Reply by ●June 16, 20122012-06-16
On 6/16/12 1:52 PM, niarn wrote:>>> Is there a reference to where the solution was first 'envisioned'? >> >> by me, at the dspguru site (and it was posted here just before Grant >> used it on the dspguru site). i dunno if anyone beat me to that idea or >> not. i s'pose some Bell Systems engineer thought of this in the mid >> 60s, but i dunno. >> > > The reason I asked is that I once worked at a place where the exact same DC > blocking algo was used. And that code went into production as early as 2001 > and probably 'envisioned' in 2000. But unfortunately I don't know if the > guy who was responsible for the code could have read your posts if they > were in the public domain in 2000.i put it here (on comp.dsp) about 12 years ago (Jun 2000): http://groups.google.com/group/comp.dsp/browse_frm/thread/661e0878e7af7813/ and, as far as i know, it didn't appear before that. (however, there are some products i have worked on, dating back to 1993, that have this trick.) i first started learning about noise shaping in 1991, maybe even earlier, from papers and presentations from Stanley Lipshitz. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●June 16, 20122012-06-16
On 6/16/12 2:16 PM, robert bristow-johnson wrote:> On 6/16/12 1:52 PM, niarn wrote: >>>> Is there a reference to where the solution was first 'envisioned'? >>> >>> by me, at the dspguru site (and it was posted here just before Grant >>> used it on the dspguru site). i dunno if anyone beat me to that idea or >>> not. i s'pose some Bell Systems engineer thought of this in the mid >>> 60s, but i dunno. >>> >> >> The reason I asked is that I once worked at a place where the exact >> same DC >> blocking algo was used. And that code went into production as early as >> 2001 >> and probably 'envisioned' in 2000. But unfortunately I don't know if the >> guy who was responsible for the code could have read your posts if they >> were in the public domain in 2000. > > > i put it here (on comp.dsp) about 12 years ago (Jun 2000): > > > http://groups.google.com/group/comp.dsp/browse_frm/thread/661e0878e7af7813/ > > and, as far as i know, it didn't appear before that.actually, it was a year before that (Apr 1999): http://groups.google.com/group/comp.dsp/browse_frm/thread/54148db1e1067f3a/ -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●June 19, 20122012-06-19
On 6/16/2012 2:29 PM, robert bristow-johnson wrote:> On 6/16/12 2:16 PM, robert bristow-johnson wrote: >> On 6/16/12 1:52 PM, niarn wrote: >>>>> Is there a reference to where the solution was first 'envisioned'? >>>> >>>> by me, at the dspguru site (and it was posted here just before Grant >>>> used it on the dspguru site). i dunno if anyone beat me to that idea or >>>> not. i s'pose some Bell Systems engineer thought of this in the mid >>>> 60s, but i dunno. >>>> >>> >>> The reason I asked is that I once worked at a place where the exact >>> same DC >>> blocking algo was used. And that code went into production as early as >>> 2001 >>> and probably 'envisioned' in 2000. But unfortunately I don't know if the >>> guy who was responsible for the code could have read your posts if they >>> were in the public domain in 2000. >> >> >> i put it here (on comp.dsp) about 12 years ago (Jun 2000): >> >> >> http://groups.google.com/group/comp.dsp/browse_frm/thread/661e0878e7af7813/ >> >> >> and, as far as i know, it didn't appear before that. > > actually, it was a year before that (Apr 1999): > > > http://groups.google.com/group/comp.dsp/browse_frm/thread/54148db1e1067f3a/I used essentially the same technique to minimize the error occasioned by discarding the low bits of integer division in an HC6811 back around 1985. Tha application was a Siemens motor protector. I called it "remainder saving" because my residuals were division remainders. I don't still have access to that code, so I can't say that it is identical. Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Reply by ●June 19, 20122012-06-19
Hi All.
Foolish me. I have been generating my value for 'A' in MATLAB, and did not click to the fact that Robert's code produces a value for 'A' in the required {16,15} format.
Given that 'A' is a very small value (ie 1 - 0.9999) however, {16,15} is not the best format to represent it. Given a fixed word length of 16-bits, the best format (according to MATLAB) is {16,23}.
In the case of my data, the input format is {22,0}. This is a result of performing prior filtering operations on 16-bit ADC values, and wanting to ensure that I maintain full-precision. So if I use an 'A' value whose format is {16,23}, I am going to end up with an accumulator that at the end of each iteration will contain data whose format is {45,23}. This comes about due to the fact that the operation:
A*y[n-1] = {16,23}*{22,0}
= {38,23}
produces an answer with a fraction length of 23, and you are having to align values whose format is {22,0} with the answer in order to add them. Providing that storage is extended appropriately, and I shift 23 places instead of 15 as per the original code, there should be no problems right?
Finally, a big thanks to Robert for putting his C implementation in the public domain. As has already been clarified, you can find it here
http://www.dspguru.com/dsp/tricks/fixed-point-dc-blocking-filter-with-noise-shaping
I also found an interesting mini-paper on the subject here (warning link to - PDF):
http://www.digitalsignallabs.com/dcblock.pdf
which includes a signal flow graph that produces an identical difference equation to that used by Robert. Suffice to say that I have learned a great deal from both - and from the replies to date.
Reply by ●June 19, 20122012-06-19
On 6/19/12 8:43 PM, Sinclair wrote:> > Foolish me. I have been generating my value for 'A' in MATLAB, and did not click to the fact that Robert's code produces a value for 'A' in the required {16,15} format.note 1: i have never, before your original post about this Rodney, seen this "{16,15}" notation for fixed point. would you consider using a more conventional convention? just so we know what you're talking about? like Q1.15 or 1Q15 or something that defines what is to the left of the imagined binary point and what is to the right? the two numbers add to be the word width. note 2: fixed-point arithmetic is essentially, *integer* arithmetic with some scaling tossed in. in that code example, i steered clear of any Q1.XX notation and just represented my integers as integers.> > Given that 'A' is a very small value (ie 1 - 0.9999) however, {16,15} is not the best format to represent it. Given a fixed word length of 16-bits, the best format (according to MATLAB) is {16,23}.MATLAB tells you that?? then MATLAB doesn't know what it's talking about. note 3: "A" is what it is. it is *scaled* from 1-pole, it is not 1-pole. nevertheless, the precise value of A is not at all critical. you want the pole to be less than 1, but very close to 1. how close simply determines how close to DC is the bottom bandedge for this high-pass filter we call a DC blocker. as long as pole less than 1 (or A is bigger than 0), this will block DC. the gain at DC is 0. it really doesn't much matter if pole is 0.995 or 0.999 or 0.9999 or whatever it rounds to with a 16-bit representation. there is no need to clutter up your code to give A more bits than that. just make sure A is greater than 0. A is a coefficient, not signal. so quantization error in A only changes the filter parameters (the cutoff frequency might be different than it would be without quantization of the coefficients). quantization error in A does not cause quantization error to the signal. you really don't need extra bits for A unless you wish to put the cutoff frequency ultra-close to DC.> > In the case of my data, the input format is {22,0}. This is a result of performing prior filtering operations on 16-bit ADC values,what prefiltering? why not DC block whatever comes out of the ADC first? (that's the main reason for this DC blocker to exist.) why not integrate the DC blocker in with whatever other filtering you're doing. then it's just one big filter. in the past, i have had to LPF filter (i wanted to approximate an integrator) incoming data *and* i needed to DC block it. the result was a simple 2nd-order BPF with the resonant frequency at a very low value. it still blocked DC (rather than integrate it, which would have been fatal).> and wanting to ensure that I maintain full-precision.well, you can never do that. the idea is to not throw away precision until you have to.> So if I use an 'A' value whose format is {16,23}, I am going to end up with an accumulator that at the end of each iteration will contain data whose format is {45,23}. This comes about due to the fact that the operation:> A*y[n-1] = {16,23}*{22,0} > = {38,23}i really have no idea what your "{xx,yy} format notation means. please consider more conventional notation.> > produces an answer with a fraction length of 23, and you are having to align values whose format is {22,0} with the answer in order to add them. Providing that storage is extended appropriately, and I shift 23 places instead of 15 as per the original code, there should be no problems right? > > Finally, a big thanks to Robert for putting his C implementation in the public domain. As has already been clarified, you can find it here > > http://www.dspguru.com/dsp/tricks/fixed-point-dc-blocking-filter-with-noise-shaping > > I also found an interesting mini-paper on the subject here (warning link to - PDF): > > http://www.digitalsignallabs.com/dcblock.pdf > > which includes a signal flow graph that produces an identical difference equation to that used by Robert. Suffice to say that I have learned a great deal from both - and from the replies to date.note 5: Tim, i just don't know what else to do. everytime this DC-blocking with noise-shaping issue comes up, i *always* point them to your algorithm, *including* writing the code so it should be easy for them to see that your alg is so much simpler and more efficient. then i *always* reiterate that "Tim's method is better." so i repeat, this alg: ______________ /* y[n] = x[n] - Quantize{ w[n] } = x[n] - (w[n] + e[n]) w[n+1] = w[n] + (1-pole)*y[n] */ short x[], y[]; double pole = 0.9999; long w, A, curr_y; unsigned long n, num_samples; A = (long)(32768.0*(1.0 - pole)); w = 0; for (n=0; n<num_samples; n++) { curr_y = (long)x[n] - (w>>15); // quantization happens here w += A*curr_y; y[n] = (short)curr_y; } ______________ is *better* and *more* efficient than the code i supplied to dspguru. that's what you should be using, Rodney. maybe i'll give up plugging your new improved alg, Tim. but it does sorta grate me to see it passed over without comment. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●June 20, 20122012-06-20
robert bristow-johnson <rbj@audioimagination.com> writes:> [...] > short x[], y[]; > double pole = 0.9999; > long w, A, curr_y; > unsigned long n, num_samples; > > A = (long)(32768.0*(1.0 - pole));Given that A is always less than 32768, why store it in a long? -- Randy Yates Digital Signal Labs http://www.digitalsignallabs.com
Reply by ●June 20, 20122012-06-20
>robert bristow-johnson <rbj@audioimagination.com> writes: >> [...] >> short x[], y[]; >> double pole = 0.9999; >> long w, A, curr_y; >> unsigned long n, num_samples; >> >> A = (long)(32768.0*(1.0 - pole)); > >Given that A is always less than 32768, why >store it in a long? >-- >Randy Yates >Digital Signal Labs >http://www.digitalsignallabs.com >If you have to cast it to a long later on.
Reply by ●June 20, 20122012-06-20
"niarn" <niaren9@n_o_s_p_a_m.gmail.com> writes:>>robert bristow-johnson <rbj@audioimagination.com> writes: >>> [...] >>> short x[], y[]; >>> double pole = 0.9999; >>> long w, A, curr_y; >>> unsigned long n, num_samples; >>> >>> A = (long)(32768.0*(1.0 - pole)); >> >>Given that A is always less than 32768, why >>store it in a long? >>-- >>Randy Yates >>Digital Signal Labs >>http://www.digitalsignallabs.com >> > > If you have to cast it to a long later on.That doesn't justify storing it as a long. -- Randy Yates Digital Signal Labs http://www.digitalsignallabs.com
Reply by ●June 20, 20122012-06-20
robert bristow-johnson <rbj@audioimagination.com> writes:> [...] > so i repeat, this alg: > > ______________ > > /* > > y[n] = x[n] - Quantize{ w[n] } > > = x[n] - (w[n] + e[n]) > > w[n+1] = w[n] + (1-pole)*y[n] > */ > > short x[], y[]; > double pole = 0.9999; > long w, A, curr_y; > unsigned long n, num_samples; > > A = (long)(32768.0*(1.0 - pole)); > w = 0; > > for (n=0; n<num_samples; n++) > { > curr_y = (long)x[n] - (w>>15); // quantization happens here > w += A*curr_y; > y[n] = (short)curr_y; > } > > ______________ > > is *better* and *more* efficient than the code i supplied to dspguru. > that's what you should be using, Rodney.I would agree it is more efficient. I don't agree it is "better" in all ways. Q15 is not the best use of 16 bits to represent A since A has lots of zeros to the right of the binary point. Something like Q19, etc. (depends on the value of A) would be more accurate. The alternative is to subtract the y[n-1] term separately and add the term p*y[n-1], where representing p by a Q15 is a very good use of all 16 bits. However, grouping your (1-p) this way does save computations. -- Randy Yates Digital Signal Labs http://www.digitalsignallabs.com
