# CORDIC

Jump to: navigation, search

This FPGA Project was finished in September 2018

I wanted to test an 24-bit audio DAC convertor board I've made up, and needed a high quality sine generator.

So I implemented CORDIC for the first time. Rather than generating one n-bit value every n+2 cycles (which is the usual way to minimize resource usage) I've decided to generate one value every cycle, as might be done in a high performance SDR environment. Why calculate 50,000 samples per second, when you can make 100,000,000!

I connected it up to a I2S DAC and the results look good, based on my old 8-bit scope - in practice that means nothing more than it is not completely broken!

The CORDIC algorithm allows you to generate high precision SIN() and COS() values using a successive approximation.

It uses a table that is the contains the series of ATAN(0.5), ATAN(0.25), ATAN(0.125), ATAN(0.0625)..., and needs a one table entry for each bit of resolution in the result.

Interestingly, both the input and output can be scaled to any unit you want - you don't need to use degrees, or radians. This makes it a good match for a binary phase accumulator.

The width of this lookup table needs to be a few bits wider than the input angle, otherwise the smallest adjustments get lost.

## Tricks and traps

Division by two using binary shifts doesn't work as expected for negative numbers - it is actually gives (x-1)/2. This gives some error if you don't make allowances for it.

For high precision rounding also becomes an issue, and getting that last bit of accuracy can be very expensive.

As the approximation gets closer to zero, you need less and less bits in the 'z' (angle) value. Most synthesis tools will not recognize this and use too many bits.

## A C hack of the algorithm

This code calculates sin() and cos() using CORDIC, and compares the results with that of the floating point library. The input angle is scaled so 2^24 is the full circle.

```#include <stdio.h>
#include <math.h>
#include <stdint.h>

#define PI             (3.14159265358979323846)
#define FULL_CIRCLE    (1<<24)
#define EIGHTH_CIRCLE  (FULL_CIRCLE/8)
#define QUARTER_CIRCLE (FULL_CIRCLE/4)
#define HALF_CIRCLE    (FULL_CIRCLE/2)
#define CORDIC_REPS    (30)
#define BITS           (30)

int64_t initial;
double angles[CORDIC_REPS];
/****************************************************************
* Calculate the values required for CORDIC sin()/cos() function
***************************************************************/
void setup(void) {
int i;
double scale = pow(0.5,0.5);

for(i = 0; i < CORDIC_REPS; i++ ) {
double angle = atan(1.0/pow(2,i+1));
angles[i]  = FULL_CIRCLE * angle / (2*PI);
scale     *= cos(angle);
printf("angle[%i] = %13.9f\n",i, angles[i]);
}
initial = (int64_t)(scale*pow(2.0,BITS));
}

/***************************************************************
* Cordic routine to calculate Sine and Cosine for angles
* from 0.0 to 2*PI
**************************************************************/
void cordic_sine_cosine(double z, int64_t *s, int64_t *c) {
int i, flip_sine_sign = 0, flip_cos_sign = 0;
int64_t x = initial,y = initial;
if(z > HALF_CIRCLE) {
z = FULL_CIRCLE-z;
flip_sine_sign = 1;
}

if(z > QUARTER_CIRCLE) {
z = HALF_CIRCLE-z;
flip_cos_sign = 1;
}
z -= EIGHTH_CIRCLE;
for(i = 0; i < CORDIC_REPS; i++ ) {
int64_t tx = (x + ((int64_t)1<<i)) >> (i+1);
int64_t ty = (y + ((int64_t)1<<i)) >> (i+1);
x -= (z > 0 ?        ty :        -ty);
y += (z > 0 ?        tx :        -tx);
z -= (z > 0 ? angles[i] : -angles[i]);
}
*c = flip_cos_sign  ? -x : x;
*s = flip_sine_sign ? -y : y;
}

/**************************************************************/
int main(int argc, char *argv[]) {
double a = 0.8, max = 0.0;
double total_e = 0.0;
int count = 0;
setup();

for(a = 0; a < FULL_CIRCLE; a++) {
int64_t s, c;
double es,ec;

cordic_sine_cosine(a, &s, &c);
es = s-(int64_t)(sin(a*(2*PI/FULL_CIRCLE))*(pow(2,BITS)));
ec = c-(int64_t)(cos(a*(2*PI/FULL_CIRCLE))*(pow(2,BITS)));

if(es > 0) total_e += es;
else       total_e -= es;

if(ec > 0) total_e += ec;
else       total_e -= ec;

if(max < es)  max = es;
if(max < -es) max = -es;
if(max < ec)  max = ec;
if(max < -ec) max = -ec;
count++;
}
printf("Error is %13.11f per calcuation out of +/-%f\n",total_e/count, pow(2,BITS));
printf("Max error is %13.11f\n",max);
return 0;
}

```

## Source

### cordic_sin_cos.vhd

