Convert to higher sample rate

glen herrmannsfeldt <gah@ugcs.caltech.edu> writes:

> Randy Yates <yates@digitalsignallabs.com> wrote:
>
> (snip)
>
>> I'm not sure we agree on the real issue here. Perhaps that's 
>> the reason for some of the confusion.
>  
>> FACT: Both large numbers (e.g., those near 1) and small numbers (those
>> near 0) have the same number of significant digits in IEEE-754 floating
>> point format.
>> 
>> FACT: The ONLY purpose the "hidden bit" serves is to allow one more
>> significant digit than you otherwise would have with a direct
>> representation.
>  
>> FACT: The number of significant digits is directly related to the
>> "relative error" (see section 1.2 in [higham]). 
>  
>> So.., while your statement above is true, I would say "So what?"! 
>> The number of significant digits in the second value (53) is 
>> greater than the number in the first value (39). So if you're 
>> using a representation that has less than 53 bits but more than 
>> 38 bits, you'll preserve the first but not the second. Of course.
>  
>> That sort of analysis doesn't prove that numbers near one are less
>> accurate than those near zero. In fact the IEEE 754 format ensures both
>> are represented with the same accuracy (or equivalently, the same
>> relative error).
>  
>> As Higham shows so well in section 1.7, the problem is that SUBTRACTIVE
>> CANCELLATION MAGNIFIES RELATIVE ERRORS. 
>  
>> In fact the same magnification could happen with numbers near zero.
>
> I believe all those are true, but the original question was for 
> the FFT, and I don't believe that is enough to explain the error
> model for the FFT.

I agree with that, Glen. This section of the thread had degenerated (at
least for me) into discussing with Robert why the hidden bit has nothing
to do with the cosine problem, or more specifically, with the error in 1
- cos(theta) when theta is very small.

> I heard from a reliable source (but I forget which one) that the
> FFT rounding is much better than the DFT rounding. All the intermediate
> roundings tend to average out in a way that one big rounding doesn't.

I have a vague intuition that is true, due to the nature of the twiddle
factors in two-by-two butterflies.

> As far as I know, though, it isn't easy to show in general how
> close one comes to the right answer. 

I can believe that!
-- 
Randy Yates
Digital Signal Labs
http://www.digitalsignallabs.com

On 4/18/14 6:09 PM, Randy Yates wrote:
> glen herrmannsfeldt<gah@ugcs.caltech.edu>  writes:
>
...
>>
>> I believe all those are true, but the original question was for
>> the FFT, and I don't believe that is enough to explain the error
>> model for the FFT.
>
> I agree with that, Glen. This section of the thread had degenerated (at
> least for me) into discussing with Robert why the hidden bit has nothing
> to do with the cosine problem, or more specifically, with the error in 1
> - cos(theta) when theta is very small.

let's review what was said:

someone said:
>> When you get to the point where sin(2*pi/(large number)) is less
>> than epsilon, then you do have to be careful. We are not so close
>> to that with IEEE double.

then Randy said:

> Seem's like you'd be in trouble at some point when sin(blah) is
> anywhere near epsilon.
>
> sin(blah)is about 4E-10 on my calc. Isn't epsilon around 1E-15 for
> doubles?

and i said:

the real problem is with the cosine of teeny values.

now, in modern times, the cosine problem might not be as much of a 
problem for the FFT as my alarmism implies, because we don't usually 
generate our coefficients with trig identities equivalent to

     e^(j*2*pi*(k+1)/N) =  e^(j*2*pi/N) * e^(j*2*pi*k/N)

now, we'll just load the twiddle coefs in from somewhere else or we'll 
take the time and generate nice N/2 coefs from the standard library (and 
the FFT won't be too fast).  but in the olden days, we would use 
something like the above to step through and generate the coefs, and 
this cosine problem is certainly a problem in that case.

back to what Randy said about sin(Tiny) as being the problem.  i am 
saying that you'll run into trouble when Tiny is larger (like the 
sqrt(Tiny)) with the cosine function of the same floating-point word 
width.  and the reason is obvious and i spelled it out.

you run into this cosine problem in IIR filtering with coefficient 
calculation when the resonant frequency is much much lower than Nyquist. 
  and you run into this cosine problem when you plot the frequency 
response of an FIR or IIR filter at frequencies much much lower than 
Nyquist.   this is because, in both of those problems there is no 
sin(omega) because the equations are even symmetry in omega and there is 
only cos(omega).

now the problem is that because of precision limitations, your 
calculated cos(omega) will not be precisely cos(omega) but will have 
some error.  so, even though it's not precisely cos(omega) for the omega 
you specified, it *is* the cosine of some *other* omega (call it "omega1").

