Verilog设计IIR低通滤波器

Verilog设计IIR低通滤波器

盘他!!!

设计一个4阶IIR低通滤波器,截止频率为2MHz, 滤波器采样频率为8Mhz,阻带衰减为40dB。量化系数为12bit。

image

IIR滤波器

IIR(无限冲激响应)滤波器,其单位脉冲冲击响应是无限长的。

差分方程为

image

可以看出,IIR存在输出对输入的反馈。

大多数IIR滤波器系数 a0 是1,传递函数为。

image

将一个4阶的直接II型结构的IIR滤波器,变成两个2阶IIR滤波器级联。

image

根据级联结构可以直接写出滤波器的差分方程。

image
$$
a_{01}y_1=b_{01}x(n)+b_{11}x(n-1)+b_{21}x(n-2) - (a_{11}y_1(n-1)+a_{21}y_1(n-2))\
a_{02}y_2=b_{02}y_1(n)+b_{12}y_1(n-1)+b_{22}y_1(n-2) - (a_{12}y_2(n-1)+a_{22}y_2(n-2))
$$
写出了差分方程,接下来需要求出对应的滤波系数。

求级联型IIR的滤波器系数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
fs=8*10^6;           %采样频率
fc=2*10^6; %低通滤波器的截止频率
Rs=40; %阻带衰减
qm=12; %滤波器系数量化位数
N=4; %滤波器阶数

%采用cheby2设计IIR低通滤波器
[b,a]=cheby2(N,Rs,fc*2/fs);
%绘制IIR低通滤波器的频率响应
freqz(b,a,128,fs);

%采用自编的E4_7_dir2cas函数完成滤波器的级联结构转换
[b0,B,A]=E4_7_dir2cas(b,a);
%将滤波器增益b0分配至第一级滤波器
B(1,:)=B(1,:)*b0;
%获取转换后的滤波器长度
S=size(B);
QB=zeros(S(1),S(2));
QA=QB;

%采用自编的E4_7_Qcoe函数完成级联滤波器的系数量化
for i=1:S(1)
[QB(i,:),QA(i,:)]=E4_7_Qcoe(B(i,:),A(i,:),qm);
end
%在命令窗口显示量化后的滤波器系数
QB,QA

image

1
2
3
4
QB =	94         	140			94
1024 162 1024
QA = 2048 -1213 268
1024 -953 585

上面的MATLAB代码对求出的IIR的滤波系数,然后再转化成级联型IIR的滤波系数,再进行量化,将含有小数的系数映射12bit。可以在FPGA中运算。

将系数带入差分方程得

image

x(n)为输入的待滤波数据,y1为一次滤波,y2位最终输出滤波后的数据。

image

Verilog设计

Verilog显示上述差分方程还是可以的,只需要进行延时和运算。其中乘法运算,直接用使用了有符号数乘号,也可用移位加来完成。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
module Cheyb2_Tap1st(
input clk,
input rst_n,
input signed [11:0] din,
output reg signed [11:0] dout
);

/*
2048y1(n) = 94[x(n)+x(n-2)] + 140x(n-1) - [-1213y1(n-1)+268y1(n-2)]
2048y2(n) = 2048[y1(n)+y1(n-2)] + 324y1(n-1) - [-1907y2(n-1)+117y2(n-2)]
*/

//-------------------------------------------------------
//first
reg signed [11:0] din_r1;
reg signed [11:0] din_r2;

wire signed [23:0] dout_r;
reg signed [12:0] dout_r1;
reg signed [12:0] dout_r2;

//wire signed [23:0] xmult0;
wire signed [23:0] xmult1;
wire signed [23:0] xmult2;
wire signed [23:0] xout;
wire signed [23:0] ymult1;
wire signed [23:0] ymult2;

