Altera FFT - output scaling. How to do this correctly?

[Altera_Forum]

Hi, 

 I do not get correct scaled results if I do everything the FFT user guide recommends.   I have instantiated a 1024 points FFT, with single output and burst architecture. The data width and twiddle factor precision is 16bits.   I input 1024 times 16bit signed integers from the ADC to the 16bit real-input of the FFT module. (The ADC has a +/-2.5V input) To the imaginary input of the FFT module I apply 16bit zeros. This should be correct, shouldn’t it? The FFT module outputs 1024 pairs of 16bit real and 16bit imaginary part.   First question: Are the real and imaginary 16bit values the mantissa in the sense of a floating point mantissa (first bit = sign, second bit = 1.digit in front of the radix point, remaining 14bits = decimal places) OR do I have to interpret these values as fixed point values.   I have a sample frequency of 3kHz and used a very low frequency of 1Hz and 9Hz to verify the scaling to not have additional falsifications, caused by the rectangular window I used for the first shot. I input samples of an ideal sine waveform with different amplitudes and the same frequency, BUT whether I scale the output back with the given exponent value or not, the calculated magnitude is not the expected one. For example: For a 9Hz sine with an amplitude of +/-2,3V which is equal to +30146/-30147 decimal (+/-2.5V ADC input) I got the following real and imaginary parts for the corresponding frequency.   Real: 3384 (0x0D38) Imaginary: -14391 (0xC7C9)   I would expect a magnitude in the range of 30146 decimal, but using these values the magnitude is 14783. Using the exponent makes the results too big and I see no relation to input signal with such big numbers. How can such a large output (31bits in my case, max exponent = 15bits) provide unity gain, if the input is much smaller (16bits). Just choosing the lower 16bits of the shifted/scaled values can not be correct, because of the sign. I am a bit confused.      Can somebody give me a hint? An example calcualtion would also be a great help.

[Altera_Forum]

First, your imaginary value being set to zero is ok. It simply means you are dealing with real input and so expect two lines in the spectrum for each input frequency. 

 second, you need to think of fixed point only(unless your fft is not).   third, you should get two symmetrical bins at nearly your amplitude value.    Your figure suggest a multiplication by 2(2 x 14783 =~ 30146).   About scaling: (why you say input is 16 bits, output is 31 ?) If altera says divide output by 2^exp i.e. multiply by 2^-exp then: division: use shift towards LSB extending sign bit and discarding those LSBs that drop off as you shift.   for mult:shift magnitude bits towards sign bit inserting leading zeros.

[Altera_Forum]

I get two symmetrical lines in the spectrum. Fine, so far. The frequency values are the expected ones.  

 second, you need to think of fixed point only(unless your fft is not). Since I'm using burst architecture, I can not chose fixed or floating point representation like with the variable streaming architecture. I assume the output is fixed point by default.    Another issue is, that referring to Alteras documentation, it seems to me that Altera does not differ between integer and fixed-point data format. Theortically that is not the same, but if +0,9999999 (fixed-point) or +32767 (integer) is just a question of interpretation, istn't it? The digital number I give the FFT module is the same, in this case 0x7FFF (16bit).    third, you should get two symmetrical bins at nearly your amplitude value.    I do not get neraly identical bins at neraly my amplitude value. I monitored the direct output signals of the FFT Megacore instance with the SignaltapII and I got 0x0D38 (real) and 0xC7C9 (imag) at the output for the frequency part I'm interested in. The exponent output shows 0x36 for the entire block. So this is not a bug in the wrapper I wrote around the FFT instance. That's why I asked for the input format, too. Maybe the input is not correct.   Input example: For a 9Hz sine signal, sampled with 3kHz and an amplitude of 2,3V the maximum input amplitude integer (16bit signed) is 30147 decimal. This should be correct.    about scaling: (why you say input is 16 bits, output is 31 ?) if altera says divide output by 2^exp i.e. multiply by 2^-exp then:  I said the output is 31 bit, because I referred to the FFT user guide, where it is mentioned that you have to extend your output by the maximum number of shifts, that is possible... in my case 15bits. Since I'am using the NiosII processor to process the results I can use the entire 31 bits for a high resolution, if I want. But where is the difference between these:   --output unscaled 0x7FFF  -- output scaled: shifting e.g. four times left (mult) 0x7FFF0 --> 0x7FFF when taking the upper 16bits    I win no new information with this multiplication, do I?    for mult:shift magnitude bits towards sign bit inserting leading zeros.   Should I scale the magnitude? I think you mean I should scale real and imginary parts. We should not mix the terms here. This is confusing.   By the way the spectrum is always mirrored at half the sample rate. The information above fs/2 is redundant and obviously not helpful. Why is this output? I could save half the memory by dropping the data above fs/2. Is this correct?

