FPGA设计实现FIR滤波器仿真验证_fpga实现fft如何仿真验证正确性

要求：设计一个15阶(长度为16)的低通线性相位FIR滤波器，采用blackman窗，截止频率500Hz,采样频率2KHz,采用FPGA全串行结构滤波器实现，滤波器系数量化位宽12bit，输入数据位宽12bit，输出数据位宽29bit,系统时钟16KHz。(数据速率刚好2KHz)

目录

1.首先采用Matlab求出符合要求的滤波器的系数，并观察12bit量化后的滤波器幅频响应曲线（代码如下）。

2.采用Matlab仿真产生待测试数据

3.编写Verilog代码实现FPGA设计滤波器结构

4.编写仿真测试文件

5.将仿真滤波器滤波后输出的数据写入txt文件中。

6.采用Matlab分析经过滤波后的数据，分析滤波性能。

7.总结

1.首先采用Matlab求出符合要求的滤波器的系数，并观察12bit量化后的滤波器幅频响应曲线（代码如下）。

%fir8serial.m
function hn = E4_7_Fir8Serial
N = 16;%滤波器长度
fs = 2000;%滤波器采样频率
fc = 500;%低通滤波器的截止频率
B = 12;%量化位数
 
%生成窗函数
w_black = blackman(N)';
%采用fir1函数设计FIR滤波器
b_black = fir1(N-1,fc*2/fs,w_black);
%量化滤波器系数
Q_black = round(b_black/max(abs(b_black))*(2^(B-1)-1));
hn = Q_black;
%转换成16进制数补码
Q_h = dec2hex(Q_black+2^B*(Q_black<0));
 
%求滤波器的幅频响应
m_black = 20*log10(abs(fft(b_black,1024)));m_black = m_black - max(m_black);%将最大值归位0
Ql_black = 20*log10(abs(fft(Q_black,1024)));Ql_black = Ql_black - max(Ql_black);
 
%设置幅频响应横坐标单位为Hz
x_f = [0:(fs/length(m_black)):fs/2];
m1 = m_black(1:length(x_f));
m2 = Ql_black(1:length(x_f));
 
%绘制幅频响应曲线
plot(x_f,m1,'-',x_f,m2,'-.');
xlabel('频率(Hz)');ylabel('幅度(dB)');
legend('未量化','12bit量化');
grid;

运行结果为

可以发现12bit量化符合要求

滤波器系数如下

2.采用Matlab仿真产生待测试数据

并通过设计的滤波器验证理论上的滤波效果同时将生成的测试数据写入txt文件中供后面Modelsim仿真时读取

代码如下：

%E4_7_NoiseAndCarrier.m 
f1 = 200 ;%信号1频率200Hz
f2 = 800 ;%信号2频率800Hz
Fs = 2000;%采样频率2KHz
N  = 12  ;%量化数据位宽12bit 输入数据12bit
 
%产生信号
t = 0:1/Fs:1;
c1 = 2*pi*f1*t;
c2 = 2*pi*f2*t;
s1 = sin(c1);%产生正弦信号
s2 = sin(c2);%产生正弦信号
s = s1 + s2;%产生两个单载波合成后的信号
%产生随机序列信号
noise=randn(1,length(t));%产生高斯白噪声序列
 
%归一化处理
s = s / max(abs(s));
noise = noise / max(abs(noise));
 
%12bit量化
Q_noise = round(noise*(2^(N-1)-1));
Q_s = round(s*(2^(N-1)-1));
 
%调用自己设计的滤波器函数对信号进行滤波
hn = fir8serial;
Filter_noise = filter(hn,1,Q_noise);
Filter_s     = filter(hn,1,Q_s);
 
%求信号的幅频响应
%滤波之前
m_noise = 20*log10(abs(fft(Q_noise,1024))); m_noise = m_noise - max(m_noise);
m_s     = 20*log10(abs(fft(Q_s,1024))); m_s = m_s - max(m_s);
%滤波之后
Fm_noise = 20*log10(abs(fft(Filter_noise,1024)));Fm_noise = Fm_noise - max(Fm_noise);
Fm_s     = 20*log10(abs(fft(Filter_s,1024)));Fm_s = Fm_s - max(Fm_s);
%滤波器本身的幅频响应
m_hn = 20*log10(abs(fft(hn,1024)));m_hn = m_hn -  max(m_hn);
 
