高斯分布的随机数生成器_box muller fpga

高斯分布的随机数生成器

	实现的过程是先查找相关高斯分布随机数在vivado实现的博客，先大概认识一下，然后到知网找相关的硕士论文，总结出最简单的高斯随机数生成的实现方法，在进行仿真验证。
1

在查阅相关论文后把高斯分布随机数生成器分为两个模块：
模块一：基于细胞自动机生成均匀随机数，具体是采用64阶细胞自动机的均匀随机数发生器来产生均匀随机数。
模块二：通过Box-Muller算法将均匀随机数转换成为高斯随机数。

模块一：通过细胞自动机生成均匀随机数

​ 关于细胞自动机的最早的思想起源于StanislawUlam，它是一种特殊的具有时间、空间和状态离散性的有限状态机。由于细胞自动机概念的复杂性，所以早期的对于细胞自动机的研究基本上是集中于理论方面，直到20世纪80年代初，由S.Wolfram46首先对细胞自动机的概念进行了简化，从而极大推动了细胞自动机理论及其应用研究。细胞自动机是D维空间中一组细胞单元组成的阵列，一个细胞自动机可用四元组表示为:
$$CA=(A_D,Z_q,f_i(o,r),B)$$
​
其中，空间结构A是由D维细胞单元构成的空间结构,状态空间-q是细胞自动机中细胞单元i的状态取值范围;邻域函数规则是细胞自动机的邻域半径r确定的第i个细胞单元的0阶邻域状态配置与其转移状态之间的映射: 边界条件B规定了某个细胞单元领域半径之内的邻域单元在超出细胞自动机空间结构时的处理方法。

​ 代码实现32阶细胞自动机的均匀随机数发生器来产生均匀随机数，通过给出32位细胞的初始值和规则号（90规则和150规则），将单个细胞例化32次，就可以得到32位的细胞自动机。

单个细胞如下

 module single_cell#(parameter init = 1'b0 , head = 1'b0, tail = 1'b0)
(
input clk_i,rst_n_i,
input ctrl,
input self,left,right,

output out
);

reg  self_next;

always @(*) begin
    case({head,tail})
             2'b10 : self_next = ctrl ? self ^ right : right;
             2'b01 : self_next = ctrl ? self ^ left  : left;
             2'b00 : self_next = ctrl ? left ^ self ^ right : left ^ right;
             default:  self_next = ctrl ? left ^ self ^ right : left ^ right;
    endcase
end

assign out = self_next; 

endmodule
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

例化过程：因为头尾细胞和中间细胞存在区别，分为三个部分分开例化，以中间细胞为例，只需将left ( uni_r[i-1]),更改为left (1’b0)即为头细胞，将.right (uni_r[i+1])更改为.right (1’b0)即为尾细胞

