RANSAC算法的FPGA实现

目录

一.算法原理

1.参数变量

2.算法流程

<1>随机读取数据

<2>拟合模型

<3>检验模型

<4>循环遍历

二.设计思路

1.设计目标

2.模块设计

<1>随机数产生模块

<2>暂存寄存器

<3>计算模块

<4>模型评估模块

<5>主体控制与输出模块

<6>整合

三.具体实现

1.Python数据生成模块

2.Verilog RANSAC算法仿真结果

四.性能分析

1.资源分析

2.速度分析

3.选择不同随机数种子下的结果

4.不同预设直线参数下的拟合结果

五.思考与改进

1.提速瓶颈

2.增加并行程度

3.算法流程改进

六.附录

1.vivado工程源码部分

<1>主体控制部分

七.参考文献

八.设计改进

1.不可综合语句

2.修正程序语句

3.时钟模块

九.设计实现分析

1.静态时序分析

十.下载验证

---

参考：

伪随机序列：LFSR Sequence、m-Sequence和Gold Code/Sequence-CSDN博客RANSAC算法(附RANSAC直线拟合C++与Python版本)-CSDN博客

RANSAC算法——看完保证你理解-CSDN博客伪随机序列：LFSR Sequence、m-Sequence和Gold Code/Sequence-CSDN博客

一.算法原理

1.参数变量

输入：数据集data、拟合模型model

中间参数：一次拟合需要的数据量n；算法最大遍历次数k；计算匹配阈值t；最小匹配数据数d

输出：完成信号done、最匹配的模型model参数

模型：二元一次线性模型，二元多次非线性模型，任意函数...

2.算法流程

<1>随机读取数据

随机数生成模块生成n个随机地址，取出存储器对应地址数据。

<2>拟合模型

由读出的n个数据拟合出相应模型

<3>检验模型

检验存储器中所有数据对于该模型的拟合程度

<4>循环遍历

每次循环评价随机拟合的模型，当遍历k次后结束拟合，得到最佳模型

二.设计思路

1.设计目标

        工程能够完成一组坐标点集的线性拟合。输入坐标数据（x,y）的值为8bit 无符号整数数据类型，输出拟合直线参数，包括斜率k和y轴截距b，为8bit 有符号整数。

2.模块设计

<1>随机数产生模块

        利用线性反馈移位寄存器（Linear Feedback Shift Register，LFSR）产生随机数，对于反馈移位寄存器和异或门构成的电路，可以按下面的关系式生成序列[1]：

        对于m位的线性反馈移位寄存器，最多可以产生2^m-1个不同状态。如果一个序列发生器正好生产这 2^m-1个不同状态之后才重复此序列。随机数产生器的初始队列成为种子（seed），随后队列对种子进行移位运算生成新队列。

        工程中采用32位随机数产生模块，每次取前14位分别作为两随机点的地址（128个数据对应7位地址）。

<2>暂存寄存器

        通过rom读取文件内数据信息，存入相关数据。Rom深度为128，对应7位地址线，宽度为16位数据前8位为x，后8位为y。采用组合逻辑电路直接输出对应地址数据。

<3>计算模块

        通过所给数据拟合出模型。输入为32位数据，前16位和后16位分别是(x1,y1)和(x2,y2)。实现的组合电路结构如下：

K=(y2-y1)/(x2-x1)    b=y2-k*x2

<4>模型评估模块

        给定参数T_RANGE，表示允许残差平方RSS的最大值（RSS=(y-k*x-b)^2，表示期望值与实际值只差的平方），当RSS<T_RANGE时，表示该点在允许范围之内。Fit信号输出为1。实现的组合电路结构如下：

<5>主体控制与输出模块

        用于控制迭代次数K，以及产生随机数产生模块时钟和rom地址更新。每次迭代重新生成随机数以产生新模型。遍历rom检验后比较最大线内点数fit_num，如果大于最大值则更新模型参数

<6>整合

NAME	ADDR_WIDTH	DATA_WIDTH		INPUTDATA_WIDTH				PARA_WIDTH		T_WIDTH
WIDTH/BITS	7	16		32				16		16
DETAILS		x(h_8b)	y(l_8b)	x1(h_8b)	y1(l_8b)	x2(h_8b)	y2(l_8b)	k(h_8b)	b(l_8b)

三.具体实现

1.Python数据生成模块

        使用python编写代码生成十六进制数据集，包含线内点和噪声点，并给出python实现ransac和最小二乘法的拟合结果。如下：

        设定直线的参数：斜率k=2，截距b=25，最小二乘法和ransac拟合参数结果如下：

        其中红色线是最小二乘法的拟合结果，绿色线是ransac在迭代1000次，允许误差为5的拟合结果。

        将产生数据以16进制形式保存，并存储于data.mem文件中在vivado中读取。

