0%

FPGA-Verilog 设计FIR滤波器

FPGA/Verilog 设计FIR滤波器

[TOC]

前言

这应该是第一次的FPGA(DSP方向)的实战(也算不上)分享.也算是小班教学的其中一节课吧.
话不多说,先给大家介绍一下这次要干啥先:学过信号与系统的可以直接跳过基础知识…
系统环境:

matlab2018a
quartus16.0
Multisim
VScode
ubuntu18.04
EP4CE15F23C8

设计目标

设计一个10阶的FIR低通滤波器,滤波器的通带截止频率是2MHz,阻带截止频率是4MHz.

基础知识

数字滤波器

滤波器是一种对信号有处理作用的器件或电路,其主要作用是让有用信号尽可能无衰减地通过,对无用信号尽可能大地衰减.而数字滤波器就是一个按预定的有限精度算法实现的,将输入的数字信号转换为所要求的输出数字信号的线性时不变系统.
如果对上面加粗的莫得认识的话,那你就直接想成是用数字电路实现的滤波器得了.

FIR滤波器

FIR(Finite Impulse Response,FIR)滤波器,中文名叫有限冲激响应滤波器,在讲冲激响应之前,我们不妨预习/回顾一下信号与系统:
1

就是一个输入信号,经过一个系统,每个系统会有他对应的系统函数来处理输入信号,最后得到输出信号.对常见的滤波器来说,他们的理想系统函数往往是这些:

  1. 低通滤波器(红线部分),也是这篇blog要实现的东西
    2
  2. others
    3

注意看横坐标,我们可以看到,输入信号是时域进来的,但是系统函数是频域来的.而他们的命名也是跟频率有关系的.
这个时候我们就需要用到傅里叶变换了,至于详情,请看:
这是一个通俗易懂的傅里叶变换简介/教程

我们可以看到频域上的乘积就是时域上的卷积:
$$ y(n) = x(n)h(n) = \sum^N_{k=0}h(k)x(n-k) $$
用Z变换来就是:
$$H(z) = \sum^{N-1}_{n=0}h(n)z^{-n}$$
突然想起并不是在教DSP,所以我们了解到这里就够了,也就是说,我们只要把h(k)求出来,然后知道他可以有滤波的效果就完事了.

这里用的是FIR的横向结构:
3

我们可以看到主要有三个部分,第一部分是对输入信号的延时,第二部分是输入信号和抽头系数的相乘,第三部分是相乘结果的累加.这也是fpga的基本结构假设

具体实现

matlab-获取FIR抽头系数

在matlab命令行窗口打入filterDesigner(善用tab键):
4
大家大概可以看见面板下面有很多东西,分别是:

  1. 指定滤波器类型,是IIR还是FIR,用什么算法生成
  2. 滤波器阶数,是特定阶数还是自动生成最小阶数的
  3. 采样频率,信号频率,截止频率
  4. 衰减倍数
    当中的折衷关系不是这里的重点,根据设计要求,最后要改成这样:
    5
    记得按下方的design filter.

右上的图就是系统函数的形状了.

接下来点[file]->[Export]导出,快捷键ctrl+E,直接导出到workspace就可以了

然后在命令行窗口打:round(Num*512),意思是放大512倍并取整,方便我们后面在fpga中做定点数乘法

6

保存下来即可

matlab-产生波形数据

