基于FPGA的FIR滤波器的实现（4）— 串行结构FIR滤波器的FPGA代码实现_verilog的滤波器量化

前言

前面已经详尽的介绍了有关FIR滤波器的matlab实现，使用matlab生成了FIR滤波器设计所需要的响应参数，本章开始讲解如何使用Verilog语言设计FIR滤波器

一、量化滤波器系数

对FPGA等硬件处理平台来讲，为防止运算时数据溢出，通常将输入数据及中间变量限制在绝对值小于1的范围内，FPGA处理的数据均是以二进制形式存放的数据，数据的小数点位置完全是由人为定义的，当定在紧靠最高位右边时，数据的表示范围即为绝对值小于1的小数。将一组数据先进行归一化处理，而后乘以一个整数因子，再进行四舍五入截尾处理，即可得到整数的量化数据。经上述处理后的数据，将小数点移至最高位的右边，即是满足要求的量化后的数据，最后利用matlab进制转换函数，便可以将整数类型的数据转换成二进制或十六进制数据，再由FPGA设计使用。

1、量化编程

eg：设计一个高通最优FIR滤波器，过渡带为1000-1500Hz，采样频率为8000Hz，通带波纹最大为0.01，阻带波纹最大为0.001，绘图比较12位、14位量化系数以及无量化时的幅度响应曲线，同时写入txt文件中。
需要使用firpm函数设计最优滤波器，而firpm函数无法确定滤波器的阶数，因此需要kaiserord函数求满足设计需求的最小滤波器阶数。
M程序：

function hn=filter_coe;
fs=8000;                    %采样频率
fc=[1000 1500];              %过滤带
mag=[0 1];                   %窗函数的理想滤波器幅度
dev=[0.001 0.01];            %波纹，前是阻带波纹，后是通带波纹
[n,wn,beta,ftype]=kaiserord(fc,mag,dev,fs); %获取凯塞窗参数
fpm=[0 fc(1)*2/fs fc(2)*2/fs 1];  %firpm函数的频段向量
magpm=[0 0 1 1 ];             %firpm函数的幅值向量

%设计最优滤波器
h_pm=firpm(n,fpm,magpm);

%滤波系数进行量化
h_pm12=round(h_pm/max(abs(h_pm))*(2^11-1)); %12位量化
h_pm14=round(h_pm/max(abs(h_pm))*(2^13-1)); %14位量化

%进制转换
q14_hex_pm=dec2hex(h_pm14+2^14*(h_pm14<0));

%将生成的滤波器系数数据写入FPGA所需的文件中
fid=fopen('D:\matlab work\fir1_1\filterCoe\filter_coe_1.txt','w');
fprintf(fid,'radix=10;\r\n');
fprintf(fid,'coefdata=\r\n');
fprintf(fid,'%8d\r\n',h_pm);fprintf(fid,';');fclose(fid);

%求滤波器的幅频响应
m_pm=20*log(abs(fft(h_pm,1024)))/log(10);m_pm=m_pm-max(m_pm);
m_pm12=20*log(abs(fft(h_pm12,1024)))/log(10);m_pm12=m_pm12-max(m_pm12);
m_pm14=20*log(abs(fft(h_pm14,1024)))/log(10);m_pm14=m_pm14-max(m_pm14);

%设置幅频响应的横坐标单位为Hz
x_f=[0:fs/length(m_pm):fs/2];

%只显示正频率部分的幅频响应
mf_pm=m_pm(1:length(x_f));
mf_pm12=m_pm12(1:length(x_f));
mf_pm14=m_pm14(1:length(x_f));

%绘制幅频响应曲线
plot(x_f,mf_pm,'-',x_f,mf_pm12,'-',x_f,mf_pm14,'--');
xlabel('频率(Hz)');ylabel('幅度(dB)');
legend('未量化','12位量化','14位量化');
grid;
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

2、频谱图

如图所示，量化位数对滤波器的阻带纹波有较大影响，且量化位数越高，影响越小。

二、串行结构的FPGA实现

FIR滤波器实际上是一个乘累加运算，且乘累加运算的次数由滤波器阶数来决定，滤波器系数均成一定的对称性，就可以将运算减少一半，减少硬件资源。
串行结构，串行实现滤波器的累加运算，将每级延时单元与相应系数的乘积结果进行累加后输出，因此整个滤波器实际上只需要一个乘法器运算单元。串行结构可分为全串行和半串行，全串行占用更少的加法器资源。

设计实例

eg：设计一个15阶低通滤波器，采用布莱克曼窗函数设计，截止频率为500Hz采样频率为2000Hz；采用FPGA实现全串行结构滤波器，系数量化位数为12bit，输入数据位宽为12bit，输出数据位宽为29bit

