基于伪随机数LFSR和Box Muller变换的可综合的高斯白噪声FPGA生成器，AD/DA回环输出采集HDMI显示（已通过matlab、仿真和板级验证）_fpga产生白噪声模块

硬件平台

基于XLINX公司生产的AX7035开发板，具有HDMI输出输出，可以满足在没有示波器条件下输入输出回环测试。项目中仅使用了ROM ip核用来存储查找表计算根号、对数、cos、sin，可以移植到其他任意开发中，但HDMI输出波形可能无法观测到，只能通过示波器显示。

设计内容

设计内容主要分为两部分：高斯分布序列产生和HDMI显示。该项目侧重点是高斯白噪声产生，我主要介绍LFSR序列发生器和Box Muller转换设计思路。

LFSR伪随机数生成

该模块产生32位均匀分布序列，循环周期是2^64 = 1.8*10^19 。

利用64位斐波那契型LFSR，反馈多项式为 x^64 + x^63 + x^61 + x^60 + 1，可通过调用函数改变初始种子seed参数以产生不同的随机序列。

同时修改参数中的MIN和MAX大小设置输出随机序列的范围，理论上可以设置为-2,147,483,648——+2,147,483,648，也就是±2^31。

module LFSR_fibonacci #(
    parameter seed = 64'd18446744070000000000,
    parameter MIN = 1,
    parameter MAX = 100
)(
  input wire clk,
  input wire rst_n,
  input wire enable,
  output wire  [31:0] rand_out
);
  wire signed [63:0] rand64;
  wire signed [31:0] rand_out1;
  wire signed  [31:0] rand_out2;
  reg [63:0] lfsr_reg;
 
  always @(posedge clk or posedge rst_n) begin
    if (!rst_n) begin
      lfsr_reg <= seed; // 初始值可以根据需要设置为其他值
    end else if (enable) begin
      // LFSR 反馈多项式为 x^64 + x^63 + x^61 + x^60 + 1
      lfsr_reg <= {lfsr_reg[62:0], lfsr_reg[63] ^ lfsr_reg[60] ^ lfsr_reg[61] ^ lfsr_reg[63]}; 
    end
  end
 
  assign rand64 = lfsr_reg;
  assign rand_out1 = rand64[47:16] + {rand64[15:0],rand64[63:48]};
  assign rand_out2 = rand64[31:0] + rand64[63:32];
  assign rand_out = (rand_out1 + rand_out2) % (MAX - MIN+1'b1) + MIN; 
endmodule

注意，我将LFSR设置为64位，而输出是经过64位随机序列变换得到，因此产生的随机序列周期是64位线性移位寄存器产生的伪随机数周期，但输出位数为32位。如果有对于资源要求比较高可以减小线性移位寄存器位数，通过查找响应位数特征多项式（不同位数有不同的特征多项式，这决定了产生随机序列的周期，使用特征多项式可以达到最大周期2……N-1）。

对于输出变换的解释，由LFSR直接产生的随机序列在一段时间内小数特别多，不符合我们对于均匀随机序列的要求，具体可编写测试文件观测到。因此，我对输出64位随机序列进行变换，这里我用两个rand_out1和rand_out2相加的结果取模，rand_out1和rand_out2分别是rand64取不同位数得到。当然，这里的方法不唯一，感兴趣的可以自行设计。

同时，需要注意的是，在输出级采用MAX和MIN对输出范围做了限制，输出的随机序列最大最小值均是有符号数！想要获得任意范围定点数可以自行另外设计输出逻辑。

更详细的LFSR设计原理和结构可以参照

https://cloud.tencent.com/developer/article/2287083?areaId=106001

Box-Muller转换产生高斯分布序列

Box-Muller transformation算法有两种形式，一种称之为基本形式：

其中U1，U2是均匀分布的随机数，Z1和Z2是满足高斯分布的随机数。可以看出，要得到高斯分布序列需要两个均匀分布序列，同时经过算数运算得到。

另一种表达形式是极坐标形式我没不使用，原因是极坐标形式Box-Muller转换虽然不涉及cos、sin运算，但引入了除法，这在FPGA设计中需要占用较多逻辑资源，我就使用基本形式对其变换。

利用使用ROM查找表，设计ln、根号、正余弦计算器。根据输入变量大小实现地址映射，考虑到设计精度等因素，对输入输出进行一定比例放大。对于cos只需要存储1/4周期的数据，其余数据均可通过地址映射得到这里不展开细讲。

module cos(
input                       clk,
input        [15:0]         x_in,
input                       enable,
 
output      [15:0]         result_out    
    );
wire [9:0]   addr_ROM;
wire [13:0]  ROM_out;
cos_ROM cos_ROM_0(
    .clka       (~clk),
    .ena        (enable),
    .addra      (addr_ROM),
 
    .douta      (ROM_out)
);
wire  [15:0]   angle_3600;
 
assign angle_3600 = (x_in[15] == 1'b1) ? ((~x_in[15:0]+1'b1) % 16'd3600) : (x_in[15:0] % 16'd3600); 
wire            sign;
wire    [15:0] addr_angle;
assign addr_angle = (angle_3600 <= 16'd900)     ?   (angle_3600           )  :
                    (angle_3600 < 16'd1800)     ?   (16'd1800 - angle_3600)  :
                    (angle_3600 < 16'd2700)     ?   (angle_3600 - 16'd1800)  :
                    (angle_3600 < 16'd3600)     ?   (16'd3600 - angle_3600)  : 16'd901;
assign sign       = (angle_3600 <= 16'd900)     ?   (1'b0                 )  :
                    (angle_3600 < 16'd1800)     ?   (1'b1                 )  :
                    (angle_3600 < 16'd2700)     ?   (1'b1                 )  :
                    (angle_3600 < 16'd3600)     ?   (1'b0                 )  : 1'b0;
assign addr_ROM = addr_angle[9:0];
reg sign_delay0;
always @(posedge clk) begin
    sign_delay0 <= sign;
end
assign result_out = sign_delay0 ? (~{2'b0,ROM_out} + 1'b1): {2'b0,ROM_out};  //
endmodule

输入是角度制，同时需要将角度放大10倍，也就是3600是一个周期。输出16位有符号数，为是结果放大10000倍。

另外，是利用流水线的计算模式，ln和cos是第一级流水线，ln计算的结果乘以-2后再取根号，因此需要把cos计算的结果延迟一个节拍，再相乘输出。

Box Muller 转换流水线代码如下

module Box_Muller_transform(
  input wire clk_data,
  input wire clk_cal,
  input wire enable,
  input wire signed [15:0] data_in_1,
  input wire signed [15:0] data_in_2,
 
  output reg signed [31:0] result_out_0,
  output reg signed [31:0] result_out_1
);
  wire signed [15:0] cos_sin_in;
  wire signed [15:0]  sqrt_in;
  wire signed [15:0] sqrt_out, ln_out, cos_out, sin_out;
  reg signed [15:0] cos_out_delay, sin_out_delay;
  
  reg signed [15:0] x_in_1, x_in_2;
  always @(posedge clk_data) begin
    if (enable) begin
      x_in_1 <= data_in_1;
      x_in_2 <= data_in_2;
    end
  end
  // instantiate sqrt module
  assign  sqrt_in = ln_out * (-16'd2);
  sqrt sqrt_inst (
    .clk(clk_cal),
    .enable(enable),
    .x_in(sqrt_in),
    .result_out(sqrt_out)
  );
  // instantiate ln module
  ln ln_inst (
    .clk(clk_cal),
    .enable(enable),
    .x_in(x_in_1),
    .result_out(ln_out)
  );
  // instantiate cos module
  assign cos_sin_in = x_in_2 * 16'd36 / 16'd10;
  cos cos_inst (
    .clk(clk_cal),
    .enable(enable),
    .x_in(cos_sin_in),
    .result_out(cos_out)
  );
   sin sin_inst (
    .clk(clk_cal),
    .enable(enable),
    .x_in(cos_sin_in),
    .result_out(sin_out)
  ); 
  // delay cos result
  always @(posedge clk_cal) begin
    if (enable) begin
      cos_out_delay <= cos_out;
      sin_out_delay <= sin_out;
    end
  end
  // calculate final result using your formula
  always @(posedge clk_data) begin
    if (enable) begin
      result_out_0 <= cos_out_delay * sqrt_out;
      result_out_1 <= sin_out_delay * sqrt_out;
    end
  end
endmodule

高斯序列产生

该模块例化了两个LFSR随机数生成模块（利用不同的种子）和一个Box Muller转换模块，将两个均匀分布序列转换为高斯分布序列。

同时，可以在模块参数中修改均值和方差，以达到产生任意分布正态分布的作用。默认为标准正态分布。

/默认产生32位有符号数，结果除以1M为实际值
module Guassian_generate #(
parameter average = 32'd0, //均值
parameter variance = 32'd0  //方差
)(
//input                   sys_clk,
input                   clk_data,
input                   clk_tra,
input                   rst_n,
input                   enable,
 
debug
//output  wire [31:0]         rand_out_1,
//output  wire [31:0]         rand_out_2,
 
output  wire signed [31:0]         rand_out1,
output  wire signed [31:0]         rand_out2
    );
    
wire [31:0] rand_out_1, rand_out_2;
LFSR_fibonacci #(
    .seed(64'd18446744070000000000),
    .MIN(32'd1),
    .MAX(32'd1000)
) Uniform_Distribution_generator_0(
    .clk(clk_data),
    .rst_n(rst_n),
    .enable(1'b1),
    .rand_out(rand_out_1)
);
LFSR_fibonacci #(
    .seed(64'd14742147901846518232),
    .MIN(32'd1),
    .MAX(32'd1000)
) Uniform_Distribution_generator_1(
    .clk(clk_data),
    .rst_n(rst_n),
    .enable(1'b1),
    .rand_out(rand_out_2)
);
//将均匀分布输出数据转化为16位输入
wire [15:0] tra_data_in_1, tra_data_in_2;
assign tra_data_in_1 = rand_out_1[15:0];
assign tra_data_in_2 = rand_out_2[15:0];
 
wire signed [31:0]         Guassian_0;
wire signed [31:0]         Guassian_1;
Box_Muller_transform Box_Muller_transform(
    .clk_data(clk_data),
    .clk_cal(clk_tra),
    .enable(1'b1),
    .data_in_1(tra_data_in_1),
    .data_in_2(tra_data_in_2),
    
    .result_out_0(Guassian_0),
    .result_out_1(Guassian_1)
);
assign rand_out1 = enable ? ((Guassian_0[31] == 1'b1) ? (~((~Guassian_0 + 1'b1) / variance ) + 1'b1) + average : Guassian_0 / variance + average) : 32'd0;
assign rand_out2 = enable ? ((Guassian_1[31] == 1'b1) ? (~((~Guassian_1 + 1'b1) / variance ) + 1'b1) + average : Guassian_1 / variance + average) : 32'd0;
endmodule

注意该模块默认产生32位有符号数，结果除以1000000为实际值。这里实际值是标准正态分布的结果，因为根据Box Muller转换，输入均匀分布的范围是0-1，我利用LFSR产生1-1000的随机数，相当与输入放大了1000倍。

仿真测试

测试LFSR伪随机数生成均匀分布序列

使用vivado内置仿真环境

matlab仿真

延长仿真时间生成10M个样本

LFSR仿真测试代码和MATLB验证代码同下正态高斯测试一样，只需稍加修改即可。

高斯序列产生测试Verilog代码

`timescale 1ns / 1ps
//
// Company: 
// Engineer: 
// 
// Create Date: 2024/01/14 12:37:28
// Design Name: 
// Module Name: Guassian_generate_tb
// Project Name: 
// Target Devices: 
// Tool Versions: 
// Description: 
// 
// Dependencies: 
// 
// Revision:
// Revision 0.01 - File Created
// Additional Comments:
// 
//
 
 
module Guassian_generate_tb;
reg clk_data;
reg clk_tra;
reg rst_n;
reg enable;
 
wire signed [31:0] rand_out1;
  Guassian_generate#(
    .average(32'd0),
    .variance(32'd0)
  )Guassian_generate_inst (
    .clk_data(clk_data),
    .clk_tra(clk_tra),
    .rst_n(rst_n),
    .enable(enable),
    
    .rand_out1(rand_out1)
  );
wire [7:0] rand;
assign rand = rand_out1[7:0];
// Clock generation
initial begin
  clk_tra = 0;
  forever #3 clk_tra = ~clk_tra;
end
 
initial begin
  clk_data = 0;
  #7;
  forever #10 clk_data = ~clk_data;
end
integer Guss0;
//integer Guss1;
initial begin
    Guss0=$fopen("randGuss0.txt");    //打开所创建的文件
//    Guss1=$fopen("randGuss1.txt");    //打开所创建的文件 
//    rand0=$fopen("rand0.txt");    //打开所创建的文件
//    rand1=$fopen("rand1.txt");    //打开所创建的文件 
end
 
always @(negedge clk_data) begin
     if(enable) begin        
       $fdisplay(Guss0,"%d",$signed(rand_out1));    //$fdisplay(dout_file,"%d",$signed(num)); //保存有符号数据
//       $fdisplay(Guss1,"%d",$signed(Z_standard1));    //$fdisplay(dout_file,"%d",$signed(num)); //保存有符号数据
//       $fdisplay(rand0,"%d",$signed(rand_out_1));    //$fdisplay(dout_file,"%d",$signed(num)); //保存有符号数据
//       $fdisplay(rand1,"%d",$signed(rand_out_2));    //$fdisplay(dout_file,"%d",$signed(num)); //保存有符号数据
     end
end
// 在仿真结束时关闭文件
initial begin
    #10000000; // 确保仿真足够长的时间
    $fclose(Guss0);
//    $fclose(Guss1);
//    $fclose(rand0);
//    $fclose(rand1);
    $stop;
end
// Testbench stimulus
initial begin
  rst_n = 0;
  enable = 0;
 
 #100
 rst_n =1;
 
 #100
 enable = 1;
  // Finish simulation
  #10000000;
end
endmodule

将测试文件生成的随机数存储在randGuss0.txt文件中，再通过matlab读取进行分析。

可以看出产生的随机序列为一个标准的正态分布！

matlab测试代码如下。注：这里没有使用绝对路径查看txt文件，'randGuss0.txt'使用相对路径保存的文件再sim目录下，请自行查找。

%生成rand64向量
vivado_data_Guss0  =  importdata('randGuss0.txt');  %read files
% 统计数据分布
minValue = min(vivado_data_Guss0);
maxValue = max(vivado_data_Guss0);
meanValue = mean(vivado_data_Guss0);
medianValue = median(vivado_data_Guss0);
stdDeviation = std(vivado_data_Guss0);
 
% 显示统计结果
fprintf('最小值：%f\n', minValue);
fprintf('最大值：%f\n', maxValue);
fprintf('平均值：%f\n', meanValue);
fprintf('中位数：%f\n', medianValue);
fprintf('标准差：%f\n', stdDeviation);
 
% 绘制直方图
figure(1)
histogram(vivado_data_Guss0, 100, 'Normalization', 'probability');
title('guss1数据分布直方图');
xlabel('数值');
ylabel('频率');

DC输出高斯分布序列AD采集，HDMI显示