因为懒得接信号源和写ad模块驱动,这次实验直接在fpga中用一个rom的ip核储存波形数据
matlab代码如下:

  1. 生成波形

    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
    %产生两个正弦信号sin(x)和sin(8x)叠加后的信号,取128个点,将信号放大,
    %转换成无符号数据,存入ROM中作为信号源
    clear all
    clc
    depth = 128;
    width = 16;
    x = 0 : 2*pi/(depth-1) :2*pi;
    y = sin(x)+sin(8*x);
    plot(x,sin(x),'r') %红色为sin(x)函数
    hold on
    plot(x,sin(8*x),'g') %绿色为sin(8x)函数
    hold on
    plot(x,y,'b') %蓝色为生成的混合信号
    grid
    y = (y/2) * 32768;%将信号放大32768倍
    b = signed2unsigned(y,width); %转换为无符号数输入

    %下面函数重新新建一个脚本文件
    %需要调用了如下函数,将有符号数转换成无符号数
    function b = signed2unsigned(a,wl)
    %This function covert an signed integer number into an unsinged integer
    %number. a is the input vector while wl means the width of input number;
    %Example: a = [-2,-1,0,1];
    %signed2unsigned(a,3); THEN return [2,3,0,1]
    k = 2^(wl)*(a<0);
    b = k + a;;
    b = fix(b+0.5);
    for i = 1:length(a)
    if (b(i) == 65536)
    b(i) = 0;
    end
    end
  2. 编写mif文件

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    %编写mif文件
    fid = fopen('sinx.mif','wt'); %将信号写入一个.mif文件中
    fprintf(fid,'WIDTH=%d;\n',width);%写入存储位宽8位
    fprintf(fid,'DEPTH=%d;\n',depth);%写入存储深度1024
    fprintf(fid,'ADDRESS_RADIX=UNS;\n');%写入地址类型为无符号整型
    fprintf(fid,'DATA_RADIX=UNS;');%写入数据类型为无符号整型
    fprintf(fid,'CONTENT BEGIN\n');%起始内容
    for num=0 : 127
    fprintf(fid,'%d:%16.0f;\n',num,b(num+1));
    end
    fclose(fid);

运行整套代码我们也可以看见:
7
红色是2MHz的正弦波,绿色是8MHz的正弦波,蓝色是混合信号.

FPGA-导入波形数据

这里我们需要例化一个ROM模块来在fpga上预存前面生成的波形文件,打开quartus,[Tools]->[IP Catalog]
查找ROM,会找到一个在[Basic Function]->[On Chip Memory]下面的ROM:1-PORT(其实只要看见这个名字就可以了)双击例化

  1. 设置好位宽和数据长度
    8

  2. 由于对ROM没有别的要求,所以看直接next到下图这个位置:
    9
    填进刚刚在matlab里面生成的数据文件.
    然后直接finish就可以了

  3. 例化这个新建的ROM,顺便加时钟

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    //例化ROM模块
    data_rom ROM_Init
    (
    .address (address_rom ), //rom的地址端口
    .clock (CLK_50M ), //rom的时钟端口
    .q (data_rom ) //rom的数据端口
    );

    //时序电路,用来给rom_addr寄存器赋值
    always @ (posedge CLK_50M or negedge RST_N)
    begin
    if(!RST_N) //判断复位
    address_rom <= 7'd0; //初始化time_cnt值
    else
    address_rom <= address_rom + 2 ; //用来给time_cnt赋值
    end
    //这里为什么是 +2 后面再讨论

这样子,每一个时钟周期下面我们在data_rom下面都可以得到一个波形数据.

  1. 顺道把抽头系数定义一下
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    localparam COEFF1     	=	63;  
    localparam COEFF2 = 39;
    localparam COEFF3 = 48;
    localparam COEFF4 = 54;
    localparam COEFF5 = 59;
    localparam COEFF6 = 60;
    localparam COEFF7 = 59;
    localparam COEFF8 = 54;
    localparam COEFF9 = 48;
    localparam COEFF10 = 39;
    localparam COEFF11 = 63;

FPGA - FIR结构设计

由前面我们讲过,FIR有三个部分,有三大部分,第一部分是对输入信号的延时,第二部分是输入信号和抽头系数的相乘,第三部分是相乘结果的累加.所以我们也可以分三部分来写这个东西

part1-移位寄存器

我们可以用移位寄存器来达到输入信号的延时.

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
//pipeline 1 
always @ (posedge CLK_50M or negedge RST_N)
begin
if(!RST_N)begin
data_shift[0] <= 0;
data_shift[1] <= 0;
data_shift[2] <= 0;
data_shift[3] <= 0;
data_shift[4] <= 0;
data_shift[5] <= 0;
data_shift[6] <= 0;
data_shift[7] <= 0;
data_shift[8] <= 0;
data_shift[9] <= 0;
data_shift[10] <= 0;
end
else begin
data_shift[10] <= data_shift[9];
data_shift[9] <= data_shift[8];
data_shift[8] <= data_shift[7];
data_shift[7] <= data_shift[6];
data_shift[6] <= data_shift[5];
data_shift[5] <= data_shift[4];
data_shift[4] <= data_shift[3];
data_shift[3] <= data_shift[2];
data_shift[2] <= data_shift[1];
data_shift[1] <= data_shift[0];
data_shift[0] <= data_rom;
end
end

