今日你学左米啊

0%

IIR滤波器的FPGA实现

发表于 更新于 分类于 Valine：

IIR滤波器的FPGA实现

生产实习著
此篇实现的是直接1型7阶切比雪夫1型低通滤波器

[TOC]

IIR原理

IIR定义

IIR（Infinite Impulse Response）无限脉冲响应数字滤波器,也叫作递归滤波器.
由差分方程:
$$ y(n) = \sum^M_{i=0}b_i x(n-i) - \sum^N_{i=1}a_i y(n-i)$$
进行Z变换:
$$
Y(z)=\sum_{i=0}^{M} b_{i} z^{-i} X(z)-\sum_{i=1}^{N} a_{i} z^{-i} Y(z)
$$
得到系统函数:
$$
H(z)=\frac{Y(z)}{X(z)}=\frac{\sum_{i=0}^{M} b_{i} z^{-i}}{1+\sum_{i=1}^{N} a_{i} z^{-i}}
$$
至于为什么叫IIR,和FIR有什么区别,大家做一下z逆变换就知道了

模拟滤波器设计

要提及到模拟滤波器设计的理由很简单,就是因为设计IIR的一般是先设计好相应的模拟滤波器,然后再利用双线性变换转换成数字滤波器,在这里我们介绍一个切比雪夫1型滤波器:
他的幅度平方函数长这个样:
$$
G_{n}(\omega)=\left|H_{n}(j \omega)\right|=\frac{1}{\sqrt{1+\epsilon^{2} T_{n}^{2}\left(\frac{\omega}{\omega_{0}}\right)}}
$$
当$\omega = \omega_0$时有:
$$
\left|H\left(\omega_{0}\right)\right|=\frac{1}{\sqrt{1+\epsilon^{2}}}
$$
其中,$T_{n}\left(\frac{\omega}{\omega_{0}}\right)$是n阶切比雪夫多项式,定义为:
$$
\begin{cases}
T_{n}\left(\frac{\omega}{\omega_{0}}\right)=\cos \left(n \cdot \arccos \frac{\omega}{\omega_{0}}\right) ; 0 \leq \omega \leq \omega_{0}\
T_{n}\left(\frac{\omega}{\omega_{0}}\right)=\cosh \left(n \cdot \operatorname{arccosh} \frac{\omega}{\omega_{0}}\right) ; \omega>\omega_{0}
\end{cases}
$$
当0<x<1时，$T_{n}(x)$在0和1之间变化 ;当x>1时，cos−是虚数，所以$T_{n}(x)$像双曲余弦一样单调地增加;∣Hn(w)∣对于0≤w/w0≤1呈现出在1和$\frac1{2\epsilon^2}$之间的波动；而对于w/w0〉1单调地减小。需要用三个参量来确定该滤波器：ε，w0和N。在典型的设计中，用容许的通带波纹来确定ε，而用希望的通带截止频率来确定w0。然后选择合适的阶次n，以便阻带的技术要求得到满足。

上述的 滤波器设计原理看不懂没关系,因为模拟滤波器设计理论已经非常成熟了,所以我们可以借助matlab来解决
详见:Matlab的IIR Filter Design

除此之外,常见的模拟滤波器还有巴特沃斯,切比雪夫2型,椭圆型滤波器等,此处不做详述

双线性变换

这是利用反正切的方法将s域投影到z域,理解思路有很多种,大家不妨多方查阅一下资料,这里直接给出转换公式:
$$H(z)=H_{a}\left.(s)\right|{s=\frac{2}{T{s}} \frac{1-z^{-1}}{1+z^{-1}}}$$

IIR系统结构

实际上从IIR的定义我们就可以直接得到直接1型:
a

我们不妨将ab调转过来,合并延时器,就可以得到直接2型:
b

