# Orthogonal signal generator using integer arithmetic

I have a problem with implementing an orthogonal signal generator (OSG) algorithm on a microcontroller using integer arithmetic. I use this algorithm for a single-phase phase-locked loop (PLL) algorithm, for which I need an orthogonal component of a grid voltage.

The OSG algorithm is defined as follows:

$$\frac{d}{dt}v_x = \hat{\omega} \cdot \bigl((v_g-v_x)-v_y\bigr)$$ $$\frac{d}{dt}v_y = \hat{\omega} \cdot v_x$$

where $v_g$ is the measured grid voltage, $\hat{\omega}$ is the estimated grid frequency, and $v_x$ and $v_y$ are estimated components, with $v_x$ being equal to $v_g$ for ideal estimation. For this purpose, let us assume that the grid frequency is known.

The numerical integrator is implemented as follows:

$$y_k = \frac{T_s}{12} \bigr( 23u_{k-1} - 16u_{k-2} + 5u_{k-3} \bigl) + y_{k-1}$$

where $T_s=50~\mu\text{s}$ is the sample time.

Now, this algorithm works fine in floating point implementation, but is poor in integer arithmetic implementation. Here I give both implementations:

Floating point implementation

Code declaration.

float w = (2*PI)*50;
float Ts = 50e-6;
float i1u1, i1u2, i1u3, i1y1;
float i2u1, i2u2, i2u3, i2y1;

float NumInt3rd(float u1, float u2, float u3, float y1) {
return (Ts/12)*(23*u1-16*u2+5*u3)+y1;
}


Main function.

float vg = floor(Input(0));

float vg_x = NumInt3rd(i1u1,i1u2,i1u3,i1y1);
float vg_y = NumInt3rd(i2u1,i2u2,i2u3,i2y1);

i1u3 = i1u2;
i1u2 = i1u1;
i1u1 = ((vg-vg_x)-vg_y)*w;
i1y1 = vg_x;

i2u3 = i2u2;
i2u2 = i2u1;
i2u1 = vg_x*w;
i2y1 = vg_y;


The Input(0) is a macro to get an input signal (sine wave with an amplitude of $2048$).

Integer arithmetic implementation

Code declaration.

int w = 643398L; // (2*PI)*50*2048
int i1u1, i1u2, i1u3, i1y1;
int i2u1, i2u2, i2u3, i2y1;

int NumInt3rd(int u1, int u2, int u3, int y1) {
int iu = 23*u1-16*u2+5*u3;
int iy = 240000L*y1;
return (iu+iy)/240000L;
}


Main function.

int vg = (int) Input(0);

int vg_x = NumInt3rd(i1u1,i1u2,i1u3,i1y1);
int vg_y = NumInt3rd(i2u1,i2u2,i2u3,i2y1);

i1u3 = i1u2;
i1u2 = i1u1;
i1u1 = ((vg-vg_x)-vg_y)*w/2048;
i1y1 = vg_x;

i2u3 = i2u2;
i2u2 = i2u1;
i2u1 = vg_x*w/2048;
i2y1 = vg_y;


Note that I've checked for possible overflows, it never occurs. Also, the interesting thing is that the same algorithm works fine for $T_s=300~\mu\text{s}$.

I'm not that experienced with integer arithmetic implementations. Can you please give me an advice how to possibly fix this? Thanks!

Here is the estimation error for both implementations (on y-axis: percentage of the estimation error). The estimation error in case of integer artihmetic implementation is around $\pm 5\%$. • You have a scaling factor of 2048, corresponding to 22 bit resolution, when the signal amplitude is also 2048 (if I see correctly). Now, you divide by 240000, which is 18bit wide, hence you have only 4 bit real accuracy in the integration case. I believe this might not be enough. Try increasing the resolution to see what happens. – Maximilian Matthé Jan 12 '17 at 8:24
• Dear Maximilian, thank you for your post. I managed to solve the problem - it was because of the rounding errors after division, which is more pronounced for smaller sample times. For example, -4/3 rounds to -1, while -5/3 also rounds to -1. Because of this, the error is constantly accumulated. I'll post an answer with a more detailed explanation. – Marko Gulin Jan 12 '17 at 13:58

I managed to find an answer to my problem.

The problem is with rounding in integer division. For example, -4/3 will be rounded to -1, just as -5/3. Because of this, the integration error is constantly accumulated. Instead of explaining, here I give a code how to fix this.

Code declaration.

// Global variables
int w = 643398L; // (2*PI)*50*2048
int i1u1=0, i1u2=0, i1u3=0;
int i2u1=0, i2u2=0, i2u3=0;
short i1y1=0, i2y1=0;

// Numerical integrator implementation
short NumInt3rd(int u1, int u2, int u3, short y1) {
int iu = 23*u1-16*u2+5*u3;
int iy = 240000L*y1;
short y0 = (((iu+iy)>>12)*2237+65536)>>17;
return y0;
}


Main function.

// Get voltage measurements (-2048 to +2048)
short vg = (short) Input(0);

// Numerical integrators, 3rd order
short vg_x = NumInt3rd(i1u1,i1u2,i1u3,i1y1);
short vg_y = NumInt3rd(i2u1,i2u2,i2u3,i2y1);

// Downsample frequency to prevent overflow
short wb = w>>5;

// Update integrator #1 states
i1u3 = i1u2;
i1u2 = i1u1;
i1u1 = ((int)((vg-vg_x)-vg_y)*wb)>>6;
i1y1 = vg_x;

// Update integrator #2 states
i2u3 = i2u2;
i2u2 = i2u1;
i2u1 = ((int)vg_x*wb)>>6;
i2y1 = vg_y;


It should be noted that instead of using integer division, which is very expensive in terms of required number of instruction cycles, I rather use bit shift operation. For example, 1/240000 can be well approximated as 2237/2^29 - the approximation accuracy is 0.00169277%.

And here is the sine wave estimation error. As you can see, the results are much better now. 