当前位置：网站首页>基于FPGA的FIR滤波器的实现（4）— 串行结构FIR滤波器的FPGA代码实现
基于FPGA的FIR滤波器的实现（4）— 串行结构FIR滤波器的FPGA代码实现
2022-08-11 05:53:00 【南邮学渣】
文章目录
前言
前面已经详尽的介绍了有关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;
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; %绘图
频率-幅度响应图
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);
频率-幅度响应图
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
仿真测试模块：
`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
仿真波形图：
在工程目录下的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;
仿真波形图：
验证通过，设计成功。
边栏推荐
- Daily sql-statistics of the number of professionals (including the number of professionals is 0)
- 亚马逊获得AMAZON商品详情 API 返回值说明
- daily sql - user retention rate for two days
- sql--Users who have purchased more than 3 times (inclusive) within 7 days (including the current day), and the purchase amount in the past 7 days exceeds 1,000
- 淘宝API常用接口与获取方式
- Find the shops that have sold more than 1,000 yuan per day for more than 30 consecutive days in the past six months
- ROS 话题通信理论模型
- daily sql - query for managers and elections with at least 5 subordinates
- mysql视图与索引
- ssh服务攻防与加固
猜你喜欢
HCIP OSPF dynamic routing protocol
ssh服务攻防与加固
拼多多API接口（附上我的可用API）
Daily sql-employee bonus filtering and answer rate ranking first
抖音获取douyin分享口令url API 返回值说明
HCIP WPN experiment
Get Pinduoduo product information operation details
李沐d2l(十)--卷积层Ⅰ
JD.com product details API call example explanation
详解BLEU的原理和计算
随机推荐
pytorch，numpy两种方法实现nms类间+类内
亚马逊获得AMAZON商品详情 API 返回值说明
HCIP BGP built adjacent experiment
详解BLEU的原理和计算
《Show and Tell: A Neural Image Caption Generator》论文解读
JD.com product details API call example explanation
那些事情是用Unity开发项目应该一开始规划好的？如何避免后期酿成巨坑？
HCIP MGRE\OSPF Comprehensive Experiment
抖音API接口
exness:黄金1800关口遇阻，静待美国CPI出炉
每日sql - 判断+聚合
Taobao API common interface and acquisition method
从 dpdk-20.11 移植 intel E810 百 G 网卡 pmd 驱动到 dpdk-16.04 中
HCIA experiment
淘宝sku API 接口（PHP示例）
torch.cat()用法
numpy和tensor增加或删除一个维度
抖音API接口
你是如何做好Unity项目性能优化的
已解决EROR 1064 (42000): You have an error in. your SOL syntax. check the manual that corresponds to yo