Hi David,

Before I answer your question(s), could I just check that we are both using the same terminology?

The fir_f16 and rfft_f16 library functions operate on data of type fract16. The fract16 data type is used to represent 16-bit fixed-point data in the range from 0.999996 to -1.0. It is different from floating-point because by convention the decimal point always appear to the left of the left-most binary digit.

Compare this to a regular 16-bit integer value which is used to represent values in the range 32767 to -32768; here the decimal point is always to the right of the right-most digit. So 0x7FFF represents the integer value 32767 and the fract16 value 0.999996 (also, for example, 0x4000 is the integer value 16384 or the fract16 value 0.5).

So when you say "in matlab, my fft's go up to values of about 100..." do you mean the floating-point value "100.0" or do you mean the integer "100" (which has the effective value of 0.00305 when assigned to a fract16)?

Let me know, okay?

Hi There,

I think that I can answer your question, but before I do could I first point out a couple of problems with your code snippet? The first problem is that you should be using the library function twidfftrad2_fr16 instead of twidfft_fr16; the function twidfft_fr16 will generate a twiddle table for the legacy functions cfftrad4_fr16, ifftrad4_fr16, and rfftrad4_fr16. By using this function in your application you should notice that the imaginary component of the FFT that you generate will not be correct. So you should use:

complex_fract16 twiddle[FFTSIZE/2];

twidfftrad2_fr16 (twiddle,FFTSIZE);

instead of:

complex_fract16 twiddle[(3*FFTSIZE)/4];

twidfft_fr16 (twiddle,FFTSIZE);

The second problem that I can see relates directly to the question that you are asking. The 6'th argument of rfft_fr16 is defined as a pointer to an integer, but you have specified it as '0'; correcting this error will not modify the FFT that rfft_fr16 generates but it will help you interpret the FFT. For example, instead of:

rfft_fr16(in,out, twiddle, 1, FFTSIZE, 0, 2);

use the following:

int block_exponent;

rfft_fr16(in,out, twiddle, 1, FFTSIZE, &block_exponent, 2);

Now to answer your question:

During the computation of an FFT there is a risk with fixed-point arithmetic that intermediate computations will overflow and generate saturated results, and this will obviously affect the quality of the generated output. To avoid the risk of overflow, the FFT function will have to scale intermediate results. For a given FFT, there are 'log2(FFTSIZE)' intermediate set of results, and one option is to unconditionally divide each set of intermediate results by 2 - this option is called static scaling and will result in an FFT that has been scaled by log2(FFTSIZE). Static scaling is cheap to implement on the Blackfin architecture but it can lead to a loss of precision particularly for larger FFT sizes. An alternative scaling method is known as dynamic scaling, and under dynamic scaling intermediate results are inspected and are only scaled if there is a risk that these intermediate results may cause overflow. Dynamic scaling is slower than static scaling but generally it results in a more accurate FFT.

The user selects which scaling method should be used via the 7'th argument of the rfft_fr16 function; this argument is called 'scale_method'. If the argument is set to either 0 (zero) or 1, then the function will perform static scaling; if it is set to 2 (as in your example) the function will perform dynamic scaling. If one is sure that that the FFT function will not generate a saturated result, then one can set the scale_method argument to 3, in which case the function will never scale intermediate results and there will be no loss of precision in the generated output.

To properly interpret the output from the rfft_fr16 function, one needs to know how many times the library function has scaled the intermediate results. This is particularly important for dynamic scaling and this is where the 6'th argument (known as the block exponent argument) is so important.

Before returning to your program, the FFT function will set the block exponent argument to the number of intermediate results that it has scaled (i.e. divided by 2). If static scaling has been selected, then the block_exponent argument will be set to log2(FFTSIZE), if no scaling has been selected then the block_exponent argument will be set to zero. If dynamic scaling has been specified then the block_exponent argument will be set to the number of times each set of intermediate results have be scaled - this will be some value between zero and log2(FFTSIZE).

So to check that MATLAB and rfft_fr16 generate similar FFTs, you should change the following code (at the end of your snippet):

for(i=0;i<FFTSIZE;i++)

{

real[i]=fr16_to_float(realf[i]);

imag[i]=fr16_to_float(imagf[i]);

}

to:

for(i=0;i<FFTSIZE;i++)

{

real[i]=fr16_to_float(realf[i])*ldexpf(1.0F,block_exponent);

imag[i]=fr16_to_float(imagf[i])*ldexpf(1.0F,block_exponent);

}

PS: I see that you are using version 5.0 of VisualDSP++, but I have a word of caution to anyone who is interested in this reply and is only

using Version 4.5 of VisualDSP++ (or earlier). The cfft_fr16, ifft_fr16, and rfft_fr16 library functions in the earlier releases only support static scaling and do not return a block exponent. The interface to these functions has been modified for version 5.0 of VisualDSP++.