CORDIC

From Hamsterworks Wiki!

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!

Cordic sim.png Cordic grab.png

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.

Contents

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;

Personal tools