目標反射回波檢測算法及其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實現 之一：算法概述】。