Trigonometric Approximations

[curious_mind]

I was referring to link http://www.ganssle.com/item/approximations-for-trig-c-code.htm for referring to code for sine/cosine approx.

In the article the domain for evaluating cosine/sine is 0 to pi/2. Can we change the domain from 0 to 0.125 instead for quarter cycle? Purpose is to make the domain from 0 to 2*pi to 0 to 1 so that it eases out fixed point conversion. How can we recalculate the coefficients if the domain is changed?

 

[KlausST]

Hi,

0 to Pi/2 (angle in radians)
Means the same as
0° to 90° (angle in degree)
and is the angular range for table values / calculations.

Using this as "source" ... it's simple to get all other sine and cosine values just by ADD / SUBTRACT. Easy math even for small microcontrollers.

Referring to 0.125 makes no sense at all. I guess you misunderstood the idea behind the range.

Klaus

 

[curious_mind]

I implemented the entire equation for the range 0 to 2pi in simulink. It works. Next wanted to convert this implementation to fixed point. The issue is that , the coefficient c3=.03679 gives the result (1,32,34). I guessed that if I could change the input range from 0 to 2 pi to 0 to 1, then all the numbers will fit within 32 bit fixed point implementation. Pl advice.

 

[FvM]

Still unclear what you want to achieve. The original approximation is designed to cover the full sine argument range (full circle). That's not possible with the range reduction you are considering. You may use a different range if your application uses only a subrange of sine function.

I neither understand the problem about fixed point coefficients. Any numerical method involves finite number resolution and rounding errors. The usual way is to calculate the error bounds achieved with a specific number representation and to decide if they fulfill the requirements. MATLAB uses double precision by deafult. To model truncation and rounding effects, you can either refer to fixed point package or implement the effects in your code manually.

 

[curious_mind]

Ok, let me explain. The range 0 to 2*pi is numerically equals to 0 to 6.283. While converting to fixed point the coefficients which are smaller in number (such as c3 specified earlier) do not fit into the range of 32 bit fixed point. Subsequently I hooked up the sine/cosine model from simulink whose input ranges are 0 to 1 and outputs are -1 to +1. When I tried converting this model to fixed point, it does work. Then I guessed that I need to do the same thing with the chebyshev model as well. Hence 0.125 corresponds to pi/2 and 1 denotes 2*pi. Internally coefficients have to be re-evaluated to respond to the said input range to generate the respective outputs. I hope I am clear.

 

[FvM]

Hence 0.125 corresponds to pi/2 and 1 denotes 2*pi.

In my math, a quarter circle corresponds to 0.25 rather than 0.125 if the full circle is 1.
Rescaling angle argument to unity is often done in digital signal processing. Coefficients are rescaled respectively but this doesn't change much to relative accuracy requirements. Main adavantage is that periodicity beyond 0 to 2pi is automatically achieved by integer overflow.

I still don't understand about your specific problem with c3. None of the coefficients is represented exactly, neither in single, double or fixed point.

--- Updated Mar 23, 2023 ---

Obviously, the coefficients cn have to be scaled by (2pi)^(2*(n-1))

 

[curious_mind]

Oops that was a mistake. Yes it is 0.25 and not 0.125.
I applied the exact equation for 3.2 digit accuracy from the link http://www.ganssle.com/approx.htm
for cos(x) where x is 0-->pi/2, cos(x)=c1+x*x(c2+c3*x*x), here c1= 0.99940307, c2=-0.495558, c3=.03679168

results from matlab
c1- -->1,32,31
c2 ---> 1,32,32
c3----> 1,32,34 ( when I simulate with this number, output is highly inconsistent)
multipliers---> 1,32,31
input--> 1,32,16 ( 0 to 1.57)

Now please suggest the method by which I can bring c3 within the range.

 

[Aussie Susan]

(Perhaps it is my age or my culture but I don't understand the nomenclature used in 1,32,34 etc.)

Part of the underlying issue may be an understanding of what fixed point representations really are.
How you convert between the floating point and fixed point representations is completely up to you. It might be guided/constrained by the capabilities of the processor you are using (8-bit, 16-bit 'words' etc and whether it easily handles operations that are larger than its basic word size).
For example, you could represent your c3 value, 0.03679168, as 368 by multiplying by 10,000 and rounding. If you used that scheme then you need to scale all other numbers by a similar amount and/or take the different scaling factors into account at each step in the operation.
I'm guessing that you are multiplying the numbers by 360000 which makes c3 scale to 13,245 (which is not quite the number I think you are representing above - hence my first statement). (In my culture, the ',' character separates each group of 3 digits and '.' is the 'decimal point.)

Also, if you are finding that c3 cannot be represented as a 32-bit number, then how are you handling c1 and c3 which are numerically greater than c3?

There are other alternative ways of solving you underlying issue (there nearly always are). For example a lookup table which will trade off more data space and much simpler and faster 'calculation. The number of entries is determined by the precision you want in the domain (which for angles is very application dependent but generally 1 degree is more than enough so there are 90 entries. (The precision you need for the range will drive the size of the values in each table entry.)

Susan

 

[curious_mind]

I have upscaled the coefficients by 16 times, this brings all the coefficients in 16 or 32 fixed point format. However the angle output is also scaled by 16 and I have to bring it down for correct output readings. unfortunately 1/16 does not fit into the 16 or 32 fixed point format. Please offer your suggestions

 

[Aussie Susan]

I am not understanding this so I'll leave it to others to help.
For example, scaling up by 16 means shifting the binary value to the left 4 bits. Similarly scaling by 1/16 means a right shift of the binary by 4 bits. I really don't understand what you are saying that 1/16 does not fit a 16 or 32 bit fixed point number.
Susan.

 

[curious_mind]

I upscaled the float value and not the integer value. For example , coefficient c3=.036789. All I did was to multiply this by 16. I did this for all the 3 coefficients. The result is the angle (float) which is 16 times more than the intended. To bring it down, I multiplied the result with 1/16. Now when I started converting the float model, 1/16 value did not fit in dynamic range of Q15 or Q31 format.

 

[FvM]

Makes no sense to me. You are doing a polynomial approximation, scaling the angle results in scaling the coefficients with different power, see post#6.

 

[curious_mind]

FvM
Where does this factor as stated below come from? And what is "n"?. I think you may also suggest good reference guide on fixed point conversion

Obviously, the coefficients cn have to be scaled by (2pi)^(2*(n-1))

 

[FvM]

n is coefficient index, respectively related power of x . I obtained the factors by simply applying 2pi scaling to angle, as claimed in post #1 and post #3.