CORDIC help page gives inaccurate results

1 次查看(过去 30 天)
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 个)

类别

Help CenterFile Exchange 中查找有关 Point Cloud Processing 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by