1. 使用matlab设计出滤波器系数

M程序：

function hn=black_fpga;
N=16;                   %滤波器阶数
fs=2000;                %采样频率
fc=500;                 %截止频率

w_black=blackman(N)';   %生成布莱克曼窗

b_black=fir1(N-1,fc*2/fs,w_black);  %使用fir1函数生成滤波器

b_black12=round(b_black/max(abs(b_black))*(2^11-1));    %量化系数为12位
hn=b_black12;

%转化成16进制数补码
Q_h=dec2hex(b_black12+2^12*(b_black12<0))

m_black=20*log(abs(fft(b_black,512)))/log(10);  %生成幅频响应
m_black=m_black-max(m_black);
m_black12=20*log(abs(fft(b_black12,512)))/log(10);    %量化系数为12的幅频响应
m_black12=m_black12-max(m_black12);

%写入文件
fid=fopen('D:\matlab work\fir1_1\filterCoe\filter_low_12.txt','w');
fprintf(fid,'%8d\r\n',Q_h);fprintf(fid,';');fclose(fid)

x_f=[0:fs/length(m_black):fs/2];    %设置幅频响应单位为Hz

%只显示正频率
m1=m_black(1:length(x_f));           
m2=m_black12(1:length(x_f));     

plot(x_f,m1,'-',x_f,m2,'-');
xlabel('频率(Hz)');ylabel('幅度(dB)'); 
legend('未量化','12位量化')
grid;                               %绘图
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

频率-幅度响应图

2. 使用matlab仿真滤波器测试数据

仿真测试数据经滤波器滤波后的输出数据。仿真生成12bit量化的采样频率为2000Hz的高斯白噪声、两个频率分别为200Hz和800Hz信号的合成信号，并将数据转换成二进制数据写入文本文件
M程序：

f1=200;     %信号1：200Hz
f2=800;     %信号2：800Hz
fs=2000;    %采样频率
N=12;       %量化位数

%产生信号
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));   %产生高斯白噪声序列

%归一化处理
noise=noise/max(abs(noise));    
s=s/max(abs(s));

%12bit量化
Q_noise=round(noise*(2^(N-1)-1));
Q_s=round(s*(2^(N-1)-1));

%调用自己设计的滤波器函数对信号进行滤波
hn=black_fpga;
filter_noise=filter(hn,1,Q_noise);
filter_s=filter(hn,1,Q_s);

%求信号的幅频响应
m_noise=20*log(abs(fft(Q_noise,1024)))/log(10);
m_noise=m_noise-max(m_noise);
m_s=20*log(abs(fft(Q_s,1024)))/log(10);
m_s=m_s-max(m_s);

%滤波后的幅频响应
Fm_noise=20*log(abs(fft(filter_noise,1024)))/log(10);
Fm_noise=Fm_noise-max(Fm_noise);
Fm_s=20*log(abs(fft(filter_s,1024)))/log(10);
Fm_s=Fm_s-max(Fm_s);

%滤波器本身的幅频响应
m_hn=20*log(abs(fft(hn,1024)))/log(10);
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('白噪声信号滤波前后的频谱');
legend('输入信号频率谱','输出信号频谱','滤波器响应');
grid;

subplot(212)
plot(x_f,mf_s,'-',x_f,Fmf_s,'-',x_f,Fm_hn,'-');
xlabel('频率(Hz)');ylabel('幅度(dB)');title('合成信号滤波前后的频谱');
legend('输入信号频率谱','输出信号频谱','滤波器响应');
grid;

%将生成的数据以十进制格式写入txt文件中
fid=fopen('D:\matlab work\fir1_1\filterCoe\noise.txt','w');
fprintf(fid,'%8d\r\n',Q_noise);
fprintf(fid,';');
fclose(fid);

fid=fopen('D:\matlab work\fir1_1\filterCoe\sin.txt','w');
fprintf(fid,'%8d\r\n',Q_s);
fprintf(fid,';');
fclose(fid);

%将生成的数据以二进制格式写入txt文件中
fid=fopen('D:\matlab work\fir1_1\filterCoe\noise_B.txt','w');
for i=1:length(Q_noise)
    B_noise=dec2bin(Q_noise(i)+(Q_noise(i)<0)*2^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:\matlab work\fir1_1\filterCoe\sin_B.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);
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
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111

频率-幅度响应图

3. 使用Verilog编写串行结构的FIR滤波器

Verilog代码：（需要先将matlab产生的noise_B和sin_B添加到工程目录下的simulation/modelsim文件夹中）