%设置幅频响应的横坐标为Hz
x_f = [0:(Fs/length(m_s)):Fs/2];
%只显示正频率部分
mf_noise = m_noise(1:length(x_f));
mf_s     = m_s(1:length(x_f));
Fmf_noise = Fm_noise(1:length(x_f));
Fmf_s     = Fm_s(1:length(x_f));
Fm_hn = m_hn(1:length(x_f));
 
%绘制幅频响应曲线
subplot(211)
plot(x_f,mf_noise,'-.',x_f,Fmf_noise,'-',x_f,Fm_hn,'--');
xlabel('频率(Hz)');ylabel('幅度(dB)');title('Matlab仿真白噪声滤波前后的频谱');
legend('输入信号频谱','输出信号频谱','滤波器响应');
grid;
 
subplot(212);
plot(x_f,mf_s,'-.',x_f,Fmf_s,'-',x_f,Fm_hn,'--');
xlabel('频率(Hz)');ylabel('幅度(dB)');title('Matlab仿真合成单频信号滤波前后频谱');
legend('输入信号','输出信号','滤波器响应');
grid;
 
%将生成的数据以10进制格式写入txt文档
fid = fopen('D:\test\fir\4_7\E4_7_Int_noise.txt','w');
fprintf(fid,'%8d\r\n',Q_noise);
fprintf(fid,';');
fclose(fid);
 
fid = fopen('D:\test\fir\4_7\E4_7_Int_s.txt','w');
fprintf(fid,'%8d\r\n',Q_s);
fprintf(fid,';');
fclose(fid);
 
%将生成的数据以二进制形式写入txt文档
fid = fopen('D:\test\fir\4_7\E4_7_Bin_noise.txt','w');
for i = 1:length(Q_noise)
	B_noise = dec2bin(Q_noise(i)+(Q_noise(i)<0)*2^N,N);%N为位数
	for j=1:N
		if B_noise(j)=='1'
			tb = 1;
		else
			tb = 0;
		end 
		fprintf(fid,'%d',tb);
	end 
	fprintf(fid,'\r\n');
end 
fprintf(fid,';');
fclose(fid);
 
