FPGA直接型IIR滤波器的实现

FPGA直接型IIR滤波器的实现

采用FPGA实现下列差分方程所述的IIR滤波器，并测试FPGA实现后的滤波效果，其中系统时钟频率为2KHZ、数据速率为2KHZ、输入数据位数为12比特。

900y(n)=13[x(n)+x(n-7)]+38[x(n-1)+x(n-6)]+74[x(n-2)+x(n-5)]+99[x[n-3]+x(n-4)]-[-1623y(n-1)+2047y(n-2)-1427y(n-3)+725y(n-4)-215y(n-5)+42y(n-6)-3y(n-7)]

IIR滤波器系数及运算字长的仿真_小小低头哥的博客-CSDN博客上回讲了如何使用MATLAB对IIR滤波器进行仿真，这回介绍如何用FPGA实现

1.分析直接型结构的实现方法

  对于零点系数的实现结构，其实可完全看成没有反馈结构的FIR滤波器，且可以利用系数对称的特点进一步减少乘法运算；单独查看系统极点的实现结构，即求取Y_out信号的过程也可以看成一个不带反馈结构的电路；整个IIR滤波器的闭环过程只在求取Y_sum的减法器，以及移位算法实现除法运算的过程中完成；在求取X_out、Y_out、Y_sum、信号的过程中，均可通过增加寄存器字长来实现全精度运算，唯一出现运算误差的环节是除法运算（以移位运算为例），以及除法运算后的截尾输出。当整个IIR滤波器稳定，且截位输出数据不出现溢出时，运算误差仅由除法运算产生。

图1 直接型IIR滤波器结构

2.零点系数的FPGA实现

  零点系数的FPGA实现可完全看成一个FIR滤波器，因此可完全采用实现FIR滤波器时所介绍的方法。需要注意的是，由于IIR滤波器的反馈结构特性，实现零点系数即极点系数的运算需要满足严格的时序要求。也就是说，要求在计算零、极点时不出现延时，这一结构特点实际上限制了系统的运行速度。为了提高系统的运行速度，零点系数的运算采用全并行结构。下面直接给出相应程序

//IIR零点系数的FPGA实现 
module ZeroParallel(
    input   rst,    //复位信号
    input   clk,    //时钟信号
    input   signed[11:0] Xin,   //输入数据x(n) 频率为2KHZ
    
    output  signed[20:0] Xout   //滤波后的输出数据
);

//将数据保存在寄存器中
reg signed[11:0] Xin_reg[6:0];  //分别保存x(n-1)-x(n-7)
reg [3:0] i,j;
always @(posedge clk or negedge rst)
    if(!rst) begin
        for(i=0;i<7;i=i+1) begin
            Xin_reg[i] <= 12'd0;
        end
    end
    else begin //每过一个时钟，说明计算完了一个输入数据对应的输出，数据位移动，载入新数据
        for(j=0;j<6;j=j+1) begin
            Xin_reg[j+1] <= Xin_reg[j]; 
        end
        Xin_reg[0] <= Xin;
    end
    
//将对称系数的输入数据相加  因为要在一个时钟周期内完成，所以这里只能用线网wire型
wire signed[12:0] Add_reg[3:0];
assign Add_reg[0] = {Xin[11],Xin} + {Xin_reg[6][11],Xin_reg[6]};//x(n)+x(n-7) 补码形式先扩展再相加
assign Add_reg[1] = {Xin_reg[0][11],Xin_reg[0]} + {Xin_reg[5][11],Xin_reg[5]};
assign Add_reg[2] = {Xin_reg[1][11],Xin_reg[1]} + {Xin_reg[4][11],Xin_reg[4]};
assign Add_reg[3] = {Xin_reg[2][11],Xin_reg[2]} + {Xin_reg[3][11],Xin_reg[3]};