不妨将系统函数定义式做一点变换,拆开每项系数:
$$
H(z)=\frac{Y(z)}{X(z)}=\frac{\sum_{i=0}^{M} b_{i} z^{-i}}{1+\sum_{i=1}^{N} a_{i} z^{-i}}
$$
to:
$$
H(z)=\prod_{i=1}^{L} \frac{\gamma_{0 i}+\gamma_{1 i} z^{-1}+\gamma_{2 i} z^{-2}}{1+\beta_{1 i} z^{-1}+\beta_{2 i} z^{-2}}=A \prod_{i=1}^{L} \frac{1+\alpha_{1 i} z^{-1}+\alpha_{2 i} z^{-2}}{1+\beta_{1 i} z^{-1}+\beta_{2 i} z^{-2}}=A \prod_{i=1}^{L} H_{i}(z)
$$
得到每一个子系统$H_i(z)$,就可以用直接12型来表征子系统,然后再串起来,比如:
c

FPGA实现

这里实现的是直接1型7阶切比雪夫1型低通滤波器.思路很简单,直接将ab两路当成两个FIR滤波器就完事了

MATLAB 获取系数

获取系数有两种方法:
1. 利用函数cheby1,cheby2和butter等函数进行设计,建议回看我上面的matlab设计iir的帮助页面
2. 利用滤波器设计工具箱 filterDesigner
此处为了直观使用filterDesigner,设计完之后是这样的:
d
直接点Design Filter是设计出来的是直接2型的SOS(second-order sections)版也就是级联型.我们需要

  1. Edit-> Convert Structure 选Direct-Form I SOS
  2. Edit-> Convert to Single Section
  3. Ctrl +e (File->Export)导出系数

如果不慎生成了SOS的话,也可以在命令行中输入:

[Num,Den] = sos2tf(SOS,G);
进行转换,其中Num(numerator)为分子,Den(Denominator)为分母

然后对系数进行量化:

1
2
Num = round(Num*2^9);      %最大值是2.8多,建议思考为什么是512
Den = round(Den*2^9);

FPGA实现FIR滤波器