```library IEEE;
use IEEE.STD_LOGIC_1164.ALL;
use IEEE.NUMERIC_STD.ALL;

entity cordic_sin_cos is
Port ( clk        : in  STD_LOGIC;
angle      : in  STD_LOGIC_VECTOR (23 downto 0);
signed_sin : out STD_LOGIC_VECTOR (23 downto 0) := (others => '0');
signed_cos : out STD_LOGIC_VECTOR (23 downto 0) := (others => '0'));
end cordic_sin_cos;

architecture Behavioral of cordic_sin_cos is
constant stages    : natural := 26;
constant calc_bits : natural := 28;
constant z_bits    : natural := 29;

-- Angles for the CORDIC algorithm in fixed decimal (eg. 0x4b90147 ~= atan(0.5))
type a_angles is array(0 to stages) of signed(z_bits-2 downto 0);
constant angles : a_angles := (
x"4b90147", x"27ece16", x"1444475", x"0a2c350", x"05175f8", x"028bd87",
x"0145f15", x"00a2f94", x"00517cb", x"0028be6", x"00145f3", x"000a2f9",
x"000517c", x"00028be", x"000145f", x"0000a2f", x"0000517", x"000028b",
x"0000145", x"00000a2", x"0000051", x"0000028", x"0000014", x"000000a",
x"0000005", x"0000002", x"0000001");

-- Adjusted to give full scale of +0x7FFFF & -0x7FFFF
constant initial      : signed(calc_bits-1 downto 0) := x"4dba76c";

-- Constants for quadrant calulations
constant quarter_turn : signed(z_bits-1 downto 0)    := "10000" & x"000000";
constant eigth_turn   : signed(z_bits-1 downto 0)    := "01000" & x"000000";
constant zero         : signed(signed_sin'high downto 0) := (others => '0');
signal quadrant_angle : signed(z_bits-1 downto 0)    := (others => '0');

-- Pipeline for holding the quadrant information
signal flip_cos : std_logic_vector(stages+1 downto 0);
signal flip_sin : std_logic_vector(stages+1 downto 0);

-- Pipeline for CORDIC X, Y & Z state variables positions
type a_x_or_y is array(0 to stages+1) of signed(calc_bits-1 downto 0);
type a_z      is array(0 to stages+1) of signed(z_bits-1 downto 0);
signal x       : a_x_or_y := (0 => initial, others => (others => '0'));
signal y       : a_x_or_y := (0 => initial, others => (others => '0'));
signal z       : a_z      := (others => (others => '0'));

begin

quadrant_angle(quadrant_angle'high-1 downto quadrant_angle'high-angle'high+1)  <= signed(angle(angle'high-2 downto 0));

process(clk)
begin
if rising_edge(clk) then
-- The Output stages
if flip_sin(stages+1) = '0' then
signed_sin <= STD_LOGIC_VECTOR (zero + y(stages+1)(calc_bits-1 downto calc_bits-signed_sin'length));
else
signed_sin <= STD_LOGIC_VECTOR (zero - y(stages+1)(calc_bits-1 downto calc_bits-signed_sin'length));
end if;

if flip_cos(stages+1) = '0' then
signed_cos <= STD_LOGIC_VECTOR (zero + x(stages+1)(calc_bits-1 downto calc_bits-signed_cos'length));
else
signed_cos <= STD_LOGIC_VECTOR (zero - x(stages+1)(calc_bits-1 downto calc_bits-signed_cos'length));
end if;

-- The actual CORDIC stages
for i in 0 to stages loop
if z(i)(z_bits-1) = '0' then
x(i+1) <= x(i) - resize(y(i)(calc_bits-1 downto i+1),calc_bits);
y(i+1) <= y(i) + resize(x(i)(calc_bits-1 downto i+1),calc_bits);
z(i+1) <= z(i) - angles(i);
else
x(i+1) <= x(i) + resize(y(i)(calc_bits-1 downto i+1),calc_bits);
y(i+1) <= y(i) - resize(x(i)(calc_bits-1 downto i+1),calc_bits);
z(i+1) <= z(i) + angles(i);
end if;
flip_sin(i+1) <= flip_sin(i);
flip_cos(i+1) <= flip_cos(i);
end loop;

-- Setup for the CORDIC algorithm, which only works on the first quadrant of a circle
case angle(angle'high downto angle'high-1) is
when "00" =>
flip_sin(0) <= '0';
flip_cos(0) <= '0';
z(0) <=                quadrant_angle - eigth_turn;
when "01" =>
flip_sin(0) <= '0';
flip_cos(0) <= '1';
z(0) <= quarter_turn - quadrant_angle - eigth_turn;
when "10" =>
flip_sin(0) <= '1';
flip_cos(0) <= '1';
z(0) <=                quadrant_angle - eigth_turn;
when others =>
flip_sin(0) <= '1';
flip_cos(0) <= '0';
z(0) <= quarter_turn - quadrant_angle - eigth_turn;
end case;
end if;
end process;
end Behavioral;

```