在FPGA上实现一个最简单的FIR滤波器电路
一、概要
在数字滤波器中,FIR滤波器有全局稳定性和线性相位特性,且结构简单,容易在FPGA上实现。在音频、视频、雷达、通信等信号处理领域应用广泛。
本文按照:FIR滤波器简介——基于MATLAB的FIR滤波器设计——FPGA内部的FIR滤波器部署 的流程讲解。
二、FIR滤波器简介
有限长单位冲激响应滤波器(Finite Impulse Response)
非递归性:当前时刻的输出只取决于输入信号,与之前的输出信号无关。
结构简单,硬件实现更容易
线性相位:保证幅度特性的同时,很容易做到严格的线性相位特性。
在视频信号中对相位特性要求非常严格,在音频信号中则对幅度要求严格
高阶数高延迟:与IIR滤波器相比,要实现相同性能,FIR滤波器需要更高的阶数,这就导致了更大的延迟。FIR滤波器的实时性较IIR滤波器差。
三、基于MATLAB的FIR滤波器设计
FIR滤波器设计时根据具体需要有很多指标,除了最基本的截止频率,还有通带衰减、阻带衰减、过度带宽、群延迟等。可以根据自己的指标需要来设计滤波器,但在这里为了简单起见,我们只考虑截止频率。
1,低通、高通滤波器
Fs = 5000; % 采样频率(Hz)
N = 63; % 滤波器阶数(通常为奇数)
Fc = 200; % 滤波器截止频率(Hz)
b = fir1(N, Fc / (Fs / 2), 'low' ); % 低通滤波器的FIR参数
b = fir1(N, Fc / (Fs / 2), 'high'); % 高通滤波器的FIR参数
2,带通滤波器
Fs = 5000; % 采样频率(Hz)
F1 = 400; % 滤波器下边界频率(Hz)
F2 = 600; % 滤波器上边界频率(Hz)
f = [0 F1 F2 Fs/2] / (Fs/2); % 归一化频率
A = [0 1 1 0]; % 带通滤波器的幅度响应(1表示通带,0表示阻带)
b_bandpass = firpm(N, f, A);
A = [1 0 0 1]; % 带阻滤波器的幅度响应(1表示通带,0表示阻带)
b_bandstop = firpm(N, f, A);
我们需要的就是最后生成的
b
b
b(FIR滤波器参数)。
注意:如果滤波器阶数为63,那么 b 的长度就是 63+1=64。
FIR滤波器时域公式:
y ( n ) = ∑ k = 0 M − 1 h ( k ) x ( n − k ) y\left(n\right)=\sum_{k=0}^{M-1}h\left(k\right)x\left(n-k\right) y(n)=∑k=0M−1h(k)x(n−k)
简单概括,就是输入信号和 b b b 作卷积运算,得到的结果就是输出信号。
举例说明:
假如FIR滤波器阶数为3,由FIR时域公式得到
y ( n ) = b ( 0 ) x ( n ) + b ( 1 ) x ( n − 1 ) + b ( 2 ) x ( n − 2 ) + b ( 3 ) x ( n − 3 ) y(n)=b(0)x(n)+b(1)x(n-1)+b(2)x(n-2)+b(3)x(n-3) y(n)=b(0)x(n)+b(1)x(n−1)+b(2)x(n−2)+b(3)x(n−3)
我们可以看出,一个FIR滤波器需要存储
n
n
n 个参数和
n
n
n 个数据,作
n
n
n 次乘法和
n
−
1
n-1
n−1 次加法。
由数字信号处理知识可知,实信号FIR滤波器的 b 参数,是对称的。
也就是说
b ( k ) = b ( N − k ) b(k)=b(N-k) b(k)=b(N−k) ( N N N是滤波器阶数)
利用滤波器抽头的对称性质,我们可以节省一半的乘法计算。
还是以3阶FIR滤波器为例
y ( n ) = b ( 0 ) x ( n ) + b ( 1 ) x ( n − 1 ) + b ( 2 ) x ( n − 2 ) + b ( 3 ) x ( n − 3 ) y(n)=b(0)x(n)+b(1)x(n-1)+b(2)x(n-2)+b(3)x(n-3) y(n)=b(0)x(n)+b(1)x(n−1)+b(2)x(n−2)+b(3)x(n−3)
= b ( 0 ) [ x ( n ) + x ( n − 3 ) ] + b ( 1 ) [ x ( n − 1 ) + x ( n − 2 ) ] =b(0)[x(n)+x(n-3)]+b(1)[x(n-1)+x(n-2)] =b(0)[x(n)+x(n−3)]+b(1)[x(n−1)+x(n−2)]
在实现FIR滤波器时,我们可以先将 n n n 个输入数据对位相加,之后再进行 n / 2 n/2 n/2 次乘法,最后对乘法结果相加得到输出数据。
四、FPGA内部的FIR滤波器部署
举例 实现一个63阶的FIR滤波器
说明:由于MATLAB导出的滤波器参数b的值特别小(比如63阶滤波器,b参数有64个,那么这64个参数的和为1)所以在FPGA中,如果不想引入浮点数运算单元,就要对滤波器参数b等比扩大,比如乘以2^16,直到所有参数都是整数或者四舍五入为整数。等到计算出最终结果后,再除以2^16,保证结果的值大小不改变。
举例:b参数为 [0.25 0.25 0.25 0.25] 同时乘以 2^2 ,得到 [1 1 1 1] ,在FPGA内计算出结果为 137 ,需要除以 2^2 ,结果为 34 。
1、模块输入输出
module fir_filter#(
parameter DATA_WIDTH = 8 ,//输入数据位宽
parameter FILTER_M = 64 //滤波器系数+1
)(
input clk ,//时钟信号
input rst_n ,//低电平有效复位信号
input i_en ,//输入使能信号 在clk时钟下拉高一周期
output reg o_finish ,//该信号拉高一个周期表示滤波器输出值更新
input signed [DATA_WIDTH-1: 0] i_data_x ,//i_en持续时间内i_data_x会被存入移位寄存器
output reg signed [DATA_WIDTH-1: 0] o_data_y //o_finish拉高后o_data_y的值更新
);
2、使能打拍部分
//模块使能打拍
reg [ 7: 0] r_en ;
always @(posedge clk or negedge rst_n) begin
if (!rst_n) begin
r_en[7:0] <= 'b0 ;
end else begin
r_en[7:0] <= {r_en[6:0], i_en} ;
end
end
相当于把输入使能信号i_en打拍了8拍,这样就可以任意锁定i_en信号到来后的任何一个时钟周期。
比如我想在 第一个周期处理乘法,第二个周期处理一部分加法,第三个周期处理位宽截取,那就可以分别检测 r_en[0]、r_en[1]、r_en[2] 信号。个人认为很好用。
3、移位寄存器部分
//W2为b参数扩大的2的幂次 比如扩大了2^10,M2就是10
//输入数据位宽W1 对位相加位宽 W1 + 1 乘法器相乘后位宽 W1 + 1 + W2 累加器位宽 W1 + 1 + W2 + log2(M)
// 8 9 17 23
reg signed [DATA_WIDTH-1: 0] r_data_x[0:FILTER_M-1] ;
integer i,j ;
//移位寄存器
always @(posedge clk or negedge rst_n) begin
if (!rst_n) begin
for (i=0; i<FILTER_M; i=i+1) begin
r_data_x[i] <= 'b0;
end
end else if (i_en) begin
r_data_x[0] <= i_data_x ;
for (j=0; j<FILTER_M-1; j=j+1) begin
r_data_x[j+1] <= r_data_x[j] ; //周期性移位操作
end
end
end
这部分没有什么要强调的。
4、对位相加部分
//FIR对称特性 对位相加 节省一半计算量
reg signed [DATA_WIDTH+1-1: 0] r_add_reg [0: (FILTER_M/2)-1] ;
always @(posedge clk or negedge rst_n) begin
if (!rst_n) begin
for (i=0; i<FILTER_M/2; i=i+1) begin
r_add_reg[i] <= 'd0 ;
end
end else if (r_en[0]) begin
for (i=0; i<FILTER_M/2; i=i+1) begin
r_add_reg[i] <= r_data_x[i] + r_data_x[FILTER_M-1-i] ;
end
end
end
我选择在i_en后的第一个周期对移位寄存器中的数据进行对位相加。
5、乘法部分
//乘法模块
reg signed [2*DATA_WIDTH: 0] r_mout [0: (FILTER_M/2)-1] ;
always @(posedge clk or negedge rst_n) begin
if (!rst_n) begin
for (i=0; i<FILTER_M/2; i=i+1) begin
r_mout[i] <= 'd0 ;
end
end else if (r_en[1]) begin
r_mout[ 0] <= r_add_reg[ 0] * 8'd1 ;
r_mout[ 1] <= r_add_reg[ 1] * 8'd1 ;
r_mout[ 2] <= r_add_reg[ 2] * 8'd1 ;
r_mout[ 3] <= r_add_reg[ 3] * 8'd1 ;
r_mout[ 4] <= r_add_reg[ 4] * 8'd2 ;
r_mout[ 5] <= r_add_reg[ 5] * 8'd2 ;
r_mout[ 6] <= r_add_reg[ 6] * 8'd3 ;
r_mout[ 7] <= r_add_reg[ 7] * 8'd4 ;
r_mout[ 8] <= r_add_reg[ 8] * 8'd5 ;
r_mout[ 9] <= r_add_reg[ 9] * 8'd6 ;
r_mout[10] <= r_add_reg[10] * 8'd7 ;
r_mout[11] <= r_add_reg[11] * 8'd8 ;
r_mout[12] <= r_add_reg[12] * 8'd10 ;
r_mout[13] <= r_add_reg[13] * 8'd11 ;
r_mout[14] <= r_add_reg[14] * 8'd13 ;
r_mout[15] <= r_add_reg[15] * 8'd15 ;
r_mout[16] <= r_add_reg[16] * 8'd17 ;
r_mout[17] <= r_add_reg[17] * 8'd19 ;
r_mout[18] <= r_add_reg[18] * 8'd21 ;
r_mout[19] <= r_add_reg[19] * 8'd23 ;
r_mout[20] <= r_add_reg[20] * 8'd25 ;
r_mout[21] <= r_add_reg[21] * 8'd27 ;
r_mout[22] <= r_add_reg[22] * 8'd29 ;
r_mout[23] <= r_add_reg[23] * 8'd31 ;
r_mout[24] <= r_add_reg[24] * 8'd32 ;
r_mout[25] <= r_add_reg[25] * 8'd34 ;
r_mout[26] <= r_add_reg[26] * 8'd35 ;
r_mout[27] <= r_add_reg[27] * 8'd37 ;
r_mout[28] <= r_add_reg[28] * 8'd37 ;
r_mout[29] <= r_add_reg[29] * 8'd38 ;
r_mout[30] <= r_add_reg[30] * 8'd39 ;
r_mout[31] <= r_add_reg[31] * 8'd39 ;
end
end
可以看到,对于63阶滤波器,我只做了32次乘法运算。
并且注意,现在的b参数之和为1024,也就是2的10次方。所以最后计算出来的结果要右移10次。
6、累加器部分
//累加器部分
reg signed [ 2*DATA_WIDTH + 2: 0] r_sum1 ;
reg signed [ 2*DATA_WIDTH + 2: 0] r_sum2 ;
reg signed [ 2*DATA_WIDTH + 2: 0] r_sum3 ;
reg signed [ 2*DATA_WIDTH + 2: 0] r_sum4 ;
reg signed [ 2*DATA_WIDTH + 2: 0] r_sum5 ;
reg signed [ 2*DATA_WIDTH + 2: 0] r_sum6 ;
reg signed [ 2*DATA_WIDTH + 2: 0] r_sum7 ;
reg signed [ 2*DATA_WIDTH + 2: 0] r_sum8 ;
reg signed [ 2*DATA_WIDTH + 4: 0] r_sum_out ;
//累加器部分
always @(posedge clk or negedge rst_n) begin
if (!rst_n) begin
r_sum1 <= 'd0;
r_sum2 <= 'd0;
r_sum3 <= 'd0;
r_sum4 <= 'd0;
r_sum5 <= 'd0;
r_sum6 <= 'd0;
r_sum7 <= 'd0;
r_sum8 <= 'd0;
r_sum_out <= 'd0;
end else if(r_en[4]) begin
r_sum1 <= (r_mout[ 0] + r_mout[ 1]) + (r_mout[ 2] + r_mout[ 3]);
r_sum2 <= (r_mout[ 4] + r_mout[ 5]) + (r_mout[ 6] + r_mout[ 7]);
r_sum3 <= (r_mout[ 8] + r_mout[ 9]) + (r_mout[10] + r_mout[11]);
r_sum4 <= (r_mout[12] + r_mout[13]) + (r_mout[14] + r_mout[15]);
r_sum5 <= (r_mout[16] + r_mout[17]) + (r_mout[18] + r_mout[19]);
r_sum6 <= (r_mout[20] + r_mout[21]) + (r_mout[22] + r_mout[23]);
r_sum7 <= (r_mout[24] + r_mout[25]) + (r_mout[26] + r_mout[27]);
r_sum8 <= (r_mout[28] + r_mout[29]) + (r_mout[30] + r_mout[31]);
end else if(r_en[5]) begin
r_sum_out <= (r_sum1 + r_sum2) + (r_sum3 + r_sum4) + (r_sum5 + r_sum6) + (r_sum7 + r_sum8);
end
end
这里一共32个数据,我分成两个周期把他们加在一起了,分别是r_en[4]和r_en[5]。
7、位宽截取并输出
//位宽截取并输出
always @(posedge clk or negedge rst_n) begin
if (!rst_n) begin
o_data_y <= 'd0;
o_finish <= 'd0;
end else if(r_en[6]) begin
o_data_y <= (r_sum_out>>10)[DATA_WIDTH-1:0];
o_finish <= 'd1;
end else begin
o_data_y <= o_data_y;
o_finish <= 'd0;
end
end
可以看到,我把累加器结果右移了10位,最后截取了低DATA_WIDTH位输出。输入数据和输出数据位宽一样,当然设计的时候也可以不一样。
在r_en[6]时钟周期,拉高了一周期的o_finish,并将o_data_y结果更新。
五、注意事项
1、Verilog位宽截取
在Verilog中,默认是无符号数,有符号数用signed来声明。
有符号数的存储和读写方式是补码形式。在补码形式中,只要数值大小不超过存储范围,都可以用一个signed变量来直接截取。
2、
本Verilog程序只是给出一个讲解示例,在使用时许多参数和位宽需要根据自己的需要来计算和修改。