always @(posedge clk or negedge rst_n)begin
if(rst_n == 1'b0)begin
din_r1 <= 12'b0;
din_r2 <= 12'b0;
end
else begin
din_r1 <= din;
din_r2 <= din_r1;
end
end

assign xmult1 = 94 * $signed(din + din_r2);
assign xmult2 = 140 * $signed(din_r1);

//assign xmult0 = {{6{din[11]}},din,6'd0}+{{7{din[11]}},din,5'd0}-{{11{din[11]}},din,1'd0}; //*94
//assign xmult1 = {{5{din_r1[11]}},din_r1,7'd0}+{{9{din_r1[11]}},din_r1,3'd0}+{{10{din_r1[11]}},din_r1,2'd0}; //*140
//assign xmult2 = {{6{din_r2[11]}},din_r2,6'd0}+{{7{din_r2[11]}},din_r2,5'd0}-{{11{din_r2[11]}},din_r2,1'd0}; //*94

//assign xout = xmult0 + xmult1 + xmult2;
assign xout = xmult1 + xmult2;

//-------------------------------------------------------
//
always @(posedge clk or negedge rst_n)begin
if(rst_n == 1'b0)begin
dout_r1 <= 12'b0;
dout_r2 <= 12'b0;
end
else begin
dout_r1 <= dout;
dout_r2 <= dout_r1;
end
end

assign ymult1 = 1213 * $signed(dout_r1);
assign ymult2 = 268 * $signed(dout_r2);

//assign ymult1 = {{2{dout_r1[11]}},dout_r1,10'd0}+{{5{dout_r1[11]}},dout_r1,7'd0}+{{6{dout_r1[11]}},dout_r1,6'd0}-{{11{dout_r1[11]}},dout_r1,1'd0}-{{12{dout_r1[11]}},dout_r1}; //*1213=1024+128+64-2-1
//assign ymult2 = {{4{dout_r2[11]}},dout_r2,8'd0}+{{9{dout_r2[11]}},dout_r2,3'd0}+{{10{dout_r2[11]}},dout_r2,2'd0}; //*268=256+8+4

assign dout_r = xout + ymult1 - ymult2;


always @(posedge clk or negedge rst_n)begin
if(rst_n == 1'b0)begin
dout <= 0;
end
else begin
dout <= $signed(dout_r) >>> 12;//除以2048
end
end

endmodule

我这里只贴出部分代码,查看更多代码,继续往下看。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
%test data
f1 = 1*10^6;%1MHz
f2 = 2.1 * 10^6;%2MHz
Fs = 8*10^6;%采样频率为8MHz
N = 12;
Len = 2000;
t = 0:1/Fs:(Len - 1)/Fs;
s = sin(2*pi*f2*t) + sin(2*pi*f1*t);
y = round(s / max(abs(s))*(2^11-1));%12bit量化
plot(t, y)
%grid;
b = signed2unsigned(y,N); %转换为无符号数输入

fid = fopen('sinx.coe','w'); %写到sin.coe文件,用来初始化sin_rom
fprintf(fid,'MEMORY_INITIALIZATION_RADIX=10;\n');
fprintf(fid,'MEMORY_INITIALIZATION_VECTOR=\n');
for i = 1:1:2000
fprintf(fid,'%d',y(i));
if i==2000
fprintf(fid,';');
else
fprintf(fid,',');
end
if i%15==0
fprintf(fid,'\n');
end
end
fclose(fid);

产生一个1Mhz和2Mhz的正弦叠加信号,给滤波模块。正弦叠加信号的coe文件产生方式MATLAB代码如上,产生方法和DDS信号发生器类似。

image

最终仿真结果如图

使用MATLAB FDA Tool生成HDL代码

还可以使用MATLAB FDA Tool生HDL代码

image

在MATLAB中输入命令,fdatool。

image

如图配置好滤波器参数

image

点击左侧边栏第三个按钮,Set Quantization Parameters

image

image

设置好如上图,点击Apply

image

然后上方,Targets——Generate HDl。

image

这个界面可配置文件类型(Verilog)、文件名、文件路径、Filter原理,根据你喜欢的样子。

image

这页配置端口列表等。

image

这页配置Test Bench

image

这页配置脚本,完后点击Generate,即可生成Verilog HDL代码。

关于FDA Tool的生成HDL 代码更详细的配置,请查阅Mathwork官网,hdlfilter和dsp_ug,这两个文档。在本订阅号这里也可以快捷获得,继续往下看。

image

生成好后找到对应文件,在顶层实例化,跑一下仿真,仿真结果如图。

MATLAB Generate HDL代码这种方式可以提高开发效率,但是机器生成的代码比人工写的耗费资源肯定多,不过一些复杂的算法为了提高开发速度,也会机器去生成。

在不久的未来,人工智能生成代码会越来越多,智能革命不同于前几次工业革命,是要取代人的脑力劳动。那时的社会会发生天翻地覆的变化,所以你是被操控AI还是被AI取代呢?最后瞎扯一点。

这个题目的相关文档和参考资料我已经上传至Github,以及 硅农小灶 知识星球。需要的朋友在微信订阅号后台,回复“IIR Filter”即可获得。

Reference

MathWorks官网——hdlfilter、 dsp_ug文档。

https://en.wikipedia.org/wiki/Infinite_impulse_response

《数字解调技术的MATLAB与FPGA实现——Xilinx/VHDL版》——杜勇

image

image