目标反射回波检测算法及其FPGA实现之三：
平方、积分电路及算法的顶层FPGA实现
前段时间，接触了一个声呐目标反射回波检测的项目。声呐接收机要实现的核心功能是在含有大量噪声的反射回波中，识别出发射机发出的激励信号的回波。我会分几篇文章分享这个基于FPGA的回波识别算法的开发过程和原码，欢迎大家不吝赐教。以下原创内容欢迎网友转载，但请注明出处： https://www.cnblogs.com/helesheng。
在本系列博文的第一篇中，根据仿真结果，我认为采用“反射回波和激励信号互相关”的结果来计算目标距离的方法具有较高性能和计算效率。在本系列的第二篇博文中，我在Cyclone系列的低成本FPGA中采用半并行的“双存储器式的卷积节”结构实现了数据的互相关/卷积/FIR滤波器计算。作为本系列的第三篇博文，我将实现互相关信号的平方和积分计算，并将所有算法在顶层文件中结合为一个整体。
从而通过寻找 
的极值点所在位置来确定目标反射回波出现的时间点。
（1）式中的
是互相关算法部分，其FPGA实现已在前文中介绍过。根据前文定义的符号，将离散化后的互相关信号记为R[k]。进一步离散化后可将（1）改写为如下FPGA能够实现的形式：
一、平方电路的实现
使用Quartus-II中的MegaWizard配置平方计算电路，其结构如下图所示。
图1 平方电路配置
二、积分电路设计
根据（2）式，要计算目标函数P[n]的值，还需要对历史上的 值求和（积分）。我们当然可以用缓冲器存储历史上的k0个 值，并在每次结果输出之前对缓冲器中的k0个值求和。如果让k0等于激励信号的长度N，则每次输出P[n]之前都需要计算N-1个加法。当N为64时（如前文所述），这几乎是不可完成的任务。我设计了下图所示的电路结构来实现64个历史数据的求和。
图2 积分器电路结构
其中一位深度为64的移位寄存器的作用是提供64个采样之前的“历史数据”。首先对当前数据
和最老的历史数据
求差（补码），再对差不断求和。这样根据加法交换律，最终求和结果就相当于加入了当前数据，减去了最老的历史数据。只要移位寄存器的初值全为0，在进行了64次操作后，从输出得到的就是64个历史数据的和（积分）了。其中用MegaWizard配置的移位寄存器结构如下图所示。
图3 移位寄存器的结构
前级减法器的结构如下图所示。
图4 减法器的结构
累加器采用硬件描述语言实现，代码如下。
1 always @ (posedge start or negedge rst_n) 2 begin 3 if(!rst_n) 4 acc_out[39:0] <= 40'd0; 5 else begin 6 acc_out[39:0] = $signed(acc_out[39:0]) + $signed(acc_in[31:0]); 7 end 8 end
其中的关键字$signed表示有符号数的加法器。
三、算法的顶层设计
为了将前述的A/D和D/A口控制电路、互相关/卷积/FIR滤波电路、平方和积分电路连接为一个系统，还需要在顶层设计文件中对上述模块电路进行例化和连接。另外顶层设计文件还将对系统的整体工作时序进行控制。我设计的顶层文件如下所示。
1 module CONV_POW_AD_DA( 2 ///////////顶层模块，负责调用ADC采集数据，卷积，然后用DAC输出卷积/FIR的结果///////// 3 input rst_n,//低电平复位信号 4 input iclk20,//外部晶体输入的20MHz 5 output sck_da,//D/A转换器的SPI时钟 6 output mosi_da,//D/A转换器的SPI数据信号 7 output cs_da,//D/A转换器的片选信号 8 output ld_da,//D/A转换器的双通道数据加载信号 9 output sck_ad,//A/D转换器的SPI时钟 10 output mosi_ad,//A/D转换器的SPI数据输出 11 input miso_ad,////A/D转换器的SPI数据输入 12 output cs_ad//A/D转换器的SPI口片选信号 13 ); 14 wire clk; 15 reg start; 16 reg[15:0] cnt;//用于产生总体周期的计时器 17 reg[11:0] tst_data;//用于产生测试数据的计数器 18 parameter CNT_NUM = 16'd2000;//100M时钟下，2000分频意味着50KHz溢出率 19 wire[11:0] ad_data;//AD转换结果的内部连线 20 wire[11:0] data_cha;//A通道数据 21 wire[11:0] data_chb;//B通道数据 22 wire[15:0] shft_data1;//在卷积的两个节之间传递的数据 23 wire[15:0] shft_data2;//在卷积的两个节之间传递的数据 24 wire[15:0] shft_data3;//在卷积的两个节之间传递的数据 25 wire[15:0] shft_data4;//在卷积的两个节之间传递的数据 26 wire[39:0] conv_res1;//第一节的卷积的结果 27 wire[39:0] conv_res2;//第2节的卷积的结果 28 wire[39:0] conv_res3;//第3节的卷积的结果 29 wire[39:0] conv_res4;//第4节的卷积的结果 30 wire[39:0] conv_res;//总的卷积的结果 31 wire[33:0] acc_sum;//最后计算每一节的累加和的结果 32 wire[31:0] shiftout;//移位寄存器的移出数据 33 wire[15:0] ac_sig;//去除直流偏置后的交流信号 34 wire[31:0] pow_sig;//信号的能量 35 wire signed[31:0] acc_in;//求和累加器的输入，也是加法/减法器的输出 36 reg signed[39:0] acc_out;//累加器输出的结果 37 wire signed[19:0] data_out;//累加结果开方后的输出 38 39 //assign data_cha[11:0] = acc_out[25:14];//通道a输出64个点能量累加结果 40 //assign data_cha[11:0] = shiftout[18:7];//通道a输出的数据 41 assign data_cha[11:0] = pow_sig[19:8];//通道a输出交流信号的平方 42 //assign data_cha[11:0] = ad_data[11:0];//通道a输出的数据 43 //assign data_cha[11:0] = shft_data2[11:0];//通道a输出的数据 44 //assign data_chb[11:0] = conv_res[26:15];//通道b输出的数据 45 //系数本身是16bit带符号的，由于FIR滤波器的系数之和为1，则16bits带符号的所有系数之和是32768， 46 //因此需要将滤波结果右移15bit，所以取结果的26:15。但由于DA参考电压为2.048V，小于AD的3.3V因此滤波结果还是略小。 47 //assign data_chb[11:0] = data_out[14:3];//通道b输出累加结果开方 48 assign data_chb[11:0] = acc_out[28:17];//通道b输出64个点能量累加结果 49 //assign data_chb[11:0] = pow_sig[20:9];//通道b输出交流信号的平方 50 //assign data_chb[11:0] = acc_sum[26:15];//通道b输出滤波/卷积结果 51 //assign data_chb[11:0] = conv_res1[26:15];//通道b输出的数据 52 53 /////以下系数配置中用到的连线///// 54 wire rdy_work;//配置完成信号，高电平表示配置完成 55 wire[7:0] wr_blk_addr;//向后级系数blk配置的地址 56 wire[15:0] data_bank;//向各级系数blk配置的数据 57 wire csh0,csh1,csh2,csh3;//每一节配置的系数blk的选通线 58 PLL20_100 i_pll20_100(//例化PLL产生时钟 59 .inclk0(iclk20), 60 .c0(clk)//pll输出的100M工作时钟 61 ); 62 63 always @ (posedge clk or negedge rdy_work) 64 begin 65 if(!rdy_work) 66 begin 67 start <= 1'd0; 68 cnt[15:0] <= 16'd0; 69 tst_data[11:0] <= 12'd0; //用作测试的计数器 70 end 71 else begin 72 //////////维护周期计数器/////////// 73 if(cnt < CNT_NUM-1) begin 74 cnt[15:0] <= cnt[15:0] + 16'd1; 75 tst_data[11:0] <= tst_data[11:0]; 76 end 77 else begin 78 cnt[15:0] <= 16'd0; 79 tst_data[11:0] <= tst_data[11:0] + 12'd1; //用作测试的计数器 80 end 81 82 ////////产生启动信号///// 83 if((cnt > 102)&(cnt <= 105))//在第105个clk产生启动信号 84 start <= 1'D1;//启动MCP4822输出状态机 85 else 86 start <= 1'D0; 87 end 88 end 89 90 init_coe_blk i_init_coe_blk(//例化系数配置模块，起到从系数池中读取系数并向系数blk中配置数据的作用，只在复位后的开始阶段有效，通过ready信号控制后续模块 91 .clk(clk), 92 .rst_n(rst_n),//整体复位信号 93 .ready(rdy_work),//配置完成信号，高电平表示配置完成 94 .wr_blk_addr(wr_blk_addr),//向后级系数blk配置的地址 95 .data_bank(data_bank),//向各级系数blk配置的数据 96 .csh0(csh0),//第一节配置的系数blk的选通线 97 .csh1(csh1),//第二节配置的系数blk的选通线 98 .csh2(csh2),//第三节配置的系数blk的选通线 99 .csh3(csh3)//第四节配置的系数blk的选通线 100 ); 101 102 CONV_SER16 i_conv_ser16_I(//例化第一个16阶卷积/fir电路 103 .clk(clk), 104 .a({4'd0,ad_data}),//AD转换结果作为数据输入 105 //.a({4'd0,tst_data[11:0]}),//AD转换结果作为数据输入 106 .en(rdy_work),//系数配置完成后才能使能 107 .coe_data_in16(data_bank),//初始化系数的数据输入端 108 .wr_coe_addr(~wr_blk_addr[3:0]),//初始化系数的地址输入端 109 ////!!!特别注意这里，乘加时数据从大地址进入乘加操作，因此地址要求补码，将系数翻转过来。相当于：fliplr();!!!!!// 110 .wr_coe_clk(clk),//初始化系数写入时钟 111 .wr_coe_en(csh0),//系数配置的使能端，由初始化模块地址译码产生，方便不同系数blk的选通 112 .start(start),//输入的启动卷积和fir的控制端 113 //.shft_out_dp_data(), 114 .shft_out_dp_data(shft_data1[15:0]), 115 //.s_latch(conv_res) 116 .s_latch(conv_res1[39:0]) 117 ); 118 CONV_SER16 i_conv_ser16_II(//例化第二个16阶卷积/fir电路 119 .clk(clk), 120 .a(shft_data1[15:0]),//由级联的第一节移出的数据 121 .en(rdy_work),//系数配置完成后才能使能 122 .coe_data_in16(data_bank),//初始化系数的数据输入端 123 .wr_coe_addr(~wr_blk_addr[3:0]),//初始化系数的地址输入端 124 ////!!!特别注意这里，乘加时数据从大地址进入乘加操作，因此地址要求补码，将系数翻转过来。相当于：fliplr();!!!!!// 125 .wr_coe_clk(clk),//初始化系数写入时钟 126 .wr_coe_en(csh1),//系数配置的使能端，由初始化模块地址译码产生，方便不同系数blk的选通 127 .start(start),//输入的启动卷积和fir的控制端 128 .shft_out_dp_data(shft_data2[15:0]), 129 .s_latch(conv_res2[39:0]) 130 ); 131 CONV_SER16 i_conv_ser16_III(//例化第3个16阶卷积/fir电路 132 .clk(clk), 133 .a(shft_data2[15:0]),//由级联的第一节移出的数据 134 .en(rdy_work),//系数配置完成后才能使能 135 .coe_data_in16(data_bank),//初始化系数的数据输入端 136 .wr_coe_addr(~wr_blk_addr[3:0]),//初始化系数的地址输入端 137 ////!!!特别注意这里，乘加时数据从大地址进入乘加操作，因此地址要求补码，将系数翻转过来。相当于：fliplr();!!!!!// 138 .wr_coe_clk(clk),//初始化系数写入时钟 139 .wr_coe_en(csh2),//系数配置的使能端，由初始化模块地址译码产生，方便不同系数blk的选通 140 .start(start),//输入的启动卷积和fir的控制端 141 .shft_out_dp_data(shft_data3[15:0]), 142 .s_latch(conv_res3[39:0]) 143 ); 144 CONV_SER16 i_conv_ser16_IV(//例化第4个16阶卷积/fir电路 145 .clk(clk), 146 .a(shft_data3[15:0]),//由级联的第一节移出的数据 147 .en(rdy_work),//系数配置完成后才能使能 148 .coe_data_in16(data_bank),//初始化系数的数据输入端 149 .wr_coe_addr(~wr_blk_addr[3:0]),//初始化系数的地址输入端 150 ////!!!特别注意这里，乘加时数据从大地址进入乘加操作，因此地址要求补码，将系数翻转过来。相当于：fliplr();!!!!!// 151 .wr_coe_clk(clk),//初始化系数写入时钟 152 .wr_coe_en(csh3),//系数配置的使能端，由初始化模块地址译码产生，方便不同系数blk的选通 153 .start(start),//输入的启动卷积和fir的控制端 154 .shft_out_dp_data(shft_data4[15:0]), 155 .s_latch(conv_res4[39:0]) 156 ); 157 PADD i_PADD (//为了方便DA输出将所有输出偏置为正数 158 .clock ( start ), 159 .data0x ( conv_res1[31:0] ), 160 .data1x ( conv_res2[31:0] ), 161 .data2x ( conv_res3[31:0] ), 162 .data3x ( conv_res4[31:0] ), 163 .result ( acc_sum[33:0] ) 164 ); 165 166 MCP4822 i_mcp4822(.clk100m(clk), 167 .rst_n(rdy_work), 168 .start(start), 169 .dac_data_a(data_cha),//通道a输出的数据 170 //.dac_data_b(dds_data12),//通道b连接DDS内容 171 //.dac_data_b(product12),//通道b连接乘法的高12位 172 //.dac_data_b(conv_data12),//通道b连接卷积的结果 173 .dac_data_b(data_chb),//通道b输出的数据 174 .cs(cs_da), 175 .mosi(mosi_da), 176 .sck(sck_da), 177 .ld(ld_da) 178 ); 179 MCP3202 i_mcp3202( 180 .clk100m(clk), 181 .rst_n(rdy_work), 182 .start(start), 183 .adc_data_a(ad_data), 184 .cs(cs_ad), 185 .mosi(mosi_ad), 186 .miso(miso_ad), 187 .sck(sck_ad) 188 ); 189 190 sub_offset i_sub_offset (//这个减法器用于从AD结果中去除直流偏置，结果是16位补码 191 .dataa ( {4'd0,acc_sum[26:15]} ),//卷积结果用于计算 192 //.dataa ( {4'd0,ad_data[11:0]} ),//AD结果是正整数，只要补零就可以得到16位补码 193 .datab ( 16'h0000 ),//直流偏置认为是1/2满幅度16'h0800，没有直流偏置则为16'h0000 194 .result ( ac_sig[15:0] )//去除直流偏置以后的结果，是16bit补码 195 ); 196 197 pow_cal i_pow_cal (//平方运算计算信号能量 198 .dataa ( ac_sig[15:0] ), 199 .result ( pow_sig[31:0] ) 200 ); 201 202 shft_reg i_shft_reg (//例化移位寄存器 203 .clock ( start ), 204 .shiftin ( pow_sig[31:0] ), 205 .shiftout ( shiftout[31:0] ), 206 .taps ( ) 207 ); 208 sub_add sub_add_inst (//减法器，用送入累加结果的数减去要从累加结果中拿出的数，结果有可能是负数 209 //减完后再送入累加器，以实现累加器内容的维护：最近64个数的和 210 .dataa ( pow_sig[31:0] ), 211 .datab ( shiftout[31:0]), 212 .result ( acc_in[31:0] ) 213 ); 214 ////////计算最近64个数的和///////// 215 always @ (posedge start or negedge rst_n) 216 begin 217 if(!rst_n) 218 acc_out[39:0] <= 40'd0; 219 else begin 220 acc_out[39:0] = $signed(acc_out[39:0]) + $signed(acc_in[31:0]); 221 end 222 end 223 ///////对累加和开方，提高小信号的分辨率 224 ip_sqrt ip_sqrt_inst ( 225 .radical ( acc_out[39:0] ), 226 .q ( data_out[19:0] ), 227 .remainder ( ) 228 ); 229 endmodule
上述顶层设计文件中大部分内容是各个功能模块电路的例化和连接语句，由于代码注释很详细，具体细节就不在这里赘述了，我主要介绍以下几点：
1、 其中只有唯一的一个always过程赋值语句，该语句维护了计数器cnt[15:0]的工作。总的同步启动信号start产生在cnt[15:0]为某些具体值的时候，因此cnt[15:0]技术的周期也决定了start信号产生的周期，也就是信号的周期。
2、正如我在本系列的准备篇“用Verilog-HDL状态机控制硬件接口”所介绍的，使用A/D和D/A进行FPGA算法调试的优势在于：可以将算法中各个中间节点上的信号通过D/A转换器呈现出来。代码中对data_cha[11:0]和data_chb[11:0]的赋值操作有大量的可选注释，正是由于这个原因预留的。
3、每个卷积节CONV_SER16的例化代码中，接口wr_coe_addr是系数存储器的初始化写入地址（关于系数存储器初始化部分的内容，请参见本系列博文第二篇“互相关/卷积/FIR电路的实现”中的描述）。但连接时wr_coe_addr却对系数初始化器的输出的地址wr_blk_addr取反，原因是每个卷积节在执行卷积时首先和输入进行乘加的是地址较大的系数存储器中的内容，因此地址应进行翻转。否则存储在系数初始化池中的完整系数序列，将在每个卷积节中发生翻转，从而无法产生正确结果。
四、测试结果
Quartus-II中针对Cyclone-I的FPGA综合的结果如下图所示。
图5 综合、布局和布线的结果
用matlab产生6正弦周期，长度为60个采样点的信号，在前后各补2个0，从而得到长度为64个采样点的正弦信号。之所以采用60个采样点包含6个正弦周期，是为了模拟5KHz的正弦信号@50KSPS的采样率。最后再为64个采样点加上blackman窗以降低泄露的影响。
图6 理想的目标反射回波信号
我将该信号下载到是德科技任意信号发生器33521中，用手动触发方式让33521产生上图所示的波形，以模拟理论回波信号。将其连接到FPGA平台的A/D输入端，将处理结果从FPGA平台的D/A输出。用示波器观测输入和处理结果如下图所示。
图7 实测处理结果
由上图可知FPGA中的算法实现了本系列博文第一篇中matlab理论仿真的输出波形，能够满足目标反射回波检测算法的要求。蓝色卷积结果的延时是由于FPGA算法的潜伏期造成的，不影响系统的实时性。
当然这里完成的只是一个算法验证性的实验，若要真正使用，还需要将A/D和D/A转换的速度进一步提高。
关于目标回波检测算法原理及实现步骤过程，请从本系列博文的开头开始浏览【目标反射回波检测算法及其FPGA实现 之一：算法概述】。