CORDIC help page gives inaccurate results
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (0)
Siehe auch
Kategorien
Mehr zu Point Cloud Processing finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!