FIR數字濾波器的FPGA實現(二)-序列FIR濾波器設計(1)
(二)FIR數字濾波器的FPGA實現-序列FIR濾波器設計
文章目錄
- (二)FIR數字濾波器的FPGA實現-序列FIR濾波器設計
- 0 序列FIR濾波器基本原理
- 1 基於移位暫存器的序列 FIR 濾波器
- 2 基於雙埠 RAM 的序列 FIR 濾波器
- 3 係數對稱的序列 FIR 濾波器的設計
- 4 兩種序列結構的 FIR 濾波器效能比較
對於FIR濾波器主要涉及到濾波器的設計和濾波器的實現，設計和實現的區別如下圖所示：
實現是 forward problem，設計是 inverse problem
What are inverse problems?
本文主要涉及到FIR濾波器的實現，在實現的過程中，h[k]都是已知的，而h[k]的求解一般是設計的過程。具體內容包括，FIR濾波器的基本原理，序列FIR濾波器設計（此設計為濾波器實現的“設計”和FIR濾波器的設計不同，自行理會），並行FIR濾波器設計，串並FIR濾波器設計，分散式FIR濾波器設計，快速卷積型 FIR 濾波器、多通道 FIR 濾波器、多頻響 FIR 濾波器。對於快速卷積型 FIR 濾波器、多通道 FIR 濾波器、多頻響 FIR 濾波器會簡單介紹，其中序列、並行、串並、分散式FIR濾波器設計會給出相應原始碼和模擬模型，如果條件允許會抽出一個原始碼在FPGA上執行，並進行實驗分析。
0 序列FIR濾波器基本原理
對於基本原理，有不懂得或者忘記得，強烈建議去MOOC上看下數字訊號處理（北交大 陳後金老師）的視訊（註冊免費觀看），很經典，老師講解的非常不錯，還有相關習題可以練習。
總結下：
長度（抽頭數）為N、階數為 N-1 的 FIR 系統的轉移函式、差分方程和單位衝激響應分別如式（2.1）、式(2.2)和式(2.3)所示。
H和Xn分別稱為係數向量和輸入向量。則式（2.2）可進一步表示為
從以上表達式可以看出，H(z)實際上是z-1的多項式，它的 N-1 個極點全部位於 z = 0處，所以一個 FIR 濾波器始終是穩定的。它的零點分佈可以處於有限 Z 平面的任何位置，當全部零點位於單位圓內時，就成為最小相位系統。僅有零點存在也正是 FIR 濾波器被稱作全零點濾波器的原因。
定義一次濾波運算：由當前輸入和過去輸入計算當前輸出。式 (2.2) 表明完成一次濾波運算需要 N 次乘法和 N-1 次加法。
FIR 濾波器的兩種結構：直接型和轉置型，兩者都是從時域角度描述的。
根據上一節介紹，畫出四抽頭FIR濾波器直接型結構。
基於式 (2.2) 的卷積形式，FIR 濾波器資料流的動態變化過程如圖 2.2 所示。這裡以 4抽頭為例說明。
根據圖 2.2，可將 FIR 濾波器視為一個數據窗，只對位於窗內的資料進行乘加操作以完成濾波運算。序列 FIR 濾波器的設計思想是隻用一個乘法器和一個累加器按時間順序依次完成一次濾波運算所需的N1次乘法和 N-1 次加法。設計的關鍵環節是延時線，這可通過移位暫存器或者單埠 RAM 來實現。
1 基於移位暫存器的序列 FIR 濾波器
1.1 基本理論
基於移位暫存器的序列 FIR 濾波器其設計思想是輸入資料的動態流動是以自身的流動而非通過地址的改變實現的。對圖 2.2 進一步分析可知，每完成一次濾波運算，輸入資料沿資料窗推進一格即總有一個新的資料進入資料窗並佔據資料窗的第一個位置，資料窗內的原有資料整體向前推移一個係數位，且總有一個數據移出資料窗。以fs表示輸入資料取樣率，它反映了輸入資料進入資料窗的速率；以fclk 表示濾波運算的工作速率，它反映了乘法器的工作速率。定義輸入資料的生命週期是指輸入資料從進入資料窗至移出資料窗所持續的時間。例如，圖2.2所示的輸入資料x(0)的生命週期為4 Tclk(其中Tclk=1/fclk)。對於N抽頭FIR濾波器，一次濾波運算需要N次乘法運算，故fs和fclk之間滿足下式所示的關係式。
輸入資料的這種流動特性可採用如圖 2.3 所示的移位暫存器實現。在 Xilinx 7 系列FPGA 中，該模型可用查詢表實現，可呼叫原語( Primitive) SRL16E。此外，也可呼叫 LogiCORE IP RAM-Based Shift Register 該結構具有很強的擴充套件性，這裡僅以深度為 4 舉例說明。
圖 2.3 中， ce 作為寫資料使能端，當其為高時將前級輸出資料寫入下一級。顯然 ce 的週期決定了資料的持續時間，反映了資料的生命週期。 din 作為輸入資料端，所有的輸入資料由此依時間順序進入移位暫存器。 addr 作為資料選擇器 MUX 的控制端，決定了輸出資料來自於哪個暫存器，顯然 addr 的變化速率決定了輸出資料率。 dout 是移位暫存器的輸出資料端。 可見， 該結構具有 “FIFO 進 RAM 出” 的特點即資料總是由初級暫存器進入， 而輸出取決於地址控制端。 需要注意的是， 最後一級延時單元通常採用暫存器實現， 其餘延時單元採用查詢表實現。 由於查詢表本身並沒有復位埠， 因此， 若該結構帶復位功能， 那麼只能復位末級暫存器。
基於移位暫存器的序列 FIR 濾波器硬體結構如圖 2.4 所示。 整個系統由移位暫存器、係數 ROM 控制單元和乘加器構成。 實際應用時， 為了流水處理通常在移位暫存器 MUX的輸出端增加一級寄存。係數存放在 ROM 中，可根據實際情況選擇分散式 RAM( Distributed RAM ) 或 Block RAM 圖中虛線框部分可用 Xilinx 7 系列 FPGA 中的 DSP48E1實現，對其進行動態配置完成乘加與乘累加運算之間的切換， 這可由圖中的 load 訊號控制。也可以將乘法運算與累加運算分別採用獨立模組實現。 濾波運算最終結果 firout 由末級捕獲暫存器的使能訊號（ 也稱為捕獲訊號） 使能輸出。
圖 2.4 所示各訊號之間的時序關係如圖 2.5 所示。
此處以 4 抽頭 FIR 濾波器為例。顯然， addr 是一個模值為 4 的計數器， 由 sclr 清零， 既是 MUX 的控制訊號又是係數 ROM 的讀地址訊號。 ce 由 addr 譯碼產生， 延時 3 個時鐘週期作為 capture 訊號。 這緣於訊號由乘加器輸入端 A 至輸出端 P 共需要 3 個時鐘週期。 由於此設計中 ce 週期為 4Tclk，將 其 延 時 3 個時鐘週期將與原始訊號重合， 故 capture 即為 ce。SRL16E 本身並沒有復位埠， 這裡的 sclr是 SRL16E 輸出埠所新增流水暫存器的復位訊號。
在設計時還要考慮加法器的位寬， 以保證累加運算不溢位。 假定輸入資料位寬為B， 濾波器係數位寬為 C，抽頭數為 N, 那麼輸入資料與係數相乘的結果位寬為 B+C， 累加 H次，故輸出資料位寬為
其中， ⌈ ⌉表示向上取整。式（2.8) 儘管可以保證不溢位，但卻造成了硬體資源的浪費。從另一個角度考慮，輸出資料的最大值為
只要保證位寬滿足ymax 的需求，即可確定加法器的位寬為：
結論:
- (1) 採用移位暫存器設計時， 移位暫存器的深度與濾波器抽頭數一致;
- (2) 輸入資料取樣率與系統時鐘頻率滿足式 (2.10) 所示關係式。
1.2 設計實現
FIR濾波器設計方法是以直接逼近所需離散時間系統的頻率響應為基礎。設計方法包括窗函式法和最優化方法(等同紋波法)，其中窗函式方法是設計FIR數字濾波器是最常用的方法之一。
在時域中，FIR濾波器的輸入/輸出就是一個輸入訊號與單位脈衝相應的卷積。離散方程為y(n)=x(n)*h(n)=∑x(k)h(n-k)=∑h(k)x(n-k)，其中y(n)為濾波輸出，x(n)為取樣資料，h(n)為濾波抽頭係數.設計FIR濾波器就是要找到N個係數。N-1階濾波器通常需要N個係數描述，通常需要N個乘法器和N-1個2輸入加法器實現。根據FIR表示式，濾波器實質上就是進行乘累加運算，乘累加的次數由濾波器階數決定。其序列結構如圖
係數表的產生我們利用matlab的FDATool來完成，也可以自己進行引數，利用FDATool比較方便。
詳細使用教程: https://blog.csdn.net/Pieces_thinking/article/details/83656785
首先開啟matlab的FDATool工具，選擇一個8階的低通濾波，通頻率為500Hz，截止頻率為750Hz。取樣頻率2000Hz。
FDATool介面左下側排列了一組工具按鈕，其功能分別如下所述：
- 濾波器轉換（TransForm Filter）
- 設定量化引數（Set Quantization Parameters）
- 實現模型（Realize Model）
- 匯入濾波器（Import Filter）
- 多速率濾波器（Multirate Filter）
- 零極點編輯器（Pole-zero Editor）
- 設計濾波器（Design Filter）
選擇其中的選擇Design Filter按鈕，進入設計濾波器介面，進行下列選擇，如圖所示。濾波器型別（Filer Type）為低通（Low Pass）
設計方法（Design Method）為FIR/IIR，分別採用Equiripple、Least-squares、Window、Constr.LeastPth-norm 、Constrained Equiripple、Constr.Band Equiripple（FIR濾波器設計）和 Butterworth、Chebyshev Type I、Chebyshev Type I、Elliptic、Maximally flat、Least Pth-norm、Constr.LeastPth-norm（IIR濾波器設計）。
在我們設定好FDATool的引數後，我們可以點選Analysis–>Filter Coefficients來觀察FIR濾波器的係數。
這裡的係數全是有符號型的小數，我們在FPGA中需要用整數作為濾波器的係數，所以我們要進行係數的歸一化，點選左下角設定量化引數(Set Quantization Parameters)，Filter arithmetic選擇Fixed-point(定點)。然後就可以匯出Xilink的.coe檔案了。Targets–>XILINK Coefficients(.COE) File.儲存後matlab自動開啟該.coe檔案。
在Matlab中編寫以下檔案並執行：
function [y,yref] = Serial_Shift_registers_FIR(x,h)
% x: input data vector
% h: fir coef [h(0),h(l),....h(N-l)]
x_len = length(x);
h_len = length(h);
y = zeros(x_len,1);
tap_delay = zeros(h_len,1);
for k=1:x_len
tap_delay = [x(k);tap_delay(1:h_len-1)];
for n=1:h_len
y(k) = y(k)+tap_delay(n)*h(n);
end
end
yref = filter(h,1,x);
輸入以下語句執行
Serial_Shift_registers_FIR([10,11,12,13,14,15,16,17,18,19,20],[252,242,37,79,37,242,252,7,252])
結果：
針對以上分析，編寫Verilog檔案如下：
`timescale 1ns / 1ps
//////////////////////////////////////////////////////////////////////////////////
// Company:
// Engineer:
//
// Create Date: 2018/11/19 09:11:33
// Design Name:
// Module Name: Serial_Shift_registers_FIR
// Project Name:
// Target Devices:
// Tool Versions:
// Description:
//
// Dependencies:
//
// Revision:
// Revision 0.01 - File Created
// Additional Comments:
//
//////////////////////////////////////////////////////////////////////////////////
module Serial_Shift_registers_FIR(
Clk,Rst_n,data_in,data_out
);
parameter word_data_in = 8;
parameter word_data_out = 2*word_data_in+1;
parameter order = 8;
parameter c0 = 8'hfc;
parameter c1 = 8'hf2;
parameter c2 = 8'h25;
parameter c3 = 8'h4f;
parameter c4 = 8'h25;
parameter c5 = 8'hf2;
parameter c6 = 8'hfc;
parameter c7 = 8'h07;
parameter c8 = 8'hfc;
input Clk;
input Rst_n;
input [word_data_in-1:0]data_in;
output wire [word_data_out-1:0]data_out;
reg [word_data_in-1:0]samples[order-1:0];
reg [word_data_in-1:0]data_in_tmp1;
reg [word_data_in-1:0]data_in_tmp2;
// wire [word_data_out-1:0]data_out=0;
integer i;
[email protected](posedge Clk)
if(!Rst_n)begin
data_in_tmp1<=8'd0;
data_in_tmp2<=8'd0;
end
else begin
data_in_tmp1<=data_in;
data_in_tmp2<=data_in_tmp1;
end
[email protected](posedge Clk)//同步復位,完成移位功能
if(!Rst_n)begin
for(i=0;i<order;i=i+1)
samples[i]<=0;
end
else begin
samples[0]<=data_in_tmp2;
for(i=1;i<order;i=i+1)
samples[i]<=samples[i-1];
end
//序列實現累加
assign data_out = c0*data_in_tmp2 + c1*samples[0] + c2*samples[1] + c3*samples[2] + c4*samples[3] + c5*samples[4] + c6*samples[5] + c7*samples[6] + c8*samples[7];
endmodule
和matlab檔案類似，只不過將濾波器係數直接新增到檔案中了，可以根據需求自己修改濾波器係數，這種結構的寫法，很方便很好理解，就是移位相乘、相加，編寫TestBench，如下：
`timescale 1ps / 1ps
//////////////////////////////////////////////////////////////////////////////////
// Company:
// Engineer:
//
// Create Date: 2018/11/19 09:51:23
// Design Name:
// Module Name: SSRFIR_tb
// Project Name:
// Target Devices:
// Tool Versions:
// Description:
//
// Dependencies:
//
// Revision:
// Revision 0.01 - File Created
// Additional Comments:
//
//////////////////////////////////////////////////////////////////////////////////
module SSRFIR_tb(
);
parameter word_data_in = 8;
parameter word_data_out = 2*word_data_in+1;
reg Clk;
reg Rst_n;
reg [word_data_in-1:0]data_in;
wire [word_data_out-1:0]data_out;
Serial_Shift_registers_FIR u1
(
.Clk(Clk),
.Rst_n(Rst_n),
.data_in(data_in),
.data_out(data_out)
); // Results
reg [15:0] cnt;
reg signed [15:0] e [0:10]; // Coefficient array
initial
begin
#0 Clk = 1'b0;
data_in <=0;
e[0]=10;
e[1]=11;
e[2]=12;
e[3]=13;
e[4]=14;
e[5]=15;
e[6]=16;
e[7]=17;
e[8]=18;
e[9]=19;
e[10]=20;
#10 Rst_n = 1'b0;
#20 Rst_n = 1'b1;
cnt <=0;
#1000000 $stop;
end
always #10
begin
Clk = ~Clk;
end
always @ (posedge Clk)
begin
if(!Rst_n)
cnt <= 1'b0;
else
begin
if(cnt == 16'd10000-1)
cnt = 1'b0;
else
cnt = cnt + 1'b1;
end
end
always @ (posedge Clk)
begin
case (cnt)
20:
data_in <= 0;
21:
data_in <=e[0];
22:
data_in <=e[1];
23:
data_in <=e[2];
24:
data_in <=e[3];
25:
data_in <=e[4];
26:
data_in <=e[5];
27:
data_in <=e[6];
28:
data_in <=e[7];
29:
data_in <=e[8];
30:
data_in <=e[9];
31:
data_in <=e[10];
32:
data_in <=1'b0;
default:
data_in<=data_in;
endcase
end
endmodule
利用Vivado和Modsim模擬，結果如下圖，和MATLAB的結果一樣。
以上設計方法，主要是為了直觀理解基於移位暫存器的序列 FIR 濾波器設計，實際使用過程中由於浪費的資源過多，很少使用，下面介紹一個優化後的濾波器設計（參考《數字濾波器的MATLAB與FPGA實現》），如下：
`timescale 1ns / 1ps
//----------------------------------------------------------
// 全序列FIR濾波器
//
//----------------------------------------------------------
module Serial_Shift_registers_FIR
(
input rst,
input clk, //系統時鐘16kHz
input signed [11:0] Xin, //輸入資料頻率2kHz
output signed [28:0] Yout
);
//----------------------------------------------------------
// 例項化乘法器和加法器IP核
//----------------------------------------------------------
reg signed [12:0] add_a,add_b;
wire signed [12:0] add_s;
c_addsub_1 add
(
.A (add_a),
.B (add_b),
.S (add_s)
);
reg signed [11:0] coe; //12bit量化濾波器係數
wire signed [24:0] Mout;
mult_gen_1 mult
(
.CLK (clk),
.A (coe),
.B (add_s),
.P (Mout)
);
//----------------------------------------------------------
// 模8計數器，控制輸入資料移位和濾波結果輸出
//----------------------------------------------------------
reg [2:0] cnt;
reg [11:0] Xin_Reg[15:0];
reg [3:0] i,j;
always @ (posedge clk or posedge rst)
if (rst) cnt <= 'd0;
else cnt <= cnt + 1'b1;
always @ (posedge clk or posedge rst)
if (rst) begin //清0
Xin_Reg[15] <= 'd0;
for (i=0; i<15; i=i+1'b1)
Xin_Reg[i] <= 'd0;
end
else begin
if (cnt == 'd7) begin
for (j=0; j<15; j=j+1'b1) //移位
Xin_Reg[j+1] <= Xin_Reg[j];
Xin_Reg[0] <= Xin;
end
end
//----------------------------------------------------------
// 計算
//----------------------------------------------------------
always @ (posedge clk or posedge rst)
if (rst) begin
add_a <= 'd0; add_b <= 'd0; coe <= 'd0;
end
else begin
if (cnt == 'd0) begin //第一個週期
add_a <= {Xin_Reg[0][11],Xin_Reg[0]};
add_b <= {Xin_Reg[15][11],Xin_Reg[15]};
coe <= 12'h000;
end
else if (cnt == 'd1) begin //第二個週期
add_a <= {Xin_Reg[1][11],Xin_Reg[1]};
add_b <= {Xin_Reg[14][11],Xin_Reg[14]};
coe <= 12'hffd;
end
else if (cnt == 'd2) begin //第三個週期
add_a <= {Xin_Reg[2][11],Xin_Reg[2]};
add_b <= {Xin_Reg[13][11],Xin_Reg[13]};
coe <= 12'h00f;
end
else if (cnt == 'd3) begin //第四個週期
add_a <= {Xin_Reg[3][11],Xin_Reg[3]};
add_b <= {Xin_Reg[12][11],Xin_Reg[12]};
coe <= 12'h02e;
end
else if (cnt == 'd4) begin //第五個週期
add_a <= {Xin_Reg[4][11],Xin_Reg[4]};
add_b <= {Xin_Reg[11][11],Xin_Reg[11]};
coe <= 12'hf8b;
end
else if (cnt == 'd5) begin //第六個週期
add_a <= {Xin_Reg[5][11],Xin_Reg[5]};
add_b <= {Xin_Reg[10][11],Xin_Reg[10]};
coe <= 12'hef9;
end
else if (cnt == 'd6) begin //第七個週期
add_a <= {Xin_Reg[6][11],Xin_Reg[6]};
add_b <= {Xin_Reg[9][11],Xin_Reg[9]};
coe <= 12'h24e;
end
else begin //第八個週期
add_a <= {Xin_Reg[7][11],Xin_Reg[7]};
add_b <= {Xin_Reg[8][11],Xin_Reg[8]};
coe <= 12'h7ff;
end
end
//----------------------------------------------------------
// 對乘法結果進行累加
//----------------------------------------------------------
reg signed [28:0] sum;
reg signed [28:0] yout;
always @ (posedge clk or posedge rst)
if (rst) begin
sum <= 'd0; yout <= 'd0;
end
else begin
if (cnt == 'd2) begin
yout <= sum; sum <= 'd0; sum <= sum+Mout;
end
else sum <= sum+Mout;
end
assign Yout = yout;
endmodule
使用MATLAB設計一個2kHz取樣，500Hz截止的15階低通濾波器（h(n)長度為16），量化位數為12bit，輸入訊號位寬也為12bit。
系統時鐘與輸入資料頻率（即取樣頻率）相同即可；但是序列結構由於需要n/2倍（此處為8倍）個週期完成一次運算，則FPGA系統時鐘需要是取樣頻率的n/2倍。
例項化全序列FIR濾波器結構所需的乘法器和加法器IP核需要自己新增。
使用MATLAB生成一個200khz+800kHz的混合頻率訊號，寫入txt檔案，。編寫Testbench讀取txt檔案對訊號濾波。
Xilinx_SerialFIR_tb.v檔案
`timescale 1 ns/ 1 ns //設定模擬時間單位:ns
module Xilinx_SerialFIR_tb();
reg [11:0] Xin;
reg clk;
reg rst;
reg clk_data; //資料時鐘，速率為系統時鐘clk的1/8；
// wires
wire [28:0] Yout;
Serial_Shift_registers_FIR u1 (
.Xin(Xin),
.Yout(Yout),
.clk(clk),
.rst(rst)
);
parameter clk_period=625; //設定時鐘訊號週期（頻率）:1.6MHz
parameter clk_period_data=clk_period*8;
parameter clk_half_period=clk_period/2;
parameter clk_half_period_data=clk_half_period*8;
parameter data_num=2000; //模擬資料長度
parameter time_sim=1250000; //模擬時間data_num*clk_period
initial
begin
//設定時鐘訊號初值
clk=1;
clk_data=1;
//設定復位訊號
rst=1;
#100 rst=0;
//設定模擬時間
#12500000 $finish;
//設定輸入訊號初值
Xin=12'd10;
end
//產生時鐘訊號
always
#clk_half_period clk=~clk;
always
#clk_half_period_data clk_data=~clk_data;
//從外部TX檔案(SinIn.txt)讀入資料作為測試激勵
integer Pattern;
reg [11:0] stimulus[1:data_num];
initial
begin
$readmemb("G:/USE/FILE/DSP/FPGA/Serial_Shift_registers_FIR_B/sin.txt",stimulus);
Pattern=0;
repeat(data_num)
begin
Pattern=Pattern+1;
Xin=stimulus[Pattern];
#clk_period_data;//資料週期為時鐘週期的8倍
end
end
endmodule
發現Yout的幅度變化很低，檢測後發現24-29bit都是符號位，將0-24bit提取為一個新的虛擬匯流排Yout_bus。明顯看到經過500Hz低通濾波器濾波後，輸入的200+800Hz訊號只剩下200Hz的單頻訊號。