so  cos(omega1)  =  cos(omega) + rounding_error

so then the filter response that you're plotting won't give you

    H( e^(j*omega) )

it gives you H( e^(j*omega1) ).

and as omega gets very small (and, with single precision, you run into 
trouble in the bottom one or two octaves in 48 kHz audio, and in 96 kHz 
audio you run into it at one octave higher) the, difference between 
omega and omega1 is significant when single-precision float is used.

now, here is the solution to the cosine problem:  everywhere that you 
see cos(x), if x is considered to get small, what you want to do is 
replace every occurance of cos(x) with 1 - 2*(sin(x/2))^2 and rework all 
of your equations so that they are functions of sin(x/2) rather than 
functions of cos(x).

now, Randy, if you say that this isn't a problem and that this problem 
with cos() does not occur before an equivalent precision problem with 
sin(), i'm just gonna throw my hands up in the air.  i've spelled it out 
as best as i can.  cos() of tiny numbers is a bigger problem than sin() 
of tiny numbers, all things being equal.  with fixed-point it is true 
because sin(x) differs from zero by about x.  but cos(x) differs from 1 
by about x^2.  and a much larger theta will your cos(theta) be 
indistinguishable from cos(0) than would sin(theta) be indistinguishable 
from sin(0).  and that's *fixed* point.

with floating point, there is a second reason why cos() of small values 
is more of a problem than sin() of the same small values and that is 
because floating point maintains good values of sin(theta) even as theta 
gets small.  but floating point does no better than fixed point (and 
even worse if you consider the cost in bits of the exponent) of the same 
bits.  your little (theta^2/2) will fall of the edge really damn quick 
and that's all of the information you have about theta.  so your cos() 
becomes the cos() of a *different* theta.

>> I heard from a reliable source (but I forget which one) that the
>> FFT rounding is much better than the DFT rounding. All the intermediate
>> roundings tend to average out in a way that one big rounding doesn't.

i think if you computed the naive DFT with a double-wide accumulator 
before rounding and storing the single precision result, that this would 
have far less quantization error (because it doesn't grow with N) than 
the basic Cooley-Tukey FFT (which grows with N for fixed-point and with 
log(N) for float) with each pass data stored as single precision.

but if your naive DFT has a single-wide accumulator (so there would be 
roundoff after ever MAC), then the quantization error would also grow 
with N.


-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

Robert,

This (sub-)thread is getting damn-near impossible to track because of
all the he-said/she-said's and the multitude of possible issues and
their causes, so I'm not going to take the time address everything you
said. But I will address this, which I think is the (or a least a key)
disagreement between us.

robert bristow-johnson <rbj@audioimagination.com> writes:
> [...]

> back to what Randy said about sin(Tiny) as being the problem.  i am
> saying that you'll run into trouble when Tiny is larger (like the
> sqrt(Tiny)) with the cosine function of the same floating-point word
> width.  and the reason is obvious and i spelled it out.

Are you talking about the 

  0.000000000000001xxxxxxxxxxxxxxxxxxxxxxxxvvvvvvvvvvvvvv

versus 

  0.111111111111111xxxxxxxxxxxxxxxxxxxxxxxxvvvvvvvvvvvvvv

stuff? 

If so, you haven't spelled out shit. I don't agree with you on this, and
_I_ have spelled out for _you_ why, by using either direct mathematical
specifics or references to them (via Higham). If you simply want to use
proof by assertion, go right ahead - I will stop beating my head against
the wall.

