CORDIC help page gives inaccurate results
2 views (last 30 days)
Show older comments
I am trying to use the code provided at this help page for computing square roots using CORDIC:
I took the provided code, added the input initialization and output normalization (as described in the article) to get this function:
function s = cordic_sqrt(v, n)
% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
% function s = cordic_sqrt(v, n)
% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
% Calculates the (approximate) square root of v over n CORDIC iterations.
%
% Taken from the Matlab help pages:
% mathworks.com/help/fixedpoint/examples/compute-square-root-using-cordic.html
% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
% Set initial values
x = v + 0.25;
y = v - 0.25;
k = 4; % Used for the repeated (3*k + 1) iteration steps
for idx = 1:n
xtmp = x * 2^(-idx);
ytmp = y * 2^(-idx);
if y < 0
x = x + ytmp;
y = y + xtmp;
else
x = x - ytmp;
y = y - xtmp;
end
if idx==k
xtmp = x * 2^(-idx);
ytmp = y * 2^(-idx);
if y < 0
x = x + ytmp;
y = y + xtmp;
else
x = x - ytmp;
y = y - xtmp;
end
k = 3*k + 1;
end
end
% Calculate normalization
An = prod(sqrt(1 - 2.^(-2*(1:n))));
% Normalize
s = x / An;
However, this gives very inaccurate results. For example:
cordic_sqrt(3/8, 30)
ans =
0.611175220934138
(whereas MATLAB's sqrt function gives sqrt(3/8) = 0.612372435695794). I notice that if I comment out this part of the code:
% if idx==k
% xtmp = x * 2^(-idx);
% ytmp = y * 2^(-idx);
% if y < 0
% x = x + ytmp;
% y = y + xtmp;
% else
% x = x - ytmp;
% y = y - xtmp;
% end
% k = 3*k + 1;
% end
then I get accurate results (identical to MATLAB's sqrt function). So what is going on here?
I feel suspicious about the quality of the code provided in the article. I don't think the definition of the normalization term is exactly correct (since it doesn't take into account these extra "if idx==k" iterations). However, this won't make any difference in my example because my number of iterations is so high (30).
0 Comments
Answers (0)
See Also
Categories
Find more on Point Cloud Processing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!