//采用移位运算及加法运算实现乘法运算
wire signed[20:0] Mult_Reg[3:0];
assign Mult_Reg[0] = {{6{Add_reg[0][12]}},Add_reg[0],2'd0}    //将13位的有符号数扩展为21位 再左移两位 *4
                      + {{7{Add_reg[0][12]}},Add_reg[0],1'd0} // *2
                      + {{8{Add_reg[0][12]}},Add_reg[0]};     // *1 相加后总的就是*7
assign Mult_Reg[1] = {{4{Add_reg[1][12]}},Add_reg[1],4'd0} + {{6{Add_reg[1][12]}},Add_reg[1],2'd0} + 
                     {{8{Add_reg[1][12]}},Add_reg[1]};        //*21    
assign Mult_Reg[2] = {{3{Add_reg[2][12]}},Add_reg[2],5'd0} + {{5{Add_reg[2][12]}},Add_reg[2],3'd0} + 
                     {{7{Add_reg[2][12]}},Add_reg[2],1'd0};   //*42 
assign Mult_Reg[3] = {{3{Add_reg[3][12]}},Add_reg[3],5'd0} + {{4{Add_reg[3][12]}},Add_reg[3],4'd0} + 
                     {{5{Add_reg[3][12]}},Add_reg[3],3'd0};   //*56      

//对滤波器系数与输入数据的乘法结果进行累加，并输出滤波后的数据
assign Xout = Mult_Reg[0] + Mult_Reg[1] + Mult_Reg[2] + Mult_Reg[3];
endmodule

  Add_reg为什么是13位？因为两个12位的Xin_reg数相加，顶多是13位。

  Mult_Reg为什么是21位，因为最大值为56*Add_reg(13位)，而56无符号数需要6位表示，有符号数需要7位表示。因此至少需要7+13=20位。

  Xout为什么是21位，如IIR系统公式所示，Xout=x(n) * 2 * (7+21+42+56)=x(n) * 252 = x(n)(12位) * 252（有符号数9位）= 12+9=21位宽。

3.极点系数的FPGA实现

  零点系数的FPGA实现可完全看成一个FIR滤波器，因此可完全采用实现FIR滤波器时所介绍的方法。极点系数的计算涉及反馈结构，如何实现呢？可以将其分解为两部分，即
$$512y(n)=Zero(n)-Pole(n)$$
式中，Zero(n)是关于x的函数（已在第二节中实现），Pole(n)是关于y的函数。

  因此可以将极点系数的运算也堪称一个典型的乘加运算。==需要注意的是，为保证严格的时序特性，乘法器IP核不能使用输入\输出带有寄存器的结构。==相关实现代码如下

//极点系数的FPGA实现
module PoleParallel(
    input rst,
    input clk,
    input signed[11:0] Yin, //输入数据 y(n)
    
    output signed[25:0] Yout //滤波后的输出数据
);

//将数据存储在寄存器中
reg signed[11:0] Yin_reg[6:0];//分别存储y(n-1)-y(n-7)
reg [3:0] i;
always @(posedge clk or negedge rst)
    if(!rst) begin
        for(i=0;i<7;i=i+1) begin
            Yin_reg[i] <= 12'd0;
        end
    end
    else begin
        for(i=0;i<6;i=i+1) begin
            Yin_reg[i+1] <= Yin_reg[i];
        end
       Yin_reg[0] <= Yin; 
    end

//乘加运算
wire signed[11:0] coe[7:0];     //滤波器系数位12比特量化数据
wire signed[22:0] Mult_Reg[6:0];//乘法器输出为23比特数据
//assign coe[0] = 12'd512;
assign coe[1] = -12'd922;
assign coe[2] = 12'd1163;
assign coe[3] = -12'd811;
assign coe[4] = 12'd412;
assign coe[5] = -12'd122;
assign coe[6] = 12'd24;
assign coe[7] = -12'd2;
//建议使用没有流水线的乘法器IP核，此处简便起见，就直接用了*号
genvar j;
generate 
    for(j=0;j<7;j=j+1) begin:en
        assign Mult_Reg[j] = coe[j+1] * Yin_reg[j];
    end
endgenerate

//对滤波器系数与输入数据的乘法结果进行累加，并输出滤波后的数据
assign Yout = {{3{Mult_Reg[0][22]}},Mult_Reg[0]} + {{3{Mult_Reg[1][22]}},Mult_Reg[1]}
                + {{3{Mult_Reg[2][22]}},Mult_Reg[2]} + {{3{Mult_Reg[3][22]}},Mult_Reg[3]}
                + {{3{Mult_Reg[4][22]}},Mult_Reg[4]} + {{3{Mult_Reg[5][22]}},Mult_Reg[5]}
                + {{3{Mult_Reg[6][22]}},Mult_Reg[6]} ;

endmodule

4.顶层文件的设计

实现了IIR滤波器零、极点系数运算后，顶层文件的设计就变得简单了。如下所示

module Direct_Filter_IIR(
    input   rst,
    input   clk,
    input   signed[11:0] din,   //数据输入频率为2KHZ
    
    output  signed[11:0] dout   //滤波后的输出数据
);

//实例化IIR滤波器零点系数及极点系数运算模块
wire signed[20:0] Xout;
ZeroParallel u_ZeroParallel(
     .rst       (rst),    
     .clk       (clk),    
     .Xin       (din),
     .Xout      (Xout)
);

wire signed[11:0] Yin;
wire signed[25:0] Yout;
PoleParallel u_PoleParallel(
    .rst        (rst),
    .clk        (clk),
    .Yin        (Yin), 
    .Yout       (Yout)
);
wire signed[25:0] Ysum;
assign Ysum = {{5{Xout[20]}},Xout} - Yout;  //求IIR公式右边的结果 先将21位的补码扩展为26位的补码
//因为IIR滤波器系数中a(1)=512,需将加法结果除以512，采用右移9比特的方法实现除法运算
wire signed[25:0] Ydiv;                   //得到y(n) 将Ysum除以a(1)=512 
assign Ydiv = {{9{Ysum[25]}},Ysum[25:9]}; //将26位的补码右移9位，为保证右移后还是26位，符号位扩展成前9位
//根据仿真结果可知，IIR滤波器的输出范围与输入范围相同，输入脉冲函数0-1 输出<1 因此可直接进行截尾输出
assign Yin = (rst?Ydiv[11:0]:12'd0);
assign dout = Yin;

endmodule  

5.FPGA测试仿真

5.1 编写M文件，生成二进制输入测试数据

  测试数据生成文件用的与Fir那一章一致，就不重复了

5.2编写Test Bench文件，用Modelsim进行仿真

`timescale 1 ns/ 1 ns
module Direct_Filter_IIR_vlg_tst();
// constants                                           
// test vector input registers
reg clk;
reg [11:0] din;
reg rst;
// wires                                               
wire [11:0]  dout;

// assign statements (if any)                          
Direct_Filter_IIR i1 (
// port map - connection between master ports and signals/registers   
	.clk(clk),
	.din(din),
	.dout(dout),
	.rst(rst)
);

parameter clk_period = 500000;     //设置时钟信号周期
parameter clk_half_period = clk_period / 2; 
parameter data_num = 2000;      //仿真数据长度
parameter time_sim = data_num * clk_period / 2; //仿真时间

initial                                                
begin                                                  
    //设置时钟信号的初值
    clk = 1;
    //设置复位信号
    rst = 0;
    #10000 rst = 1;
    //设置仿真时间
    #time_sim $finish;
    //设置输入信号的初值
    din = 12'd10;    
end

//产生时钟信号                                                    
always  #clk_half_period clk = ~clk;
    
//从外部文本文件读入数据作为测试激励‘
integer Pattern;
reg signed[11:0] stimulus[1:data_num];
initial    
begin                                                                  
       //$readmemb("Bin_noise.txt",stimulus);
        $readmemb("E4_8_Bin_s.txt",stimulus);
        Pattern = 0;
        repeat(data_num) begin
            Pattern = Pattern + 1;
            din = stimulus[Pattern];
            #clk_period;
        end
end    

//将仿真数据dout写入外部文本文件中(Sout.txt)
integer file_out;
initial 
begin
    @(posedge rst);	//等待复位完成   
    //file_out = $fopen("Noiseout.txt");
    file_out = $fopen("Sout.txt");
    if(!file_out) begin
        $display("could not open file");
        $finish;
    end
end
wire rst_write;
wire signed[11:0] dout_s;
assign dout_s = dout;       //将dout转换成有符号数据
assign rst_write = clk & (rst);    //产生写入时钟信号，复位信号状态时不写入数据
always @(posedge rst_write)
    $fdisplay(file_out,"%d",dout_s);
    
endmodule   

仿真结果如图2所示

图2 Modelsim 仿真结果

初步判断结果还可以

5.3编写M文件，用MATLAB分析FPGA仿真后的数据

直接上结果了

图3 MATALB 读取FPGA输出数据的频谱仿真结果

图4 MATALB 读取FPGA输出数据的时域仿真结果

经过比较，FPGA仿真的数据与MATLAB直接仿真的信号相同，合成的信号也变成了单频的信号，从而可知，整个IIR滤波器的设计正确无误。