fixed FFT point implementation woes

Hi,
   I am trying to create a mex file that does fixed point FFT (512
point?). I am generally ok with creating a C code in single precision
to compute the the result and have ensured that it works fine as i get
the correct result. (i read it back and plotted it in matlab.)
    the problem comes in when i try to convert it to fixed point. My
input data is 24 bits integer. In addition, I take the sine and cosine
twiddles to also be 24 bits... I do the cos(2*pi*n/N) and sin(2*pi/N)
for n from 0 to 255 the multply the result by ((2^23)-1) and store it
in a look up table.

    What i really dun quite understand is how the scaling factors
work. In FFT books, they generally say that the output from each scale
should maybe be scaled by 1 bit to pre-empt bit growth problems. Also,
for unscaled bit-growth, it generally grows by log2(N) for an N point
FFT. But how does that work?
    What i see here is the following:

1) 24 bit input data. If i dun convert it to anything, it becomes Q.24
format right?
2) 24 bit input data x 24 bit twiddle factor (Q.24 also). this becomes
48 bits. to prevent overflow, i scale by 24 say.

But this is too much and the output frequency spectrum becomes a mess.
I would forsee plenty of precision loss. If i scale by 1 from the
output of each stage, then the FFT gets clipped.

any ideas?what did i go wrong here?  :(

     could someone kindly help me with this? thanks for your help.

On Aug 4, 4:07=A0am, chrisde...@gmail.com wrote:
> Hi,
> =A0 =A0I am trying to create a mex file that does fixed point FFT (512
> point?). I am generally ok with creating a C code in single precision
> to compute the the result and have ensured that it works fine as i get
> the correct result. (i read it back and plotted it in matlab.)
> =A0 =A0 the problem comes in when i try to convert it to fixed point. My
> input data is 24 bits integer. In addition, I take the sine and cosine
> twiddles to also be 24 bits... I do the cos(2*pi*n/N) and sin(2*pi/N)
> for n from 0 to 255 the multply the result by ((2^23)-1) and store it
> in a look up table.
>
> =A0 =A0 What i really dun quite understand is how the scaling factors
> work. In FFT books, they generally say that the output from each scale
> should maybe be scaled by 1 bit to pre-empt bit growth problems. Also,
> for unscaled bit-growth, it generally grows by log2(N) for an N point
> FFT. But how does that work?
> =A0 =A0 What i see here is the following:
>
> 1) 24 bit input data. If i dun convert it to anything, it becomes Q.24
> format right?
> 2) 24 bit input data x 24 bit twiddle factor (Q.24 also). this becomes
> 48 bits. to prevent overflow, i scale by 24 say.

remember that when you do "fixed point" arithmetic, you are really
doing INTEGER arithmetic where *you* are knowledgable about where the
binary point lies.  now i would expect your data is instead Q1.23 so
that there are 23 fractional bits and on bit left of the binary point
for sign.  the implied range is then from -1.000000 oto +0.99999988


>
> But this is too much and the output frequency spectrum becomes a mess.
> I would forsee plenty of precision loss. If i scale by 1 from the
> output of each stage, then the FFT gets clipped.
>
> any ideas?what did i go wrong here? =A0:(

you might need to consider clipping in each pass of the FFT.  perhaps
you need to divide by 2 in the output of each butterfly, then your DFT
has a 1/N factor in front of the summation.  but then in the iDFT
there is no divide by two and your data amplitude grows.  or perhaps
you split the difference and divide by two only in every other pass of
the FFT.  then your DFT definition (and the iDFT) both have 1/sqrt(N)
in the definition.  this sorta thing you have to check out with the
nature of your data and see what works best.

or you can try "block floating-point" arithmetic.

r b-j

robert bristow-johnson wrote:
(snip)

> remember that when you do "fixed point" arithmetic, you are really
> doing INTEGER arithmetic where *you* are knowledgable about where the
> binary point lies.  now i would expect your data is instead Q1.23 so
> that there are 23 fractional bits and on bit left of the binary point
> for sign.  the implied range is then from -1.000000 oto +0.99999988

That is pretty much true for assembly and C programming.

There are programming languages that do fixed point arithmetic
where the compiler keeps track of the radix point.
(Not necessarily binary.)  You still have to be a little careful
to know where overflow might occur.

-- glen

chrisdekoh@gmail.com wrote:
(snip)

>     the problem comes in when i try to convert it to fixed point. My
> input data is 24 bits integer. In addition, I take the sine and cosine
> twiddles to also be 24 bits... I do the cos(2*pi*n/N) and sin(2*pi/N)
> for n from 0 to 255 the multply the result by ((2^23)-1) and store it
> in a look up table.

>     What i really dun quite understand is how the scaling factors
> work. In FFT books, they generally say that the output from each scale
> should maybe be scaled by 1 bit to pre-empt bit growth problems. Also,
> for unscaled bit-growth, it generally grows by log2(N) for an N point
> FFT. But how does that work?

First you have to be able to do the multiply.  If you multiply two
numbers, each with 23 bits after the binary point the product has
46 bits after the binary point.  You probably don't want all those
bits, but C doesn't have a good way to get the right bits.

If you have the bits, it is best to expand by one bit for each
stage, such that you avoid overflow but still keep all the
significant bits.  You can figure out at the end which bits
you really want.  Otherwise, you need to know more about the
actual data.  In most cases, bit growth will be somewhat
slower, except for the DC (0) term, which can fairly easily
accumulate large values.

-- glen