2.Verilog RANSAC算法仿真结果

        工程拟合结果直线参数为：斜率k=2(0x02),b=24(0x18)，完全符合设定值。符合误差范围的点数最多为47个（0x2F），并在第450代（0x01C2）时随机产生的地址拟合出了最佳参数。此时的残差平方和RSS为2086（0x826）。

        由于拟合和评估部分皆为组合逻辑电路，算法完成的速度取决于迭代次数K和时钟频率CLK_FRE。仿真采用100MHZ时钟，K取1000次，每次迭代遍历128个点需要128个时钟周期，故算法每次执行需要：128*K/CLK_FRE秒。本次仿真算法在0.00128s（128*1k/1M）左右结束，done信号置高电平有效，故算法频率为781.25HZ。

四.性能分析

1.资源分析

        设计使用到的查找表、IO、缓存器等资源均较少，满足设计的资源限制。

2.速度分析

        下表给出不同迭代次数和时钟频率下的算法频率F=CLK_FRE/(128*K)：

F=CLK_FRE/(128*K)
CLK_FRE(/M)	K	frequency(/Hz)
100	1000	781.25
100	10000	78.125
100	500	1562.5
80	1000	625
15	1000	117.1875

3.选择不同随机数种子下的结果

        在迭代次数k=1000次下，不同随机数种子下最终结果出现的时机（其中k是代数，para是拟合参数）：

seed	0x12345678	0xFFFFFFFF	0x11111111	0xcd62f912	0x63e8c090
K	450	403	13	524	360
PARA	0x0218	0x0218	0x0218	0x0218	0x0218

4.不同预设直线参数下的拟合结果

        在seed=0x12345678下，不同预设直线参数下的拟合结果如下，其中k，b为设定值，ransac_k、ransac_b为拟合值，fit_num为检测为线内的点数（实际为128）：

k	2	3	4	2.5
ransac_k	2	3	4	2
b	25	10	77	16
ransac_b	24	11	51	41
fit_num	47	64	25	39

五.思考与改进

1.提速瓶颈

        在算法的实现流程中，计算和评估模块均采用组合逻辑电路，而 控制部分电路采用循环遍历的方式，依次选取rom中数据对模型进行检验，大大降低了算法的速度。

2.增加并行程度

        在检验阶段，设计一寄存器组，提前将rom中所有数据存储至寄存器组后并行输出，连接至同等数量的evaluate_kernel模块同时进行检验。实现由之前的循环遍历校验转化为同时并行检验：

3.算法流程改进

<1>工程中模型的拟合采用直接随机选取n个数据进行拟合的方法，进而确定线内点。可以改进为最后增加一个步骤，由已经确定的线内点重新使用最小二乘法拟合参数，减小拟合参数的残差平方和RSS。

<2>线内点的确定可改进为投票式，即由k个模型共同投票符合要求的片内点，选取票数最高的点作为真正的片内点[2]。

六.附录

1.vivado工程源码部分

<1>主体控制部分