[Altera_Forum]

--- Quote Start ---   By the way the spectrum is always mirrored at half the sample rate. The information above fs/2 is redundant and obviously not helpful. Why is this output? I could save half the memory by dropping the data above fs/2. Is this correct?  --- Quote End ---     yes and no. If you don't want phase information then yes.  amplitude is symmetrical around dc(and around Fs/2 as well) but phase may not be. The issue becomes relevant if you frequency shift your sinusoid to higher point e.g. move from dc to a carrier at 1MHz then you get two side bands, one from each side of your baseband (around dc or each side of Fs/2).     --- Quote Start ---   Another issue is, that referring to Alteras documentation, it seems to me that Altera does not differ between integer and fixed-point data format. Theortically that is not the same, but if +0,9999999 (fixed-point) or +32767 (integer) is just a question of interpretation, istn't it? The digital number I give the FFT module is the same, in this case 0x7FFF (16bit).   --- Quote End ---     integers are represented as fixed point. If you think of .1, .2 ...etc then they are first scaled to integers to be represented on the binary ladder of 16bits. Scaling is up to your flavour,possible signal range and some other factors.     --- Quote Start ---   --output unscaled 0x7FFF  -- output scaled: shifting e.g. four times left (mult) 0x7FFF0 --> 0x7FFF when taking the upper 16bits    I win no new information with this multiplication, do I?      --- Quote End ---     The issue of scaling makes sense only if there is room to scale, otherwise altera got it wrong !  Frankly I don't understand why altera left scaling for the field engineer.     --- Quote Start ---   for mult:shift magnitude bits towards sign bit inserting leading zeros.   Should I scale the magnitude? I think you mean I should scale real and imginary parts. We should not mix the terms here. This is confusing.     --- Quote End ---     magnitude bits as opposed to sign bit, context obvious     --- Quote Start ---   I do not get neraly identical bins at neraly my amplitude value. I monitored the direct output signals of the FFT Megacore instance with the SignaltapII and I got 0x0D38 (real) and 0xC7C9 (imag) at the output for the frequency part I'm interested in. The exponent output shows 0x36 for the entire block. So this is not a bug in the wrapper I wrote around the FFT instance. That's why I asked for the input format, too. Maybe the input is not correct.   about scaling: (why you say input is 16 bits, output is 31 ?) if altera says divide output by 2^exp i.e. multiply by 2^-exp then:   I said the output is 31 bit, because I referred to the FFT user guide, where it is mentioned that you have to extend your output by the maximum number of shifts, that is possible... in my case 15bits. Since I'am using the NiosII processor to process the results I can use the entire 31 bits for a high resolution, if I want. But where is the difference between these:   --- Quote End ---     The exponent has to be a small value and having 15 bits is unreasonable. Do you mean maximum exponent value can reach 2^15 -1 = 32767  Do you then mean you shift 32767 bits. Altera IP gurus must limit the exponent to a realistic value in harmony with output bit width.

[Altera_Forum]

I mean the exponent is 6bits signed and can range from 1 to -15, referring to page A-2 of the FFT user guide.  

 The user guide tells to scale the output with the exponent value of the processed block. Doesn't that mean that I have scale every block of 1024 points I process? But if not, what is the exponent output for?   The question is: Do I have to use the exponent value? And how do I have to apply it exactly. The goal is to get the amplitude spectrum.   Right now I calculate the magnitude this way:   mag = sqrt(16bitreal^2 + 16bitimag^2) ...(the result is too small)   I also tried:   mag = sqrt((16bitreal^2 * 2^(-exponent)) + (16bitimag^2 * 2^(-exponent)) ) ... the result is too large, compared to the input range.   Is this correct?

[Altera_Forum]