part2-乘法器

这里本来想例化乘法器ip来写的,但是鉴于很懒,所以就直接在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

always @ (posedge CLK_50M or negedge RST_N)
begin
if(!RST_N)
mul_data[0] <= 0;
else
mul_data[0] <= data_shift[0] * COEFF1;
end

always @(posedge CLK_50M or negedge RST_N)begin
if(!RST_N)
mul_data[1] <= 0;
else
mul_data[1] <= data_shift[1] * COEFF2;
end

………………………………………………………………………

always @(posedge CLK_50M or negedge RST_N)begin
if(!RST_N)
mul_data[9] <= 0;
else
mul_data[9] <= data_shift[9] * COEFF10;
end
always @(posedge CLK_50M or negedge RST_N)begin
if(!RST_N)
mul_data[10] <= 0;
else
mul_data[10] <= data_shift[10] * COEFF11;
end

part3-加法器

这里给出两种写法,顺道说明一下fpga里面的些许技巧

  1. 类c写法

    1
    2
    3
    4
    5
    6
    7
    always @(posedge CLK_50M or negedge RST_N)begin
    if(!RST_N)begin
    dout_r <= 0;
    end
    else
    dout_r <= mul_data[0]+mul_data[1]+mul_data[2]+mul_data[3]+mul_data[4]+mul_data[5]+mul_data[6]+mul_data[7]+mul_data[8]+mul_data[9]+mul_data[10];
    end
  2. 例化ip写法
    用的是parallel_add的ip,其实只要自己手动优化一下,也可以不用ip.因为他只能做8的倍数的,但是我的滤波器只有10阶.虽然有点浪费资源,但是简单演示的话,把其他位置0就完事了.例化过程较为简单,就不介绍了

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    wire [19:0] add_result;
    add unit_add(
    .clock (CLK_50M),
    .data0x (mul_data[0][21:6]),
    .data10x (mul_data[1][21:6]),
    .data11x (mul_data[2][21:6]),
    .data12x (mul_data[3][21:6]),
    .data13x (mul_data[4][21:6]),
    .data14x (mul_data[5][21:6]),
    .data15x (mul_data[6][21:6]),
    .data1x (mul_data[7][21:6]),
    .data2x (mul_data[8][21:6]),
    .data3x (mul_data[9][21:6]),
    .data4x (mul_data[10][21:6]),
    .data5x (16'b0),
    .data6x (16'b0),
    .data7x (16'b0),
    .data8x (16'b0),
    .data9x (16'b0),
    .result (add_result)
    );

至此,整个系统的各个关键部件就搭建完成了.

仿真结果

仿真波形

10
第一个波形是输入信号,
第二个波形是加法器的输出,
第三个波形是用上面第一种加法器写出来的波形,
第四个波形是用上面第二种加法器写出来的波形.

可以看出基本有那么一点点的滤波的作用,然而噪声还是非常的大,其中最主要的原因就在于,这个FIR滤波器只有十阶,再加上FIR本身的滤波特性,莫得办法啦

而事实上由于阶数太少了,以至于滤波效果还没有收到字长效应的影响,这也是让我始料未及的…

方案对比

写这个最主要的原因是为了回答一下那些说用类c写法来直接写fpga的大兄弟们…因为真的很多人学了之后都这样子…

方案一:类c写法

我们来看一下他的资源占用和最高运行速率:
11
12

总共用了857个逻辑单元,最高运行速率是59.2MHz

方案二:例化ip写法

13
14

总共用了797个逻辑单元,最高运行速率是190.59MHz.
还千万不要忘记了,这个ip是16个输入位的,而我们只用了其中的11个!

所以结果就是,我用的资源比你少,速度可以比你快,出来的波形还一样.虽然类C写法确实带来了不少的方便,但是他确实敌不过速度的制约,速度的制约出问题之后,显然时序约束就是天荒夜谈了.

结语

这一个应用或许是fpga在dsp上面较为简单的一个应用了,用fpga写dsp,还是要有扎扎实实的dsp基础才行.

这个也志在跑通一个流程吧,FIR其实有很多东西,考虑到篇幅,我都没有写出来.不过也莫得关系了.

期末愉快

如果你觉得有丶收获的话