vhdl/lib/radio/cordic/eval_cordic.m

function eval_cordic(Nmax, doplot)
N = 100;

arctan_tbl = (atan(2.^-(0:N-1))/pi)'
%gain_tbl = (1./cumprod(sqrt(1+2.^-(2*(0:N-1)))))'

% Rotation mode
disp('Rotation:')
for i=1:N,
    phi(i) = 2*pi*(i-N/2)/N;
    [x(i),y(i)] = cordic_rot(0, 1, phi(i), Nmax, 0);
end;
close all;
figure
plot(1:N, x, 1:N, y);
legend('x','y')
grid;

% Vector mode
disp('Magnitude and phase:')
for i=1:N,
    [mag(i), phi(i)] = cordic_vec(x(i), y(i), Nmax, 0);
end;
figure
plot(1:N, mag, 1:N, phi);
legend('Magnitude','Phase')
grid;


% ------------------------------------------------------------
function [mag, phi] = cordic_cart2polar(x, y, Nmax, doplot)
    [mag, phi] = cordic_vec(x, y, Nmax, doplot);
  
% ------------------------------------------------------------
function [x, y] = cordic_polar2cart(mag, phi, Nmax, doplot)
    [mag, phi] = cordic_vec(x, y, Nmax, doplot);
  
% ------------------------------------------------------------
function [xout, yout] = cordic_rot(xin, yin, phi, N, doplot)

    % Initialize
    arctan_tbl = atan(2.^-(0:N-1));
    gain_tbl = 1./cumprod(sqrt(1+2.^-(2*(0:N-1))));

    coeffs = arctan_tbl';
    gains = gain_tbl';

    % Calculate initial rotation
    if (phi < (-pi/2))
        z = phi + pi;
        x = -xin;
        y = -yin;
    else
        if (phi > pi/2)
            z = phi - pi;
            x = -xin;
            y = -yin;
        else
            x = xin;
            y = yin;
            z = phi;
        end;
    end;

    % Loop
    for i=1:N,
        xi(i) = x;
        yi(i) = y;
        zi(i) = z;
    
        S = 2^(-(i-1));
        if (z < 0)
			t = x + y*S;
			y = y - x*S;
			z = z + arctan_tbl(i);
			x = t;
        else
			t = x - y*S;
			y = y + x*S;
			z = z - arctan_tbl(i);
			x = t;
        end;
        
    end;

    % Result
    k = gain_tbl(i);
    xout = k*x;
    yout = k*y;

    % Output
    if(doplot)
        close all;
        figure(1);
        plot(1:N, xi, 1:N, yi, 1:N, zi)
        legend('x','y','z');
        title('Convergence');
        grid;

        figure(2);
        plot(xi.*gain_tbl, yi.*gain_tbl, '-x')
        title('Polar plot');
        grid;
    end;

% ------------------------------------------------------------
function [mag, phi] = cordic_vec(xin, yin, N, doplot)

    % Initialize
    arctan_tbl = atan(2.^-(0:N-1));
    gain_tbl = 1./cumprod(sqrt(1+2.^-(2*(0:N-1))));

    coeffs = arctan_tbl';
    gains = gain_tbl';
    
    % Calculate initial rotation
    if (xin < 0)
        x = -xin;
        y = -yin;
        if (yin < 0)
            z = -pi;
        else
            z = pi;
        end;
    else
        x = xin;
        y = yin;
        z = 0;
    end;
  
    % Loop
    for i=1:N,
        xi(i) = x;
        yi(i) = y;
        zi(i) = z;
    
        S = 2^(-(i-1));
        if (y < 0)
			t = x - y*S;
			y = y + x*S;
			z = z - arctan_tbl(i);
			x = t;
        else
			t = x + y*S;
			y = y - x*S;
			z = z + arctan_tbl(i);
			x = t;
        end;
        
    end;

    % Result
    k = gain_tbl(i);
    mag = k*x;
    phi = z;
    
    % Output
    if(doplot)
        close all;
        figure(1);
        plot(1:N, xi, 1:N, yi, 1:N, zi)
        legend('x','y','z');
        title('Convergence');
        grid;

        figure(2);
        plot(xi.*gain_tbl, yi.*gain_tbl, '-x')
        title('Polar plot');
        grid;
    end;