可以见我前面的FPGA/Verilog 设计FIR滤波器
我这里为了重定义系数和后续好改成SOS,做了一些小功夫和去除了累加器,自己写了个加法树.

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
module fir_7
#(parameter WIDTH=4'd12 , c0 = 12'b0 , c1 = 12'b0 , c2 = 12'b0 , c3 = 12'b0 , c4 = 12'b0 , c5=12'b0 , c6=12'b0 ,c7=12'b0 )
(
input rst_n ,
input clk   ,
input signed [WIDTH-1'b1:0] din,
output signed [2*WIDTH+2:0] dout
);
integer i=0;
integer j=0;

reg signed [WIDTH-1'b1:0] din_delay[7:0];
always@(posedge clk or negedge rst_n)
begin
	if(!rst_n)
		begin
		for(i=0;i<8;i =i+1)
			din_delay[i] <= 'b0;
		end
	else
		begin
		for(i=1;i<8;i =i+1)
			din_delay[i] <= din_delay[i-1];
		din_delay[0] <= din;
		end
end

reg signed [2*WIDTH-1'b1:0] din_mult[7:0];
always @ (posedge clk or negedge rst_n)
begin
	if(!rst_n)
		begin
		for(j=0;j<8;j =j+1)
			din_mult[j] <= 'b0;
		end
	else
		begin
		din_mult[0] <= din_delay[0]*c0;
		din_mult[1] <= din_delay[1]*c1;
		din_mult[2] <= din_delay[2]*c2;
		din_mult[3] <= din_delay[3]*c3;
		din_mult[4] <= din_delay[4]*c4;
		din_mult[5] <= din_delay[5]*c5;
		din_mult[6] <= din_delay[6]*c6;
		din_mult[7] <= din_delay[7]*c7;
		end
end

reg signed [2*WIDTH:0] din_mult1[3:0];
reg signed [2*WIDTH+1:0] din_mult2[1:0];
reg signed [2*WIDTH+2:0] din_mult3;

always @ (posedge clk or negedge rst_n)
begin
	if(!rst_n)
		begin
			din_mult1[0] <= 'b0;din_mult1[1] <= 'b0;din_mult1[2] <= 'b0;din_mult1[3] <= 'b0;
			din_mult2[0] <= 'b0;din_mult2[1] <= 'b0;din_mult3	 <= 'b0;
		end
	else
		begin
		din_mult1[0] <= din_mult[0] + din_mult[1];
		din_mult1[1] <= din_mult[2] + din_mult[5];
		din_mult1[2] <= din_mult[3] + din_mult[6];
		din_mult1[3] <= din_mult[4] + din_mult[7];
		
		din_mult2[0] <= din_mult1[0] + din_mult1[1];
		din_mult2[1] <= din_mult1[2] + din_mult1[3];
		
		din_mult3    <= din_mult2[0] + din_mult2[1];
		end
end
assign dout = din_mult3;
		
endmodule

FPGA根据系统框图连线

这里比较简单,就给出个关键部分吧:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
fir_7 
#(.WIDTH(4'd12),
 .c0(12'd1), .c1(12'd7), .c2(12'd21), .c3 (12'd36),
 .c4(12'd36), .c5(12'd21), .c6(12'd7) ,.c7(12'd1))
zero (
.rst_n(RST_N),
.clk(CLK_50M),
.din(Yin),
.dout(Yout)
);

fir_7 
#(.WIDTH(4'd12),
 .c0(12'd256) , .c1(-12'd510), .c2(12'd875), .c3 (-12'd971),
 .c4(12'd838), .c5(-12'd527), .c6(12'd231) , .c7(-12'd61))
poly (
.rst_n(RST_N),
.clk(CLK_50M),
.din(din),
.dout(Xout)
);

最终生成的框图是这样的:
e

MATLAB生成仿真数据

由于上几次的博客生成仿真数据的脚本不怎么好用,所以找了一个更好的jio本:

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
 %=============设置系统参数==============%
%f1=1e6;        %设置波形频率
f2=200e3;
f3=1000e3;
Fs=10e6;        %设置采样频率
L=4096;         %数据长度
N=12;           %数据位宽
%=============产生输入信号==============%
t=0:1/Fs:(1/Fs)*(L-1);
%y1=sin(2*pi*f1*t);
y2=0.5*sin(2*pi*f2*t);
y3=0.5*sin(2*pi*f3*t);
%y4=y1+y2+y3;
y4= y2+y3;
y_n=round(y4*(2^(N-2)-1));      %N比特量化;如果有n个信号相加，则设置（N-n）
%=================画图==================%
a=10;           %改变系数可以调整显示周期
stem(t,y_n);
axis([0 L/Fs/a -2^N 2^N]);      %显示
%=============写入外部文件==============%
fid=fopen('sinx.txt','w');    %把数据写入sin_data.txt文件中，如果没有就创建该文件  
for k=1:length(y_n)
    B_s=dec2bin(y_n(k)+((y_n(k))<0)*2^N,N);
    for j=1:N
        if B_s(j)=='1'
            tb=1;
        else
            tb=0;
        end
        fprintf(fid,'%d',tb);
    end
    fprintf(fid,'\r\n');
end
fprintf(fid,';');
fclose(fid);

写testbeach

这里只给出读取数据部分:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
integer i;   //数组坐标
reg signed [11:0] stimulus[1:data_num];  //数组形式存储读出的数据
initial 
begin
    RST_N = 1'b1;
	#60 RST_N = 1'b0;
	#60 RST_N = 1'b1; 
    $readmemb("sinx.txt", stimulus);  //将txt文件中的数据存储在数组中
    i = 0;
    repeat(data_num*10) begin   //重复读取数组中的数据
        i = i + 1;
        din = stimulus[i%data_num]; 
        #PERIOD;         //每个时钟读取一次
    end
	 $stop;
end

仿真结果

tg
后仿在1.2V 0度模型下最高运行速率180Mhz(板子最高200M)

结语

因为生产实习和电赛,最近不得不备一些这种小东西,会一直更新到7.26号,各位看官有兴趣的可以追踪一下.

如果csdn再不加入对verilog的支持,我就举家搬迁到简书或者自己搭一个博客了.

如果你觉得有丶收获的话