How CORDIC sin/cos in HDL works.

[hamster_nz]

Somebody asked me how this CORDIC code works to calculate SIN() and COS() from a phase angle, using only addition, subtraction, and division by a power of two. So I've made the visual analogy that I've been thinking about of a while...

Have a look at the first picture.

How it's constructed: It starts with a right angled triangle with two sides of equal length, so a 45 degree angle.

On its hypotenuse is another right angled triangle, and its short side is 1/2th the length of the first triangle's hypotenuse. It has an angle of 26.56 degrees.

On that triangle's hypotenuse is another right angled triangle, and it's short side 1/4th of the previous triangle's hypotenuse. It has an angle of 13.50 degrees

On that triangles's hypotenuse is another right angled triangle, and it's short side 1/8th of the previous trangle's hypotenuse. It has an angle of 7.10 degrees

On that triangles's hypotenuse is another right angled triangle, and it's short side 1/16th of the previous trangle's hypotenuse. It has an angle of 3.57 degrees

So laid out like in the first image, it has a total angle of 45.0 + 26.56 + 13.59 + 7.13 + 3.57 = 95.85 degrees.

So far pretty boring, but the key to CORDIC is that you can fold the paper!

Because of this you can get 32 different angles between -95.85 degrees and +95.85 degrees, with the angles of (+/- 45.0 +/- 26.56 +/- 13.59 +/- 7.13 +/- 3.57).

The second photo has an example of

45.0 - 26.56 + 13.59 + 7.13 + 3.57 = 42.73 degrees.

That is what this bit of the code is doing. It is deciding which way to add the triangle (clockwise or counter-clockwise) to get closer to the desired angle,  then calculates where the corner of this new triangle will be.

Code: [Select]

               if remainder(i) < 0 then
                   remainder(i+1) <= remainder(i) + angles(i);
                   x(i+1)   <= x(i) + y(i)(19 downto i+1);
                   y(i+1)   <= y(i) - x(i)(19 downto i+1);
               else
                   remainder(i+1) <= remainder(i) - angles(i);
                   x(i+1)   <= x(i) - y(i)(19 downto i+1);
                   y(i+1)   <= y(i) + x(i)(19 downto i+1);
               end if;


If the error between the current angle and the angle we want is < 0, then it adds (y(i), -x(i)) scaled by 1/2^n onto (x(i),y(i)) to get (x(i+1),y(i+1))

If the error between the current angle and the angle we want is >= 0, then it adds (-y(i), x(i)) scaled by 1/2^n onto (x(i),y(i))  to get (x(i+1),y(i+1))

Both paths also updates the error angle using a lookup table 'angles', ready for the next, narrower triangle to be added.

As the process repeats, the error value gets "driven to zero", and the last x() and y() values are a scaled cosine and sine of the 'angle' phase angle input.

My code is just a little bit different, as it first decides what quadrant the angle is in, then works from the diagonal of that quadrant to the desired angle, so starts with the 26.56 degree triangle. It also start of at (+/-30000,+/-30000) so the output range for the sin() and cos() values will be close to +/-32767.

Full code for the module is here:

Code: [Select]

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

entity sin_cos is
    Port ( clk     : in  STD_LOGIC;
           de_in   : in  std_logic;
           angle   : in  STD_LOGIC_VECTOR (19 downto 0);
           de_out  : out std_logic := '0';
           sine    : out STD_LOGIC_VECTOR (19 downto 0)  := (others=>'0');
           cosine  : out STD_LOGIC_VECTOR (19 downto 0)  := (others=>'0'));
end sin_cos;