module FirFullSerial(
	input wire clk,
	input wire rst_n,				
	input signed [11:0]Xin,		//数据输入频率为2khz
	output signed [28:0]Yout	//滤波后的输出数据
);

	reg signed [11:0]coe;		//滤波器为12bit量化数据
	
	wire signed [12:0]add_s;	//输入为12bit量化数据，两个对称系数相加需要13bit存储
	wire signed [24:0]mout;

	mult mult_inst(
		.clock(clk),
		.dataa(coe),
		.datab(add_s),
		.result(mout)
	);
	
	//实例化有符号数相加IP核，对输入数据进行1位符号位扩展，输出结果为13bit数据
	reg signed[12:0]add_a;
	reg signed[12:0]add_b;
	
	adder adder_inst(
		.dataa(add_a),
		.datab(add_b),
		.result(add_s)
	);
	
	//3位计数器，计数周期为8，为输入数据速率
	reg [2:0]count;
	always @(posedge clk or negedge rst_n)
	if(!rst_n)
		count <= 3'b0;
	else
		count <= count+1'b1;
		
	//将数据存入移位寄存器Xin_reg中
	reg [11:0]Xin_reg[15:0];
	reg [3:0]i,j;
	always @(posedge clk or negedge rst_n)
	if(!rst_n)
		begin
			for(i=0;i<15;i=i+1)
				Xin_reg[i] = 1'b0;
		end
	else
		begin
			if(count == 7)
				begin
					for(j=0;j<15;j=j+1)
						Xin_reg[j+1] <= Xin_reg[j];
					Xin_reg[0] <= Xin;
				end
		end
		
	//将对称叙述的输入数据相加，同时将对应的滤波器系数送入乘法器
	//需要注意的是，下面程序只是用了一个加法器及一个乘法器资源
	//以8倍数据速率调用乘法器IP核，由于滤波器长度为16，系数具有对称性，
	//故可在一个数据周期内完成所有8个滤波器系数与数据的乘法运算
	//为了保证加法运算不溢出，输入、输出数据均扩展为13bit
	always @(posedge clk or negedge rst_n)
	if(!rst_n)
		begin
			add_a <= 13'd0;
			add_b <= 13'd0;
			coe <= 12'd0;
		end
	else
		begin 
			if(count==3'd0)
				begin
					add_a <= {Xin_reg[0][11],Xin_reg[0]};
					add_b <= {Xin_reg[15][11],Xin_reg[15]};
					coe <= 12'h000;		//c0
				end
			else if(count==3'd1)
				begin
					add_a <= {Xin_reg[1][11],Xin_reg[1]};
					add_b <= {Xin_reg[14][11],Xin_reg[14]};
					coe <= 12'hffd;		//c1
				end
			else if(count==3'd2)
				begin
					add_a <= {Xin_reg[2][11],Xin_reg[2]};
					add_b <= {Xin_reg[13][11],Xin_reg[13]};
					coe <= 12'h00f;		//c2
				end
			else if(count==3'd3)
				begin
					add_a <= {Xin_reg[3][11],Xin_reg[3]};
					add_b <= {Xin_reg[12][11],Xin_reg[12]};
					coe <= 12'h02e;		//c3
				end
			else if(count==3'd4)
				begin
					add_a <= {Xin_reg[4][11],Xin_reg[4]};
					add_b <= {Xin_reg[11][11],Xin_reg[11]};
					coe <= 12'hf8b;		//c4
				end
			else if(count==3'd5)
				begin
					add_a <= {Xin_reg[5][11],Xin_reg[5]};
					add_b <= {Xin_reg[10][11],Xin_reg[10]};
					coe <= 12'hef9;		//c5
				end
			else if(count==3'd6)
				begin
					add_a <= {Xin_reg[6][11],Xin_reg[6]};
					add_b <= {Xin_reg[9][11],Xin_reg[9]};
					coe <= 12'h24e;		//c6
				end
			else
				begin
					add_a <= {Xin_reg[7][11],Xin_reg[7]};
					add_b <= {Xin_reg[8][11],Xin_reg[8]};
					coe <= 12'h7ff;		//c7
				end
		end
		
	//对滤波器系数与输入数据的乘法结果进行累加，并输出滤波后的数据
	//考虑到乘法器及累加器的延时，需要计数器为2时对累加器清零，同时输出滤波器结果数据
	//类似的时延长度一方面可通过精度计算获取，但更好的方法是通过行为仿真查看
	reg signed[28:0]sum;
	reg signed[28:0]yout;
	always @(posedge clk or negedge rst_n)
	if(!rst_n)
		begin
			sum <= 29'd0;
			yout <= 29'd0;
		end
	else
		begin
			if(count==3'd2)
				begin
					yout <= sum;
					sum =29'd0;
					sum = sum+mout;
				end
			else
				sum = sum+mout;
		end
	
	assign Yout = yout;

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
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
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146

仿真测试模块：

`timescale 1ns/1ns
module FirFullSerial_tb;

	reg clk;
	reg rst_n;		
	reg [11:0]Xin;
	wire [28:0]Yout;
	wire clk_data;   //数据时钟，速率为时钟的八分之一

	FirFullSerial FirFullSerial_inst(
		.clk(clk),
		.rst_n(rst_n),				
		.Xin(Xin),		//数据输入频率为2khz
		.Yout(Yout)	//滤波后的输出数据
	);
	
	parameter clk_period = 62500;		//设置时钟信号周期/频率：16KHz
	parameter clk_period_data = clk_period*8;
	parameter data_num = 2000;			//仿真数据长度
	parameter time_sim = data_num*clk_period;	//仿真时间
	
	initial clk=1'b1;
	always #(clk_period/2) clk=~clk;
	
	initial begin
		rst_n=1'b0;
		#20000 rst_n = 1'b1;
		#time_sim $stop;
		Xin = 12'd0;
	end
	
	reg [2:0]cn_clk = 3'd0;
	always @(posedge clk) 
		cn_clk <= cn_clk + 1'b1;
	
	assign clk_data = cn_clk[2];
	
	//从外部文件读入数据作为测试激励
	integer Pattern;
	reg [11:0]stimulus[1:data_num];
	initial begin
		$readmemb("noise_B.txt",stimulus);
		//$readmemb("sin_B.txt",stimulus);
		Pattern = 0;
		repeat(data_num)
			begin
				@(posedge clk_data)
				Pattern = Pattern + 1;
				Xin = stimulus[Pattern];
			end
	end
		
	//将仿真数据dout写入外部文件中
	integer file_out;
	initial begin
		file_out = $fopen("noise_out.txt");
		//file_out = $fopen("sout.txt");
		if(!file_out)
		begin
			$display("could not open file!");
			$stop;
		end
	end
	
	wire rst_write;
	wire signed [28:0]dout_s;
	assign dout_s = Yout;
	assign rst_write = clk_data & (rst_n);
	always @(posedge rst_write)
		$fdisplay(file_out,"%d",dout_s);

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
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

仿真波形图：

在工程目录下的simulation/modelsim文件夹里便产生了sout和noise_out。

4. 使用matlab将产生的程序进行仿真验证

M程序：

%E4_7_NoiseAndCarrierOut.M
f1=200;       %信号1频率为200Hz
f2=800;       %信号2频率为800Hz
Fs=2000;      %采样频率为2KHz
N=12;         %量化位数

%从文本文件中读取数据
%测试输入数据分别放在Noise_in和S_in变量中
fid=fopen('C:\matlab work\fir1_1\filterCoe\noise.txt','r');
[Noise_in,N_n]=fscanf(fid,'%lg',inf);
fclose(fid);

fid=fopen('C:\matlab work\fir1_1\filterCoe\sin.txt','r');
[S_in,S_n]=fscanf(fid,'%lg',inf);
fclose(fid);

%滤波后的输出结果数据分别放在Noise_out和S_out变量中
fid=fopen('C:\matlab work\fir1_1\filterCoe\noise_out.txt','r');
[Noise_out,N_count]=fscanf(fid,'%lg',inf);
fclose(fid);

fid=fopen('C:\matlab work\fir1_1\filterCoe\sout.txt','r');
%fid=fopen('C:\matlab work\fir1_1\filterCoe\E4_7_Sout.txt','r');
[S_out,S_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(150:length(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=black_fpga;
m_hn=20*log10(abs(fft(hn,1024))); m_hn=m_hn-max(m_hn);

%设置幅频响应的横坐标单位为Hz
x_f=[0:(Fs/length(out_noise)):Fs/2];
%只显示正频率部分的幅频响应
mf_noise=out_noise(1:length(x_f));
mf_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_noise,'-',x_f,mf_hn,'--');
xlabel('频率(Hz)');ylabel('幅度(dB)');title('FPGA仿真白噪声信号滤波前后的频谱');
legend('输入信号频谱','输出信号频谱','滤波器响应');
grid;
subplot(212);
plot(x_f,mf_in_s,'--',x_f,mf_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; 
t_in_noise=Noise_in(1:length(t));
t_in_s=S_in(1:length(t));
t_out_noise=Noise_out(1:length(t));
t_out_s=S_out(1:length(t));
figure(2);
subplot(211);
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;
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
82

仿真波形图：


验证通过，设计成功。