fid = fopen('D:\test\fir\4_7\E4_7_Bin_s.txt','w');
for i=1:length(Q_s)
	B_s = dec2bin(Q_s(i)+(Q_s(i)<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);

运行结果：

我们可以发现当输入信号为白噪声时输出800Hz以后的频率衰减已经达到-40dB了

当输入信号为200Hz和800Hz的叠加信号时，滤波后的输出只剩200Hz的信号了，800Hz被滤除了。

3.编写Verilog代码实现FPGA设计滤波器结构

/*
*
*@Author: X-Z
*@Date:2023-09-17 00:45:22
*@Function:全串行FIR数字滤波器 加法器位宽13bit,量化12bit扩展一位 乘法器输入位宽12bit和13bit输出数据位宽25bit
*/
 
/*
    设计一个15阶(长度2为16)的低通线性相位的FIR滤波器，采用布莱克曼窗，截止频率500Hz,采样频率2KHZ,
    采用FPGA实现全串行结构的滤波器，系数量化位数12Bit,输入数据位宽12bit，输出数据位宽29bit
    系统时钟16KHz
*/
 
module  FirFullSerial(
    input           wire                clk    ,//FPGA系统时钟频率16KHz
    input           wire                rst_n  ,
    input           wire signed [11:0]  Xin    ,//数据输入频率2KHz
    output          wire signed [28:0]  Yout   //滤波器输出数据
);
 
//例化有符号数乘法器IP核
    reg  signed [11:0]   coe  ;//滤波器为12bit量化数据
    wire signed [12:0]   add_s;//输入为12bit量化滤波器系数，两个对称项系数相加需要13bit位宽
    wire signed [24:0]   Mout ;
mul u_mult(
	.aclr    (~rst_n),
	.clock   (clk   ),
	.dataa   (coe   ),
	.datab   (add_s ),
	.result  (Mout  )
    );
 
//例化有符号数加法器IP核,对输入数据进行1位数据位的扩展，输出结果为13bit数据
    reg signed [12:0]  add_a;
    reg signed [12:0]  add_b;
add u_adder(
	.dataa  (add_a),
	.datab  (add_b),
	.result (add_s)
    );
 
//3位计数器，计数周期为8，刚好为数据输入的速率2KHZ
    reg [2:0] cnt;
    always @(posedge clk or negedge rst_n)begin
        if(!rst_n)
            cnt <= 3'd0;
        else
            cnt <= cnt + 1'd1; 
    end
 
//将数据移入移位寄存器Xin_Reg中
    reg  [11:0] Xin_Reg [15:0];//16个12bit位宽的存储器，用来存储滤波器系数
    reg  [3:0] i,j;
    always @(posedge clk or negedge rst_n)begin
        if(!rst_n)begin//初始化寄存器为0
            for(i=0;i<15;i=i+1)
                Xin_Reg[i] <= 12'd0;
        end
        else begin
            if(cnt==3'd7)begin
                for(j=0;j<15;j=j+1)
                    Xin_Reg[j+1] <= Xin_Reg[j];
                Xin_Reg[0] <= Xin;     
            end
        end
    end
 
//将对称系数的输入数据相加，同时将对应滤波器的系数送人乘法器
//需要注意的是以下程序只是用了一个加法器和一个乘法器资源
//以8倍速率调用乘法器IP核，由于滤波器长度为16，系数具有对称性，所以
//可以在一个时钟周期内完成所有8个滤波器系数与对应数据的乘法运算
//为保证加法运算不溢出，输入输出数据均扩展为13bit
//乘法器使用了1级流水线，加法器没有使用流水线会立即输出相加结果
/*
    12bit量化滤波器系数
    0          -3          15          46        -117        -263         590        2047     
    2047       590        -263        -117          46          15          -3       0
*/
 
    always @(posedge clk or negedge rst_n)begin
        if(!rst_n)begin
            add_a <= 13'd0;
            add_b <= 13'd0;
            coe   <= 12'd0;//滤波器系数
        end
        else begin
            case(cnt)
                3'd0 : begin
                    add_a <= {Xin_Reg[0][11],Xin_Reg[0]};//将数据扩展一位符号位
                    add_b <= {Xin_Reg[15][11],Xin_Reg[15]};
                    coe   <= 12'd0;//c0
                end
                3'd1 : begin
                    add_a <= {Xin_Reg[1][11],Xin_Reg[1]};
                    add_b <= {Xin_Reg[14][11],Xin_Reg[14]};
                    coe   <= -12'd3;//c1
                end
                3'd2 : begin
                    add_a <= {Xin_Reg[2][11],Xin_Reg[2]};
                    add_b <= {Xin_Reg[13][11],Xin_Reg[13]};
                    coe   <= 12'd15;//c2
                end
                3'd3 : begin
                    add_a <= {Xin_Reg[3][11],Xin_Reg[3]};
                    add_b <= {Xin_Reg[12][11],Xin_Reg[12]};
                    coe   <= 12'd46;//c3
                end
                3'd4 : begin
                    add_a <= {Xin_Reg[4][11],Xin_Reg[4]};
                    add_b <= {Xin_Reg[11][11],Xin_Reg[11]};
                    coe   <= -12'd177;//c4
                end
                3'd5 : begin
                    add_a <= {Xin_Reg[5][11],Xin_Reg[5]};
                    add_b <= {Xin_Reg[10][11],Xin_Reg[10]};
                    coe   <= -12'd263;//c5
                end
                3'd6 : begin
                    add_a <= {Xin_Reg[6][11],Xin_Reg[6]};
                    add_b <= {Xin_Reg[9][11],Xin_Reg[9]};
                    coe   <= 12'd590;//c6
                end 
                3'd7 : begin
                    add_a <= {Xin_Reg[7][11],Xin_Reg[7]};
                    add_b <= {Xin_Reg[8][11],Xin_Reg[8]};
                    coe   <= 12'd2047;//c7
                end
                default : begin
                    add_a <= 13'd0;
                    add_b <= 13'd0;
                    coe   <= 12'd0;//滤波器系数
                end
            endcase
        end
    end
 
//对滤波器系数与输入数据的乘法结果进行累加，并输出滤波后的数据
//考虑到乘法器及累加器的延时，需要计数器为2时对累加器清零，同时输出
//滤波器结果数据，类似的时延长度一方面可以通过精确计算获得，但更好的方法是
//通过行为仿真查看
    reg signed [28:0] sum ;
    reg signed [28:0] yout; 
    always @(posedge clk or negedge rst_n)begin
        if(!rst_n)begin
            sum  <= 29'd0;
            yout <= 29'd0;
        end
        else begin
            if(cnt==3'd2)begin
                yout <= sum ;//同时cnt==1时最后一个数运算完成但cnt==2时才有效；输出yout有效时cnt==3
                sum  <= Mout;//ct==2时输出第一个乘法器运算结果
            end
            else 
                sum <= sum + Mout; //经过1级流水线，cnt==1时全部加完cnt==2时输出累加结果有效
        end
    end
 
    assign Yout = yout;
 
endmodule

4.编写仿真测试文件

`timescale 1ns/1ps
module  FirFullSerial_tb();
 
    reg         clk     ;
    reg         rst_n   ;
    reg [11:0]  Xin     ;
    wire        clk_data;//数据时钟，速率为系统时钟clk速率的1/8
    wire [28:0] Yout    ;
 
//例化滤波器模块
FirFullSerial u_FirFullSerial(
    /*input           wire                */ .clk  (clk  )  ,//FPGA系统时钟频率16KHz
    /*input           wire                */ .rst_n(rst_n)  ,
    /*input           wire signed [11:0]  */ .Xin  (Xin  )  ,//数据输入频率2KHz
    /*output          wire signed [28:0]  */ .Yout (Yout )  //滤波器输出数据
);
 
    parameter  clk_period = 62500;//设置时钟信号周期(频率)：16KHz
    parameter  clk_period_data = clk_period * 8;
    parameter  clk_half_period = clk_period / 2;
    parameter  clk_half_period_data = clk_period_data / 2;
    parameter  data_num = 2000;//仿真数据长度
    parameter  time_sim = data_num * clk_period_data;//仿真时间
 
    initial begin
        clk      <= 1'b1;
        rst_n = 1'b0;
        #(5*clk_period)
        rst_n = 1'b1;
        //设置仿真时间
        #time_sim $finish;
        Xin = 12'd10;//输入信号赋初值
    end
 
    //产生时钟信号
    always #clk_half_period clk = ~clk;//16KHz
 
    reg [2:0] cn_clk = 0;
    always @(posedge clk) cn_clk <= cn_clk + 3'd1;
    assign clk_data = cn_clk[2];//时钟周期与系统时钟对齐；2KHz
 
    //从外部txt文件(SinIn.txt)中读入数据作为测试数据
    integer pattern;
    reg [11:0] stimuls [1:data_num];
    initial begin
        //文件必须放在"工程目录\simulation\modelsim"路径下
        $readmemb("E4_7_Bin_noise.txt",stimuls);
        //$readmemb("E4_7_Bin_s.txt",stimuls);
        pattern = 0;
        repeat(data_num)begin
            @(posedge clk_data);
            pattern = pattern + 1;
            Xin = stimuls[pattern];
        end
    end
 
    //将仿真数据dout写入外部txt文件中(out.txt)
    integer file_out;
    initial begin
        //放在"工程目录\simulation\modelsim"路径下
        file_out = $fopen("E4_7_Noiseout.txt");
        //file_out = $fopen("E4_7_sout.txt");
        if(!file_out)begin
            $display("could not open file!");
            $finish;
        end
    end
 
    wire rst_write;
    wire signed [28:0] dout_s;
    assign dout_s = Yout;//将dout转换成有符号数据
    assign rst_write = clk_data & rst_n;//产生写入时钟信号，复位状态下不写入数据
    always @(posedge rst_write)
        $fdisplay(file_out,"%d",dout_s);
 
endmodule 

5.将仿真滤波器滤波后输出的数据写入txt文件中。

Modelsim仿真结果：

通过仿真我们可以发现大致实现了滤波性能，但还是不够直观，因此我们采用Matlab进行分析滤波性能。

6.采用Matlab分析经过滤波后的数据，分析滤波性能。

代码如下：

%E4_7_NoiseAndCarrierOut.m 
f1 = 200 ;%信号1频率200Hz
f2 = 800 ;%信号2频率800Hz
Fs = 2000;%采样频率2KHz
N  = 12  ;%量化位数
 
%从文本文件中读取数据
%测试输入数据分别放在Noise_in和S_in变量中
fid = fopen('D:\test\fir\4_7\E4_7_Int_Noise.txt','r');
[Noise_in,N_n] = fscanf(fid,'%lg',inf);
fclose(fid);
 
fid = fopen('D:\test\fir\4_7\E4_7_Int_s.txt','r');
[S_in,N_n] = fscanf(fid,'%lg',inf);
fclose(fid);
 
%滤波后的输出结果数据放在Noise_out和S_out变量中
fid = fopen('D:\test\fir\4_7\FirFullSerial\prj\simulation\modelsim\E4_7_Noiseout.txt','r');
[Noise_out,N_count] = fscanf(fid,'%lg',inf);
fclose(fid);
 
fid = fopen('D:\test\fir\4_7\FirFullSerial\prj\simulation\modelsim\E4_7_sout.txt','r');
[S_out,N_count] = fscanf(fid,'%lg',inf);%以浮点数形式全部读取文件中的数据
fclose(fid);
 
%归一化处理
Noise_out = Noise_out / max(abs(Noise_out));
S_out     = S_out / max(abs(S_out));
Noise_in  = Noise_in / max(abs(Noise_in));
S_in      = S_in /max(abs(S_in));
 
%求信号的幅频响应
out_noise = 20*log10(abs(fft(Noise_out,1024)));out_noise = out_noise - max(out_noise);
out_s     = 20*log10(abs(fft(S_out,1024)));out_s = out_s - max(out_s);
in_noise  = 20*log10(abs(fft(Noise_in,1024)));in_noise = in_noise - max(in_noise);
in_s      = 20*log10(abs(fft(S_in,1024)));in_s = in_s - max(in_s);
 
%滤波器本身的幅频响应
hn = fir8serial;
m_hn = 20*log10(abs(fft(hn,1024)));m_hn = m_hn - max(m_hn);
 
%设置幅频响应横坐标单位为Hz
x_f = [1:(Fs/length(out_noise)):Fs/2];
%只显示正频率部分的幅频响应
mf_out_noise = out_noise(1:length(x_f));
mf_out_s     = out_s(1:length(x_f));
mf_in_noise  = in_noise(1:length(x_f));
mf_in_s      = in_s(1:length(x_f));
mf_hn        = m_hn(1:length(x_f));
 
%绘制幅频响应特性曲线
figure(1);
subplot(211);
plot(x_f,mf_in_noise,'--',x_f,mf_out_noise,'-',x_f,mf_hn,'-.');
xlabel('频率(Hz)');ylabel('幅度(dB)');
title('FPGA仿真白噪声滤波前后幅频响应曲线');
legend('输入信号频谱','输出信号频谱','滤波器响应');
grid;
subplot(212);
plot(x_f,mf_in_s,'--',x_f,mf_out_s,'-',x_f,mf_hn,'-.');
xlabel('频率(Hz)');ylabel('幅度(dB)');title('FPGA仿真单频叠加信号滤波前后幅频响应曲线');
axis([0 1000 -100 0]);
legend('输入信号幅频响应','输出信号频谱','滤波器响应');
grid;
 
%绘制时域波形
%设置显示数据的范围
t = 0:1/Fs:50/Fs;t=t*1000;
%length(t)
t_in_noise  = Noise_in(1:length(t));
t_in_s      = S_in(1:length(t));
t_out_s     = S_out(1:length(t));
t_out_noise = Noise_out(1:length(t));
 
 
figure(2);
subplot(211);
%length(hn);
plot(t,t_in_noise,'--',t,t_out_noise,'-');
xlabel('时间(ms)');ylabel('幅值');title('FPGA仿真白噪声滤波前后的时域波形');
legend('输入信号波形','输出信号波形');
grid;
subplot(212);
plot(t,t_in_s,'--',t,t_out_s,'-');
xlabel('时间(ms)');ylabel('幅值');
title('FPGA仿真单频叠加信号滤波前后时域波形');
legend('输入信号波形','输出信号波形');
grid;
 
 
 
 
 
 
 
 
 
 
 

运行结果为：

---

我们发现采用FPGA滤波前后的波形和Matlab仿真波形图基本一致，唯一不同的是Matlab的滤波效果更好一些，FPGA衰减没有Matlab仿真的大，这是因为FPGA中存在字长效应以及数据处理时存在截尾效应。因此可以发现理论与实际的差距，但最终设计满足要求。

注：自学完杜勇的FIR滤波器所写，能力有限，不喜勿喷。