%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Creative Commons
% Attribution-Noncommercial 2.5 India
% You are free:
% to Share — to copy, distribute and transmit the work
% to Remix — to adapt the work
% Under the following conditions:
% Attribution. You must attribute the work in the manner 
% specified by the author or licensor (but not in any way 
% that suggests that they endorse you or your use of the work). 
% Noncommercial. You may not use this work for commercial purposes. 
% For any reuse or distribution, you must make clear to others the 
% license terms of this work. The best way to do this is with a 
% link to this web page.
% Any of the above conditions can be waived if you get permission 
% from the copyright holder.
% Nothing in this license impairs or restricts the author's moral rights.
% http://creativecommons.org/licenses/by-nc/2.5/in/

% Checked for proper operation with Octave Version 3.0.0
% Author	: Krishna
% Email		: krishna@dsplog.com
% Version	: 1.0
% Date		: 9 December 2007
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Using CORDIC for phase and magnitude computation
clear

% creating the reference angles
Y_k = ones(32,1) + j*2.^(-1*[0:1:31])'; % values of Y
alpha_k = angle(Y_k)*180/pi; % reference angles, used for finding the cumulative phase

X = 7+3*j; % complex number for which magnitude is to be computed

K = 16; % number of iterations
Z = X; % assigning the output as the input

% first rotation - for shifting to first/fourth quadrant
if imag(Z) <=0
% rotate by 90 degrees
Z = -1*imag(Z) + j*real(Z); % same as Z = Z*exp(j*pi/2);
thetaHat = 90;
else
% rotate by -90 degrees
Z = imag(Z) - j*real(Z); % same as Z = Z*exp(-j*pi/2);
thetaHat = -90;
end

% running for different values of k
for k = 0: K-1

s = -1*sign(imag(Z)); % difference
if s == 0 % to handle when sign(imag(Z)) = 0
s = 1;
end
Z = Z*[1 + s*j*2^(-1*k)];
thetaHat = thetaHat + s*alpha_k(k+1);

% dumping the variables sign_v(k+1) = s;
thetaHat_v(k+1) = thetaHat;
real_Z(k+1) = real(Z);
end

% comparing the absolute value obtained via cordic
Z_abs_ideal = abs(X);
Z_abs_cordic = Z/1.64676025786545;
err_abs = Z_abs_ideal - real(Z_abs_cordic)

Z_angle_ideal = angle(X)*180/pi;
Z_angle_cordic = -1*thetaHat;
err_angle = Z_angle_ideal - Z_angle_cordic