If you want to move forward, see my previous response, which hasn't
changed. For the sake of clarity, I'll summarize here:

  The representations of sin(eta) and cos(eta) in IEEE 754 FP, when eta
  is very close to zero, have the same (or very similar - same order of
  magnitude) relative error, or equivalently, the same number of
  significant digits.

  Also, the problem of subtractive cancellation can occur (NOT "DOES
  occur in your pet problem") with either representation.

By the way, I think you may be confusing issues with the _computation_
of sin(eta) and cos(eta) with their representation. What I'm talking
about is their representation, i.e., if you knew the exact result of
sin(eta) and cos(eta) to a hundred decimal places and cast those to
IEEE-754 FP, what would the error be?

And by the way, I could give a rat's ass at this point about what
applications (e.g., IIR filtering) cause the proported problem to
appear. You need to show what the problem IS, specifically, without
waving a hand, and preferably with a little mathematics, and without
asserting "this comes up in XYZ application" (which is really
secondary at this point of analysis).

And stay on target. You wander like the steering on a late-60's Ford.
-- 
Randy Yates
Digital Signal Labs
http://www.digitalsignallabs.com

PS: Robert, if you don't have access to the Higham text and want to read
it, let me know and I'll get a scan to you.

--RY

Randy Yates <yates@digitalsignallabs.com> writes:

> Robert,
>
> This (sub-)thread is getting damn-near impossible to track because of
> all the he-said/she-said's and the multitude of possible issues and
> their causes, so I'm not going to take the time address everything you
> said. But I will address this, which I think is the (or a least a key)
> disagreement between us.
>
> robert bristow-johnson <rbj@audioimagination.com> writes:
>> [...]
>
>> back to what Randy said about sin(Tiny) as being the problem.  i am
>> saying that you'll run into trouble when Tiny is larger (like the
>> sqrt(Tiny)) with the cosine function of the same floating-point word
>> width.  and the reason is obvious and i spelled it out.
>
> Are you talking about the 
>
>   0.000000000000001xxxxxxxxxxxxxxxxxxxxxxxxvvvvvvvvvvvvvv
>
> versus 
>
>   0.111111111111111xxxxxxxxxxxxxxxxxxxxxxxxvvvvvvvvvvvvvv
>
> stuff? 
>
> If so, you haven't spelled out shit. I don't agree with you on this, and
> _I_ have spelled out for _you_ why, by using either direct mathematical
> specifics or references to them (via Higham). If you simply want to use
> proof by assertion, go right ahead - I will stop beating my head against
> the wall.
>
> If you want to move forward, see my previous response, which hasn't
> changed. For the sake of clarity, I'll summarize here:
>
>   The representations of sin(eta) and cos(eta) in IEEE 754 FP, when eta
>   is very close to zero, have the same (or very similar - same order of
>   magnitude) relative error, or equivalently, the same number of
>   significant digits.
>
>   Also, the problem of subtractive cancellation can occur (NOT "DOES
>   occur in your pet problem") with either representation.
>
> By the way, I think you may be confusing issues with the _computation_
> of sin(eta) and cos(eta) with their representation. What I'm talking
> about is their representation, i.e., if you knew the exact result of
> sin(eta) and cos(eta) to a hundred decimal places and cast those to
> IEEE-754 FP, what would the error be?
>
> And by the way, I could give a rat's ass at this point about what
> applications (e.g., IIR filtering) cause the proported problem to
> appear. You need to show what the problem IS, specifically, without
> waving a hand, and preferably with a little mathematics, and without
> asserting "this comes up in XYZ application" (which is really
> secondary at this point of analysis).
>
> And stay on target. You wander like the steering on a late-60's Ford.

-- 
Randy Yates
Digital Signal Labs
http://www.digitalsignallabs.com

On 4/19/14 10:17 PM, Randy Yates wrote:
> Robert,
>
> This (sub-)thread is getting damn-near impossible to track because of
> all the he-said/she-said's and the multitude of possible issues and
> their causes, so I'm not going to take the time address everything you
> said. But I will address this, which I think is the (or a least a key)
> disagreement between us.
>
> robert bristow-johnson<rbj@audioimagination.com>  writes:
>> [...]
>
>> back to what Randy said about sin(Tiny) as being the problem.  i am
>> saying that you'll run into trouble when Tiny is larger (like the
>> sqrt(Tiny)) with the cosine function of the same floating-point word
>> width.  and the reason is obvious and i spelled it out.
>
> Are you talking about the
>
>    0.000000000000001xxxxxxxxxxxxxxxxxxxxxxxxvvvvvvvvvvvvvv
>
> versus
>
>    0.111111111111111xxxxxxxxxxxxxxxxxxxxxxxxvvvvvvvvvvvvvv
>
> stuff?
>
> If so, you haven't spelled out shit.

Randy, denial ain't just a river in Egypt.

so break out your MATLAB or Octave or whatever and consider the 
single-precision floating-point case and the cosine function. consider a 
sample rate of 96 kHz.  what are the values of

    96000/(2*pi)*acos(1 - n*2^(-24))

for n = 0, 1, 2, 3, 4, up to about 10.  "acos" is "arccos".

can you see a problem?  if so, can you tell us if the sine function 
suffers a corresponding problem to anywhere close to the same degree?

-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

robert bristow-johnson <rbj@audioimagination.com> writes:
> [...]
> Randy, denial ain't just a river in Egypt.

Robert, darkness ain't just experienced at night.
-- 
Randy Yates
Digital Signal Labs
http://www.digitalsignallabs.com