`timescale 1ns / 1ps
//
// Engineer: switch_swq
// Create Date: 2024/05/20 10:48:12
//
module evaluate
#(
    parameter ADDR_WIDTH = 7,
    parameter DATA_WIDTH = 16,
    parameter INPUTDATA_WIDTH = 32,
    parameter PARA_WIDTH = 16,  
    parameter T_WIDTH = 18,  
    parameter T_RANGE = 10,
    parameter K_RANGE = 1000,
    parameter DATA_FILE = "data.mem"
)
(
    input wire clk_in,
    input wire rst_in,
    output reg done,
    //output reg [PARA_WIDTH-1:0]best_para
    output wire [5:0]sel, 
    output wire [7:0]seg          
);

wire clk,rst;
wire [31:0]random_num;
wire [PARA_WIDTH-1:0] para;
wire [DATA_WIDTH-1:0] data;
wire [INPUTDATA_WIDTH-1:0]input_data;
wire fit;
wire [T_WIDTH-1:0]RSS;
reg random_clk;
reg [ADDR_WIDTH-1:0] addr,fit_count,fit_max;
reg [T_WIDTH-1:0]rss_num,rss_min;
reg [15:0]k_num;
reg [PARA_WIDTH-1:0]best_para;
wire [23:0]output_data;
wire [5:0]point;

assign rst=~rst_in;

always@(posedge clk or posedge rst)
if(rst==1'b1)
    addr<=1'b0;
else
    addr<=addr+1;
    
always@(posedge clk or posedge rst)
if(rst==1'b1)
    k_num<=1'b0;
else if(k_num<=K_RANGE && addr==1'b0)
    k_num<=k_num+1;
    
always@(posedge clk or posedge rst)
if(rst==1'b1)
    random_clk<=1'b0;
else if(k_num<=K_RANGE && addr==1'b0)
    random_clk<=1'b1;
else 
    random_clk<=1'b0;
   
always@(posedge clk or posedge rst)
if(rst==1'b1)
    done<=1'b0;
else if(k_num>K_RANGE)
    done<=1'b1;
else 
    done<=1'b0;

always@(posedge clk or posedge rst)
begin
    if(rst==1'b1)begin
        rss_num<='b0;
        fit_count<='b0;
        fit_max<='b0;        
        
        best_para<='b0;
     end
     else if(k_num<=K_RANGE)begin   
        if(fit==1'b1)begin
            fit_count<=fit_count+1;
            rss_num<=rss_num+RSS;
        end     
        
        if(addr==1'b0)begin
            if(fit_count>fit_max)begin
                best_para<=para;
                fit_max<=fit_count;
                rss_min<=rss_num;
            end
            rss_num<='b0;
            fit_count<=1'b0;
        end
    end
end

assign output_data={8'b0,best_para};
assign point=6'b0;

clk_wiz_0 clk_wiz_0_inst
   (
    // Clock out ports
    .clk_out1(clk),     // output clk_out1
   // Clock in ports
    .clk_in1(clk_in)      // input clk_in1
);

blk_mem_gen_0 blk_mem_gen_0_inst1 (
  .clka(clk),    // input wire clka
  .addra(addr),  // input wire [6 : 0] addra
  .douta(data)  // output wire [15 : 0] douta
);

blk_mem_gen_0 blk_mem_gen_0_inst2 (
  .clka(clk),    // input wire clka
  .addra(random_num[6:0]),  // input wire [6 : 0] addra
  .douta(input_data[INPUTDATA_WIDTH-1:INPUTDATA_WIDTH/2])  // output wire [15 : 0] douta
);

blk_mem_gen_0 blk_mem_gen_0_inst3 (
  .clka(clk),    // input wire clka
  .addra(random_num[13:7]),  // input wire [6 : 0] addra
  .douta(input_data[INPUTDATA_WIDTH/2-1:0])  // output wire [15 : 0] douta
);

evaluate__kernel
#(
    .DATA_WIDTH(DATA_WIDTH),
    .PARA_WIDTH(PARA_WIDTH),
    .T_WIDTH   (T_WIDTH   ),
    .T_RANGE   (T_RANGE   ) 
)
evaluate__kernel_inst
(
    .data(data),
    .para(para),
    .RSS(RSS),
    .fit(fit)
);

calculate_kernel
#(
    .INPUTDATA_WIDTH(INPUTDATA_WIDTH),
    .PARA_WIDTH(PARA_WIDTH)
)
calculate_kernel_inst
(
    .input_data(input_data),
    .para(para)
);

random_num random_num_inst
(
    .clk   (random_clk)  ,
    .rst   (rst   )  ,
    .data  (random_num)
);

seg_dynamic seg_dynamic_inst
(
    .sys_clk  (clk      )   ,
    .sys_rst_n(rst      )   ,
    .data     (output_data)   , 
    .point    (point),
    .sel      (sel      )   ,
    .seg      (seg      )    
);

endmodule

七.参考文献

[1]束礼宝,宋克柱,王砚方.伪随机数发生器的FPGA实现与研究[J].电路与系统学报,2003(03):121-124.

[2]江洁,凌思睿.一种投票式并行RANSAC算法及其FPGA实现[J].电子与信息学报,2014,36(05):1145-1150.

八.设计改进

1.不可综合语句

        原代码中夹杂着不可综合语句，如rom中的initial语句。现将数据存储器使用xilinx block memory generator ip核替换原本不可综合的rom模块。

2.修正程序语句

（1）原程序出现阻塞赋值和非阻塞赋值语句混用情况，现对evaluate、evaluate_kernel、calculate_kernel相关语句进行了修正。（说明：控制模块evaluate和随机数生成模块random_num采用时序逻辑，计算和评估模块_kernel仍使用组合逻辑）

（2）原程序的控制模块evaluate中语句较为混乱，现依据信号划分语句，使结构更为清晰。

3.时钟模块

        添加xilinx clocking wizard ip核代替原硬件时钟直连，方便后续进行时序分析和调试。

九.设计实现分析

1.静态时序分析

        在vivado设计实现中，观察“Timing Summary”内的时钟域内路径“intra-clock paths”以确定时钟频率。可以看到路径上的时序要求“requirement”为30.303ns，数据信号在时钟周期内稳定下来的最长时间在path1上，如下图：

        该路径是从存储器blk_mem_gen0至寄存器fit_count_reg，即对应从随机从存储器中取数至模型遍历评估一个数结束。

        按照信号稳定的最长时间为30.303ns，计算出的最大频率约为33.00MHz。在时钟约束中设置时钟的输出频率为33MHz，implement后的结果如下：

        可以看到系统的最差负时序裕量“WNS”和最差保持时序裕量“WHS”均为正；总的负时序裕量“TNS”和保持时序裕量“THS”为零，即所有的时序路径都能满足这些要求。

十.下载验证

        在工程中添加数码管显示模块，并生成比特流对工程进行下载验证。使用开发板芯片型号为xc7a35tfgg484-2，资源使用情况如下：

        实际运行情况如下：