module mata
  #(parameter INTA_VEC = 32'b0100_1000_0001_0010_0100_1000_0001_0010,
    parameter RULE_VEC   = 32'b0000_1100_0100_0111_0000_1100_0000_0110,
    N = 32   )

   (
     input clk_i,rst_n_i,
     input pause_i,        //暂停
     output [N-1:0] out_1  //均匀随机数
   );

  reg  [N-1:0] uni_r;
  wire [N-1:0] uni_next;
  assign out_1 = uni_r;

  always @(posedge clk_i or negedge rst_n_i)
  begin
    if(~rst_n_i)
      uni_r <= INTA_VEC;
    else if(pause_i)
      uni_r <= uni_r;
    else
      uni_r <= uni_next;
  end

  generate
    genvar i;
    for (i = 0; i < N; i = i + 1 )
    begin
      if( i==0 )
      begin//第一位
        single_cell
          #(.init (INTA_VEC[i]),.head(1'b1 ),.tail (1'b0))
          single_cell_first_num (
            .clk_i (clk_i ),
            .rst_n_i (rst_n_i ),
            .ctrl (RULE_VEC[i] ),
            .left ( 1'b0),
            .right (uni_r[1]),
            .self (uni_r[i]),
            .out  ( uni_next[i])
          );
      end

      else if(i == N - 1)
      begin//最后一位
        single_cell
          #(.init (INTA_VEC[i] ),.head(1'b0 ),.tail (1'b1))
          single_cell_num (
            .clk_i (clk_i ),
            .rst_n_i (rst_n_i ),
            .ctrl (RULE_VEC[i] ),
            .left ( uni_r[i-1]),
            .right (1'b0),
            .self (uni_r[i]),
            .out  ( uni_next[i])
          );
      end


      else
      begin //中间位
        single_cell
          #(.init (INTA_VEC[i] ),.head(1'b0 ),.tail (1'b0))
          single_cell_mid_num (
            .clk_i (clk_i ),
            .rst_n_i (rst_n_i ),
            .ctrl (RULE_VEC[i] ),
            .left ( uni_r[i-1]),
            .right (uni_r[i+1]),
            .self (uni_r[i]),
            .out  ( uni_next[i])
          );
      end
    end
  endgenerate

endmodule
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78

对32阶的细胞自动机进行仿真

给出32位细胞的初始值和规则号，给定时钟，采集数据写入new_data_b.txt文本，方便接下来在matlab上进行仿真。

module tb_celluar_auto();
  // Parameters
  localparam  INTA_VEC = 32'b0100_1000_0001_0010_0100_1000_0001_0010;
  localparam  RULE_VEC = 32'b0000_1100_0100_0111_0000_1100_0000_0110;
    
  localparam  N = 32;
  localparam  period = 10; //100Mhz

  // Ports
  reg clk_i = 0;
  reg rst_n_i = 1;
  reg pause_i = 0;
  wire [N-1:0] out_1;

  mata 
  #(
    .INTA_VEC(INTA_VEC ),
    .RULE_VEC(RULE_VEC ),
    .N (N)
  )
  mata_dut (
    .clk_i (clk_i ),
    .rst_n_i (rst_n_i ),
    .pause_i (pause_i ),
    .out_1  ( out_1)
  );

parameter M = 200_000;
reg [N-1:0] men [M:1];
integer index = 1;
initial begin //复位
    begin
      #(period*2) rst_n_i = 1'b0;
      #(period*10) rst_n_i = 1'b1;

      #(period*201_000);
      $writememb("new_data_b.txt",men);
      $finish;
    end
  end

  initial //采集数据
    begin
        #(period*100);
        forever begin
            #(period);
            men[index] = out_1;
            index = (index >= M) ? 1 : index + 1;
        end
    end
  
  always //生成时钟
    #(period)  clk_i = ! clk_i ;

endmodule
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

将txt文本的数据导入到matlab上可以看到仿真结果基本符合在0到1之间均匀分布。
附上matlab的代码

clc;
fid = fopen('D:\BaiduNetdiskDownload\new_data_b.txt');
data = textscan(fid, '%b');
fclose(fid);
data = cell2mat(data);
data = double(data) / double(max(data));
%29 -> [0,1] -> data / max(data);
 %h = histogram(data,16);
h = histogram (data, 100,'BinLimits',[0,1]);
1
2
3
4
5
6
7
8
9

模块一结束下面进行模块二，Box-Muller算法部分

​ Box-Muller算法是最早的产生高斯白噪声的算法之一，它是利用均匀随机数来分别计算出高斯随机数的幅度和相位值从而产生高斯随机数的算法， Box-Muller算法可以同时将两个均匀随机数转换成为两个高斯随机数。其隐含的原理非常深奥，但结果却是相当简单。假设x和x为相互独立的零均值的服从标准正态分布的高斯随机变量，那么可知R=√x2+x20=arctan(x2/x)分别服从瑞利分布和均匀分布。利用这一事实，可以通过将一对服从瑞利分布的随机变量和一对服从均匀分布的随机变量进行变换，产生出一对高斯变量来。具体算法如下:

​ 其中，u1和u2是在[01]上均匀分布的均匀随机变量，u1用来产生高斯随机数的幅度值，u2则产生高斯随机数的相位值。因为算法涉及自定义对数函数和三角函数，对于其中的自定义对数函数我采用查找表的方式实现，首先在matlab里生成对数函数对应自变量和因变量映射关系存放在ROM IP核里，通过查找自变量得到因变量。

使用matlab产生数据并生成.coe文件

clear all;
clc;

N = 16;
out_N = 16;
x = linspace(0.000001,1,pow2(N));
y = sqrt(-2*log(x));
y_fix = round(y/max(y)*(2^out_N));
y_fix(1) = y_fix(1) - 1;
y_qua = y_fix;
width = 16;
depth = pow2(N);

plot(x,y_qua);

fid = fopen("f.coe","w");   % 创建coe文件;
fprintf(fid,"memory_initialization_radix=16;\n");
fprintf(fid,"memory_initialization_vector=\n");

for num = 0 : (depth-1)
    fprintf(fid,"%02X,\n",y_qua(num+1));
end

fseek(fid,-2,1);
fprintf(fid,";");

fclose(fid);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

添加ROM IP核
在IP Catalog里搜索ROM，选择Bloke Memory Generator
在IP核配置界面basic选择单端输入rom即Single Port ROM

将matlab生成的coe文件导入到ROM IP核中，并对IP核进行命名,注意.coe文件的路径不能有中文，否者会标红。

添加DDS IP核
三角函数用DDS IP核实现，下图包含IP核配置情况

接下来再添加乘法器的ip核，并且命名。

在IP Sources里可以看到已添加的IP核

例化IP核，完成Box-Muller算法部分。

module box_muller(
  input clk_i ,rst_n_i,//复位
  input [63:0] uni_i,

  output [31:0] gau1_o,
  output [31:0] gau2_o
    );

    wire [15:0] uni_g= uni_i [63:48];
    wire [15:0] uni_f= uni_i [31:16];
    wire [15:0] fun_f_o;
    wire [15:0] g1_sin_o;
    wire [15:0] g2_cos_o;
    wire m_axis_data_tvalid;
    //wire [15:0] g1_sin_o,g2_cos_o;

  dds_g12 dds_g12_inst (
  .aclk(clk_i),                                // input wire aclk
  .aresetn(rst_n_i),                          // input wire aresetn
  .s_axis_phase_tvalid(1'b1),  // input wire s_axis_phase_tvalid
  .s_axis_phase_tdata(uni_g),    // input wire [15 : 0] s_axis_phase_tdata
  .m_axis_data_tvalid(m_axis_data_tvalid),    // output wire m_axis_data_tvalid
  .m_axis_data_tdata({g1_sin_o,g2_cos_o})      // output wire [31 : 0] m_axis_data_tdata
          );
  bram_fun_f bram_fun_f_inst (                           //ROM
  .clka(clk_i),    // input wire clka
  .addra(uni_f),  // input wire [15 : 0] addra
  .douta(fun_f_o)  // output wire [15 : 0] douta
);

  multi_box_muller multi_box_muller_inst_1 (             //两个乘法器
    .CLK(clk_i),  // input wire CLK
    .A(g1_sin_o),      // input wire [15 : 0] A
    .B(fun_f_o),      // input wire [15 : 0] B
    .P(gau1_o)      // output wire [31 : 0] P
);

  multi_box_muller multi_box_muller_inst_2 (
    .CLK(clk_i),  // input wire CLK
    .A(g2_cos_o),      // input wire [15 : 0] A
    .B(fun_f_o),      // input wire [15 : 0] B
    .P(gau2_o)      // output wire [31 : 0] P
);

endmodule
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

接下来进行仿真，例化Box-Muller算法模块。

module tb_box_muller();
  // Parameters
  // Ports
  reg clk_i = 0;
  reg rst_n_i = 1;
  reg [63:0] uni_i = 64'd0;
  wire [15:0] g1_sin_o;
  wire [15:0] g2_cos_o;
  wire [15:0] fun_f_o;
  
  localparam  period = 10; 

  box_muller 
  box_muller_dut (
    .clk_i (clk_i ),
    .rst_n_i (rst_n_i ),
    .uni_i (uni_i ),
    .g1_sin_o (g1_sin_o ),
    .g2_cos_o (g2_cos_o ),
    .fun_f_o  ( fun_f_o)
  );

  initial begin
    begin
        #(period*2) rst_n_i = 1'b0;
        #(period*10) rst_n_i = 1'b1;
        #(period*201_000);
      $finish;
    end
  end

  always
  #(period) uni_i = uni_i +  64'd1;

  always
  #(period/2)  clk_i = ! clk_i ;

endmodule
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

然后进行Box-Muller算法模块的仿真，仿真结果如图可以看到自定义对数函数和三角函数均已实现

实现高斯随机数的生成

创建顶层模块Top例化细胞自动机和Box-Muller算法模块

module gauss_top(
    (* io_buffer_type = "none" *)input clk_i ,rst_n_i,//复位

    (* io_buffer_type = "none" *)output [31:0] gau1_o,
    (* io_buffer_type = "none" *)output [31:0] gau2_o
    );

    wire [31:0] out_1;
    wire [31:0] out_2;
    mata
    #(
      .INTA_VEC(32'b0000_1100_0100_0111_0000_1100_0000_0110 ),
      .RULE_VEC(32'b0000_1100_0100_0111_0000_1100_0000_0110 ),
      .N ( 32 )
    )
    mata_u1 (
      .clk_i (clk_i ),
      .rst_n_i (rst_n_i ),
      .pause_i (1'b0 ),
      .out_1  ( out_1)
    );

    mata 
    #(
      .INTA_VEC(32'b0100_0110_0000_1001_1011_1011_1101_0101 ),
      .RULE_VEC(32'b0100_0110_0000_1001_1011_1011_1101_0101 ),
      .N ( 32 )
    )
    mata_u2 (
      .clk_i (clk_i ),
      .rst_n_i (rst_n_i ),
      .pause_i (1'b0 ),
      .out_1  ( out_2)
    );

    box_muller 
    box_muller_dut (
      .clk_i (clk_i ),
      .rst_n_i (rst_n_i ),
      .uni_i ( {out_1,out_2} ),
      .gau1_o (gau1_o ),
      .gau2_o  ( gau2_o)
    );
endmodule
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

将细胞自动机模块与Box-Muller算法模块例化到tb_top, 将采集仿真的数据g1 g2并且导入到new_data_g1.txt new_data_g2.txt文本里。

module tb_top();
// Parameters
localparam  period = 10; //100Mhz
// Ports
  reg clk_i = 0;
  reg rst_n_i = 1;
  wire [31:0] gau1_o;
  wire [31:0] gau2_o;
  gauss_top 
  gauss_top_dut (
    .clk_i (clk_i ),
    .rst_n_i (rst_n_i ),
    .gau1_o (gau1_o ),
    .gau2_o  ( gau2_o)
  );

  parameter M = 200_000;
  parameter N = 32;

  reg [N-1:0] mem1 [M:1];
  reg [N-1:0] mem2 [M:1];
  integer index = 1;
  initial begin //复位
    begin
      #(period*2) rst_n_i = 1'b0;
      #(period*10) rst_n_i = 1'b1;

      #(period*201_000);
      $writememb("new_data_g1.txt",mem1);
      $writememb("new_data_g2.txt",mem2);
      $finish;
    end
  end

  initial //采集数据
  begin
      #(period*100);
      forever begin
          #(period);
          mem1[index] = gau1_o;
          mem2[index] = gau2_o;
          index = (index >= M) ? 1 : index + 1;
      end
  end

  always
  #(period/2)  clk_i = ! clk_i ;
endmodule
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

最后将txt文本里的数据导入到maltab上验证一下，可以看到仿真结果基本是符合了高斯分布的随机数。

附上matlab验证代码

clear all;
clc;

fid = fopen('D:\FPGA\Vivado projects\winter exercise\project_matlab\new_data_g1.txt');
data = textscan(fid,'%bs32');
fclose(fid);
data = cell2mat(data);

data = double(data)/pow2(31);
num_point = 1000;

% h = histogram(data,16);
h = histogram(data,num_point,'BinLimits',[-1,1]);

counts = hist(data,num_point,'BinLimits',[-1,1]);
counts = counts/max(counts);

f = fittype('a*exp(-((x-b)/c)^2)');
x = linspace(-1,1,num_point);
y = counts;
plot(x,y,'.')

[cfun,gof] = fit(x(:),y(:),f);
yy = cfun.a*exp(-((x-cfun.b)/cfun.c).^2);
hold on;
plot(x,yy,'r','LineWidth',2);
A = cfun.a;%幅值
mu = cfun.b;%期望
sigma = cfun.c;%标准差
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

![img](https://img-blog.csdnimg.cn/2a0f859d02a14e5a9e2c8fab0dc6d367.png?x-oss-process=image/watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBAd2VpeGluXzUzNzEyMTQ4,size_19,color_FFFFFF,t_70,g_se,x_16