architecture Behavioral of sin_cos is
    type a_int   is array(0 to 15) of integer;
    constant angles : a_int := ( 77376, 40884, 20753, 10417,
                                  5213,  2607,  1304,   652,
                                   326,   163,    81,    41,
                                    20,    10,     5,     3);
    type a_value is array(0 to angles'high+1) of signed(19 downto 0);
    signal de_sr : std_logic_vector(angles'high+1 downto 0) := (others => '0');
    signal remainder : a_value   := (others => (others => '0'));
    signal x         : a_value  := (others => (others => '0'));
    signal y         : a_value := (others => (others => '0'));
begin

    de_out <= de_sr(de_sr'high);
    sine   <= std_logic_vector(x(x'high));
    cosine <= std_logic_vector(y(y'high));

process(clk)
    begin
        if rising_edge(clk) then
            for i in 0 to angles'high loop
               de_sr(i+1) <= de_sr(i);
               if remainder(i) < 0 then
                   remainder(i+1) <= remainder(i) + angles(i);
                   x(i+1)   <= x(i) + y(i)(19 downto i+1);
                   y(i+1)   <= y(i) - x(i)(19 downto i+1);
               else
                   remainder(i+1) <= remainder(i) - angles(i);
                   x(i+1)   <= x(i) - y(i)(19 downto i+1);
                   y(i+1)   <= y(i) + x(i)(19 downto i+1);
               end if;
            end loop;
           
            de_sr(0) <= de_in;
            if de_in = '1' then
                case angle(angle'high downto angle'high-1) is
                    when "00"   => x(0) <=  to_signed( 300000,20);  y(0) <=  to_signed( 300000,20); remainder(0) <= signed(angle) - x"20000";
                    when "01"   => x(0) <=  to_signed(-300000,20);  y(0) <=  to_signed( 300000,20); remainder(0) <= signed(angle) - x"60000";
                    when "10"   => x(0) <=  to_signed(-300000,20);  y(0) <=  to_signed(-300000,20); remainder(0) <= signed(angle) - x"A0000";
                    when others => x(0) <=  to_signed( 300000,20);  y(0) <=  to_signed(-300000,20); remainder(0) <= signed(angle) - x"E0000";
                end case;
            end if;
           
        end if;
    end process;
end Behavioral;


[BrianHG]

A simpler VHDL fixed radius square version (equal X&Y radius) of my SystemVerilog ellipse generator I wrote for NockieBoy's 8bit GPU & my DDR3 graphics demo.

[BrianHG]

Since 'rstofer' requested it, here is my code:  https://github.com/BrianHGinc/BrianHG-DDR3-Controller/blob/main/BrianHG_DDR3_GFX_source/ellipse_generator.sv

Lines 326 through 341 step a single value at a time, calculating 0-45 degrees.

The ry2 & rx2 hold the computed x&y radius squared where 'p' is the error offset.
'x' and 'y' coordinates just incs and decs 1 depending on the polarity of the 'p' error accumulated every iteration.

Since my algorithm was designed to draw pixels, it is called 2 times (once with inv enabled & disabled drawing 0-45 deg, then 90-45 deg where it is actually plotting the 0-45 backwards with X&Y coords backwards.) x 4 times to build the 4 quadrants of a complete ellipse.

Note that it has been a year since I looked at this code.


Obviously when drawing a true circle, all quadrants within the 45 degrees can all simultaneously be derived from that first 45 degree.

hamster_nz solution is so much easier as he has a fixed squares & radius.  It looks a more like Bresenham's circle algorithm which my ellipse engine is based on.

[BrianHG]

Functional Ellipse generator testbench: https://github.com/BrianHGinc/SystemVerilog-TestBench-BPM-picture-generator

[rstofer]

Since 'rstofer' requested it, here is my code:  https://github.com/BrianHGinc/BrianHG-DDR3-Controller/blob/main/BrianHG_DDR3_GFX_source/ellipse_generator.sv

Thanks!  This is my first exposure to something significant written in System Verilog.  I have a lot to learn!

[BrianHG]

Since 'rstofer' requested it, here is my code:  https://github.com/BrianHGinc/BrianHG-DDR3-Controller/blob/main/BrianHG_DDR3_GFX_source/ellipse_generator.sv

Thanks!  This is my first exposure to something significant written in System Verilog.  I have a lot to learn!

Start with the BMP generator version testbench as the DDR3 is a monster of a complete project.

The BMP generator will show you how to make a testbench to test a SystemVerilog module where my ELLI module is being demonstrated.

The DDR3 project also has test-benches for each section except the ellipse generator.