原理参考论文:Floating-Point_Inverse_Square_Root_Algorithm_Based_on_Taylor-Series_Expansion
主要代码如下:
clc
close all
clear
%%
a_init = [1.125 1.375 1.625 1.875]; %% 区间门限
a = a_init(1); %%
x = 1.12; % 1~2
y1 = 1/sqrt(x);
%% 基于泰勒展开式求开根号倒数
y2 = 1/sqrt(a)*(1-(x-a)/(2*a)+3*(x-a)^2/(8*a^2)-5*(x-a)^3/(16*a^3)+35*(x-a)^4/(128*a^4));
%% 误差对比
sub_out = y1-y2
%% 生成FPGA 数组系数
N = 20; %% 小数部分量化位数20位,整数部分量化2位
a_init_bin = zeros(4,N);
out = zeros(5,N);
for kk=1:4
a = a_init(kk);
A(1) = 1/sqrt(a);
A(2) = A(1)*1/(2*a);
A(3) = A(1)*3/(8*a^2);
A(4) = A(1)* 5/(16*a^3);
A(5) = A(1)* 35/(128*a^4);
a_init_bin(kk,:) = dectobin(a-1,N); %% 小数转定点数
B = 1/sqrt(2) * A;
for k1=1:5
% out(k1,:) = dectobin(A(k1),N); %% even number
out(k1,:) = dectobin(B(k1),N); %% odd number
end
switch kk
case 1
init0 = out ;
case 2
init1 = out ;
case 3
init2 = out ;
case 4
init3 = out ;
end
end
%% 测试部分:matlab浮点计算结果和FPGA定点计算结果对比
xin = [0.625 1.12 1.412 1.566 1.851 3.851 7.851]; % 选取7个值,可代表四个区间,以及奇偶
len = length(xin);
for hh=1:len
x = xin(hh);
x_out = 1/sqrt(x);
sqrt_matlab_out(hh) = x_out;
end
%% FPGA定点计算结果: compare fpga with matlab
load sqrt_out.txt
fpga_out = sqrt_out.';
for hh=1:len
fpga_bin = DEC_BIN(fpga_out(hh),N+3); %% 十进制转二进制数
bit_decimal = fpga_bin ;
bit1 = 3; %整数位宽
bit2 = N; %小数位宽
decimal = 0;
integer = 0;
%% 二进制数转整数
for j = 1 : bit1
integer = integer + bit_decimal(j) * (2^(bit1-j));
end
%% 二进制数转小数
for j = 1 : bit2
decimal = decimal + bit_decimal(j+bit1) * (2^(-j));
end
fpga_out_dec(hh) = integer + decimal;
end
fpga_sub_matlab = fpga_out_dec - sqrt_matlab_out;
figure
plot(fpga_sub_matlab);
title('fpga和matlab计算误码')
xlabel('point')
ylabel('Approximate errors')
从上图可以看出,matlab计算误码和FPGA计算误差10-6量级,如果位宽再扩大,以及划分更多区间,误差会进一步减低。
如需核心代码可私信。