#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <iostream>
#include <iomanip>
#include "awgn.h"

using namespace std;
// fixed-point header file
// Qm.f
// m: the number of bits for integer
// f: the number of bits for fraction 

typedef signed short fx_t;

#define M 3  // m
#define F 2  // f



#define F_MAX ( ( 1 << (M+F) ) - 1 )
#define F_MIN ( -1 * ( 1<<(M+F) ))
#define Fl2Fx(x) (fx_t) (x >= 0) ? ( (x * (1 << F)) + 0.5 ) : ( (x * (1 << F)) - 0.5 )    //floating to fixed-point (using integer)
#define Fx2Fl(x)  (((float)x)/(1 << F))


fx_t sat_check(fx_t x)
{
    if (x > F_MAX)
        return F_MAX;
    else if (x < F_MIN)
        return F_MIN;
    else
        return x;
}

int main(void){

    float R = 0.5;
    double snr = 3.0;
    double sigma;

    double max = 0;


    const int N = 100; //length

    float sample[N] = { 0 };
    fx_t f_sample[N] = { 0 };  //fixed point
    fx_t q_sample[N] = { 0 };  // quantized

    int q_bit = 4; // the number of bit for quatization
    int q_max = pow(2, q_bit) - 1;

    //fixed-point 에서의 maximum 값, overflow 처리 필요.

    cout << "Max. Fixed-Point :" << F_MAX << endl;
    cout << "Min. Fixed-Point :" << F_MIN << endl;
    //getchar();
    sigma = pow(10, -0.05*snr) / sqrt(2 * R);

    //cout << "Orig.\t\tFixed\t\tFx_int\t\tFx_fl" << endl;
    for (int i = 0; i < N; i++)
    {
        float temp;
        sample[i] = 1 + gaussian(0, sigma);

        f_sample[i] = Fl2Fx(sample[i]);
        //temp = (int)(sample[i] * (1 << F));
        //temp = (float)temp / (1 << F);
        //cout << setw(8) << sample[i] << "\t" << temp << "\t\t";

        //printf("%d\t\t", f_sample[i]);
        //printf("%f\n", Fx2Fl(f_sample[i]));
        //putchar('\n');
    }
    cout << "F_samp\t\tOrig.\t\tSat" << endl;
    fx_t temp = 0;
    for (int i = 0; i < N; i++)
    {

        temp += f_sample[i];
        cout << f_sample[i]<<"\t\t"<<temp << "\t\t" << sat_check(temp) << endl;
        getchar();
    }
    cout << endl;
    cout << "maximum:  " << max << endl;

    return 0;
}
fixed_point_proto.md

fixed-point in C

original material: artist-embedded.org/EmbeddedControl Slides
reference 1: fixedpt.html
reference 2: Q_format

1. Fixed-point Representation

• x: real number
• X: fixed-point number
• N: wordlength
• m: integer (excluding sign bit)
• f: number of fraction bit
• “Q-format” : Qm.n

0/1101/011
sign bit/4bit integer/3bit fraction

2. Conversion to and from fixed-point

• real to fixed

  • Multiply the floating point number by 2^f
  • Round to the nearest integer
    $$X = round(x \dot{} 2^f)$$

• fixed to real

  $$x = X \dot{} 2^{-f}$$

example) 13.4 to Q4.3 format

3. Range of fixed-point representation

• negative number: 2’s complement

  • N=8, 2^(-8) ~ 2^(8-1)

    binary representation	decimal
    00000000	0
    00000001	1
    00000010	2
    …	…
    01111111	127
    10000000	-128
    10000001	-127
    …	…
    11111111	-1

• range of Qm.f [ref]

  $$[-2^m , 2^m - 2^f]$$

4. Arithmetic operations of fixed-point

• Satuation check

int16_t sat16(int32_t x)
{
    if (x > 0x7FFF) return 0x7FFF;
    else if (x < 0x8000) return 0x8000;
    else return (int16_t)x;
}

• Addition

int16_t q_add_sat(int16_t a, int16_t b)
{
    int16_t result;
    int32_t tmp;

    tmp = (int32_t)a + (int32_t)b;
    if (tmp > 0x7FFF)
        tmp = 0x7FFF;
    if (tmp < -1 * 0x8000)
        tmp = -1 * 0x8000;
    result = (int16_t)tmp;

    return result;
}

• Subtraction

int16_t q_sub(int16_t a, int16_t b)
{
    int16_t result;
    result = a - b;
    return result;
}

• Multiplication

// precomputed value:
#define K   (1 << (f - 1))    // f: fraction of fixed-point

int16_t q_mul(int16_t a, int16_t b)
{
    int16_t result;
    int32_t temp;

    temp = (int32_t)a * (int32_t)b; // result type is operand's type
    // Rounding; mid values are rounded up
    temp += K;
    // Correct by dividing by base and saturate result
    result = sat16(temp >> Q);

    return result;
}

• Division

int16_t q_div(int16_t a, int16_t b)
{
    int16_t result;
    int32_t temp;

    // pre-multiply by the base (Upscale to Q16 so that the result will be in Q8 format)
    temp = (int32_t)a << Q;
    // Rounding: mid values are rounded up (down for negative values).
    if ((temp >= 0 && b >= 0) || (temp < 0 && b < 0))
        temp += b / 2;
    else
        temp -= b / 2;
    result = (int16_t)(temp / b);

    return result;
}