A further note on shifting: 

 imagine your data is 15bits  if data = Ox00ff (255) shift 4 bits towards sign bit to become Ox0ff0 (or Ox0ff&four internal bits if available). = 4080 = 16 x 255   or shift 4 bits towards LSB to become Ox000f = 15 =~ 255/16   if your data = Ox7FFF (your example) = -1 shift 4 bits towards sign bit to get Ox7FF0 = -16 or 4 bits towards LSB to get Ox7FFF = -1 (this is the odd case, woudn't divide)

[Altera_Forum]

--- Quote Start ---  I mean the exponent is 6bits signed and can range from 1 to -15, referring to page A-2 of the FFT user guide.    The user guide tells to scale the output with the exponent value of the processed block. Doesn't that mean that I have scale every block of 1024 points I process? But if not, what is the exponent output for?   The question is: Do I have to use the exponent value? And how do I have to apply it exactly. The goal is to get the amplitude spectrum.   Right now I calculate the magnitude this way:   mag = sqrt(16bitreal^2 + 16bitimag^2) ...(the result is too small)   I also tried:   mag = sqrt((16bitreal^2 * 2^(-exponent)) + (16bitimag^2 * 2^(-exponent)) ) ... the result is too large, compared to the input range.   Is this correct?  --- Quote End ---     yes you apply exponent per block as it varies. Your scaling formula is correct but how do you implement in hardware. do you apply shift or use multipliers in hardware? or do you get results wrong on paper anyway?

[Altera_Forum]

When you say exponent is Ox36 then make sure this is not used as 36 but as -10 in 2's complement...

[Altera_Forum]

Now it seems to work. But I don't know the correct reason why. Let me describe. 

 First of all I use a Nios II processor to handle the data of the FFT. I don't scale anything in hardware. Performance is not a problem right know.    I wrote a wrapper around the FFT instance so that the real and imaginary outputs are merged to a 32bit bus and the exponent can be read separatetly.    here my C-code for the Scaling:   samples_u32arr[] is an static array holding all 1024 32bit values (16bit real+16bit imag) i read from the fft module. number_of_samples = 1024.   for(samplecounter_u16 = 0; samplecounter_u16 < number_of_samples; samplecounter_u16++) { real_16arr[samplecounter_u16] = (int16_t)samples_u32arr[samplecounter_u16] & 0xffff; imag_16arr[samplecounter_u16] = (int16_t)((samples_u32arr[samplecounter_u16] >> 16) & 0xffff);     real_32arr[samplecounter_u16] = ((int32_t)real_16arr[samplecounter_u16]) << -exponent_s8; imag_32arr[samplecounter_u16] = ((int32_t)imag_16arr[samplecounter_u16]) << -exponent_s8;     real_16 = (int16_t)((real_32arr[samplecounter_u16] >> 9) & 0xffff); imag_16 = (int16_t)((imag_32arr[samplecounter_u16] >> 9) & 0xffff);   result_u16arr[samplecounter_u16] = (uint16_t)sqrt(real_16*real_16 + imag_16*imag_16); }   The strange thing is, that I have to shift the data right 9 times, after I shifted it left with the value of the exponent, which is changing with the input data. I don't know why exactly the number 9!!! But if I do shift 9 times right and choose the lower 16bit the output amplitude corresponds to the input amplitude, with small deviations.   I tested the scaling with a sine and a rectangular waveform with different amplitudes at the same frequency of 9Hz (3kHz sample rate). So there are no additional deviations caused by the window function (in this case just a rectangle).

[Altera_Forum]

I am not sure about the 9 bits shift but here is my suspicion: 

If your fft output bitwidth is 25 bits while your input bitwidth is 16 bits then I guess the IP may introduce 2^9 scaling to fill up all output bits. Thus you need to divide back by 2^9. You can just take the 16 MSBs instead of shifting right then taking lower 16 bits.

[Altera_Forum]

--- Quote Start ---  I am not sure about the 9 bits shift but here is my suspicion: If your fft output bitwidth is 25 bits while your input bitwidth is 16 bits then I guess the IP may introduce 2^9 scaling to fill up all output bits. Thus you need to divide back by 2^9. You can just take the 16 MSBs instead of shifting right then taking lower 16 bits.  --- Quote End ---     My output is not 25bits wide. It is 16bits + a maximum of 15bits. So chosing the MSBs only is not a good idea

[Altera_Forum]

Ignore my above remarks since your fft output is 16 bits. 

 I understand from your equations that you are converting the 16 bit output into 32 bit bus then doing the two stage shift(by exponent then by 9). Your work-around of 9 bits shift then taking 16 LSBs doesn't make sense to me.   To avoid calculation problems of software Let us go back to your simple example: You said for your input sine of 30146 amplitude you get fft output as:   Real: 3384 (0x0D38) Imaginary: -14391 (0xC7C9)   fair enough. what is your exponent value in above example? It should be -1 so that output scales to 2^ -(-1) = 2   otherwise there is something wrong.  

[Altera_Forum]

--- Quote Start ---  Your work-around of 9 bits shift then taking 16 LSBs doesn't make sense to me.   --- Quote End ---     This does not makes sense to me, too. That's why I'm not happy with this solution.     --- Quote Start ---  To avoid calculation problems of software Let us go back to your simple example: You said for your input sine of 30146 amplitude you get fft output as:   Real: 3384 (0x0D38) Imaginary: -14391 (0xC7C9)   fair enough. what is your exponent value in above example? It should be -1 so that output scales to 2^ -(-1) = 2   otherwise there is something wrong.    --- Quote End ---     The exponent within this example is -10. Yes, there seems to be something wrong, but what is the question.

[Altera_Forum]

I should add, that I have monitored the exponent output of the FFT IP-core. The value I read is correct.  

Then the question raises: Is the input correct?

[Altera_Forum]

You siad your inputis 30146, I assume you checked that in signaltap. 

 This matlab code can be used as reference:   % generate a sinusoid at 9Hz, Fs at 3Khz x = round(30146*cos(2*pi*[0:1023]*9/3000));   % fft scaled for unity gain y = 1/512*fft(x);   plot(abs(y)); max(abs(y))

[Altera_Forum]

Mmmmhhh I have no Mathlab at my company. But I can check this at home this evening.  

 Referring to the input problem. I simulated the IP-core (this time with a quad-output architecture, for the reason of performance) with my test data, I attached it in form of two text files which are required by the simulation.  Within the simulation, I verfied again that the sine amplitude is 30146 (first screenshot) and that the exponent is also -10 (second screenshot). (I also checked this in signaltapII)    The real bin# 3 (=> 3x 2.92Hz frequency resolution = about 9Hz) has a value of 0x0D3A (with single output architecture this was 0x0D38). The imaginary bin# 3 has a value of 0xC7C7 (with single output architecture this was 0xC7C9). Look at the screenshots I attached.    So the results depend on my input OR on the FFT IP-core itself and not on my processor system or my wrapper.    What could be the reason for this behavior?    P.S.: Thank KAZ for spending time on helping me.

[Altera_Forum]

I checked your inputs in matlab. All is ok if you scale the output by 2. 

for example amplitude given by matlab fft for your input is 29569 which is pretty close to double your first example result.   Now you proved that the problem is indeed to do with fft IP scaling, that is what you started with anyway(-10 should be -1). You now either cautiously use your workaround of rescaling (exp=exp+9) and wait for other people on the form to respond with their experience or ask Altera for help.

[Altera_Forum]

--- Quote Start ---  I checked your inputs in matlab. All is ok if you scale the output by 2. for example amplitude given by matlab fft for your input is 29569 which is pretty close to double your first example result.   Now you proved that the problem is indeed to do with fft IP scaling, that is what you started with anyway(-10 should be -1). You now either cautiously use your workaround of rescaling (exp=exp+9) and wait for other people on the form to respond with their experience or ask Altera for help.  --- Quote End ---     I will ask Altera, because waiting is not a solution.    The last idea I have refers to the input signal. Are x(n), the samples, really the real bins OR are they the magnitude? If it would be the magnitude, the real and imaginary bins would not be correct. But then another question would raise... Where is the phase information?!    By the way kaz, did you ever implement an Altera FFT? How does it work at your project?

[Altera_Forum]

First I did altera fft in 2001 !! no problem then except that the output bit resolution wasn't well done by altera, later they improved that. I did the scaling manually by bit shifting in vhdl as speed was crucial and couldn't depend on software speed. My exponent was correct at that time. 

 your input(x): xreal(n) is only real, ximag is zero so the magnitude of your input is same whether you view your input as x=x(n) or x=complex (xreal,ximag).  I don't get what you say about magnitude?? but don't mix up between magnitude of a complex signal and amlpitude of a sinusoid. The maximum |magnitude| of your input(if you think of it as complex) is the amplitude of your sinusoid.   your output(y): y is complex i.e. (y= complex(yreal, yimag)) and both elements together determine amplitude and phase(abs(y),angle(y)) in frequency domain. Your phase info is not relevant here, phase is only relevant when comparing two sinusoids at same frequency e.g. ouput/input of a system then phase tells you what this system is doing to phase at various frequencies.

[Altera_Forum]

--- Quote Start ---  The maximum |magnitude| of your input(if you think of it as complex) is the amplitude of your sinusoid.   your output(y): y is complex i.e. (y= complex(yreal, yimag)) and both elements together determine amplitude and phase(abs(y),angle(y)) in frequency domain.  --- Quote End ---     That's exactly what I mean and that's how I learned/understood it. Ok I understood everthing correct, the input signal seems to be correct (you also looked at it), and the simulation outputs the same results as real module. Than there must be some problem with the IP-core. I opened a service request. Maybe this will result in a solution of the problem.    So far I'm using my workaround with the additional shifting by 9.

[Altera_Forum]

Hello Kaz, 

 I got the answer from Altera support. The guy who answered me, said, that my workaround can not work for every signal input, as you said and as I thought. Further he said I should find the maximum of my input samples, either by scaling up to a certain maximum or by just looking htrough the samples. Then if I have read the samples at the FFT output I have to calculate the magnitude of the vector (sqrt(RE²+IM²)) of each bin, serach through the block for the maximum magnitude nad compare it to the maximum input amplitude. With the factor between these two maxima I can scale up every calculated magnitude.    This is really frustrating and also slow, because I first have to calulate the magnitude before scaling up. Doing this in hardware is not so easy for me, because i only now the basics of VHDL. It's just more complicated.    Thanks for the help kaz.