一、利用MATLAB实现滤波

1.设计要求

（1）在MATLAB中，调用Filter Design & Analysis Tools工具，产生一个低通滤波器，以16bits定点量化滤波器系数，并产生XILINX的COE文件。

（2）在MATLAB中，产生两个正弦信号，数据位宽16bits，频率分别为3MHz和19MHz，数据长度为4096点。两个正弦信号合路（相加），送入该低通滤波器中，观察其输出信号。

2.任务一：设计滤波器

（1）在MATLAB中选择APP-滤波器设计工具，滤波器参数设置如下：

由于两路正弦波分别为3MHz和19MHz，所以将滤波器的通带截止频率设置为8MHz，将阻带开始频率设置为15MHz。

滤波器的采样频率应该大于信号最高频率的两倍，即38MHz，本文中取200MHz。注意：在Vivado中，FIR ip核的采样频率也要设置为200MHz，两处需要统一。

一般来说，数字滤波器的阶数越高，其复杂度越高，滤波器带内越平滑，具体取值取多少应在工程中取舍，本文中取15。

（2）点击左侧第三个小图标：设置量化参数，将滤波器算法指定为定点，将字长设置为16（由设计要求可知）。

点击输入/输出，可以看到该滤波器的输入范围是[-1,1]，输出小数长度是33。在后续输入时要注意信号的范围，如果不在要求的范围内会出错。

（3）点击任务栏目标-XILINX系数(.COE)文件，保存至相应目录下，后续在Vivoda中设置ip核时需要导入。

2.任务二

（1）代码

<code>clc

clear all

close all

%% 生成两个正弦信号

depth = 4096; % 存储器的单元数4096

widths = 16; % 数据宽度为16位

Fs = 200e6; % 采样频率为200MHz

t = 0:1/Fs:(4096-1)/Fs; % 时间向量

f1 = 3e6; % 第一个正弦信号频率为3MHz

f2 = 19e6; % 第二个正弦信号频率为19MHz

signal1 = sin(2*pi*f1*t)+1;

signal2 = sin(2*pi*f2*t)+1;

%% 合并两个信号

combined_signal = (signal1 + signal2)/4;

%% 调用低通滤波器设计函数

Hd = lowpass_filter;

%% 对合并信号进行滤波

filtered_signal = filter(Hd, combined_signal);

%% 绘制所有信号

figure;

subplot(4,1,1);

plot(t, double(signal1));

title('第一个正弦信号','Fontname','SimSun');

xlabel('时间 (秒)','Fontname','SimSun');

ylabel('幅度','Fontname','SimSun');

subplot(4,1,2);

plot(t, double(signal2));

title('第二个正弦信号','Fontname','SimSun');

xlabel('时间 (秒)','Fontname','SimSun');

ylabel('幅度','Fontname','SimSun');

subplot(4,1,3);

plot(t, double(combined_signal));

title('合并后信号','Fontname','SimSun');

xlabel('时间 (秒)','Fontname','SimSun');

ylabel('幅度','Fontname','SimSun');

subplot(4,1,4);

plot(t, double(filtered_signal));

title('滤波后信号','Fontname','SimSun');

xlabel('时间 (秒)','Fontname','SimSun');

ylabel('幅度','Fontname','SimSun');

上文中提到，滤波器的输入范围为[-1,1]，所以需要将合路信号变换到[-1,1]，这里选择先给两路信号加上幅值为1的直流分量,此时合路信号的范围为[0,4]，再除以4，即[0,1]，符合滤波器输入要求。这样做的原因在于：后续需要通过MATLAB产生正弦波数据表ROM，为了后续在Vivado中计算方便，产生时将信号扩大，并且添加了直流分量将信号搬移为正值，为了统一输入，此处做同样处理，即添加直流分量。

（2）结果

由图1-5可知，滤波后的信号频率与3MHz信号频率基本相同，但是幅值发生了变化，这是因为滤波器的幅度响应不是常数，所以3MHz的信号也有所衰减。

二、利用Vivado实现滤波

1.设计要求

（1）在Vivado软件中，采用FIR Compiler进行配置，并输入MATLAB产生的滤波器系数，产生一个低通滤波器。

（2）在Vivado软件中，产生一个正弦信号表：一个完整的正弦信号周期，共4096个点，精度为16bits，以ROM方式存放。

（3）例化两个正弦信号ROM，通过控制地址的步进量的方式，产生两个正弦信号输出，频率近似为3MHz和19MHz。

（4）两个正弦信号相加后，输入低通滤波器中，在Modelsim软件中，以HDL仿真方式观察其输出信号。

2.任务三：设计FIR

（1）IP核设置

①选择IP Catalog-FIR Compile

②在Filter Options中将Select Source设置为COE File，点击Cofficient File右侧的第一个框上传刚才MATLAB产生的.coe文件。

其他参数设置如下：

在MATLAB中设置的滤波器的采样频率为200MHz，在此要保持一致，即Input Sampling Frequency设置为200MHz。时钟频率Clock Frequency设置为200MHz，这里时钟频率和开发板保持一致即可，也应该和测试文件中的clk频率一致。输入数据位宽为16位，输出数据位宽为35位。

下图中输出数据为m_axis_data_tdata[39:0]，是40位，和Summary中Output Width不符。原因应该是这里设置成40位防止数据溢出，实际上是35位。

（2）代码

程序中认为输入的采样数据始终有效，因此将s_axis_data_tvalid永远置1。

<code>`timescale 1ns / 1ps

module LowBandPassFilter(

input clk,

input [15:0] wave_in,

output [39:0] wave_out,

output m_tvalid,

output s_tready

);

fir_compiler_0 fir1 (

.aclk(clk),

.s_axis_data_tvalid(1'b1),

.s_axis_data_tdata(wave_in),

.s_axis_data_tready(s_tready),

.m_axis_data_tvalid(m_tvalid),

.m_axis_data_tdata(wave_out)

);

endmodule

3.任务四：设计ROM

（1）MATLAB产生正弦波格式的COE文件

ROM的地址数为4096，数据位宽为16位，产生的信号幅值为10进制数。产生正弦信号表时，采用初始相位为0，周期为4096的波形，并且将正弦信号乘次方并乘（

2

15

2^{15}

215-1）再加

2

15

2^{15}

215，将数据变换到[1,

2

16

2^{16}

216-1]范围内。回忆上文中提到，使用MATLAB产生滤波器进行滤波时，也要先加上直流分量，把信号搬移到正值，上文就是为了与此处保持一致。MATLAB中滤波器的输入要求为[-1,1]，信号搬移后，需要除以4满足滤波器的输入要求，所以在Vivado中，合路信号也需要右移2位（即除4），后文会提到这一点。

depth = 4096; %存储器的单元数4096

widths = 16; %数据宽度为16位

fidc = fopen('sine.coe','wt'); %给文件起名字

fprintf(fidc , 'memory_initialization_radix=%d;\nmemory_initialization_vector=\n',10);

for x = 1 : depth-1

fprintf(fidc,'%d,\n',round((2^(widths-1)-1)*sin(2*pi*(x-1)/4096)+2^(widths-1)));

end

fprintf(fidc,'%d;\n',round((2^(widths-1)-1)*sin(2*pi*(depth-1)/4096)+2^(widths-1)));

fclose(fidc);

产生的COE文件内容如下：

memory_initialization_radix=10;

memory_initialization_vector=

32768,

32818,

32869,

32919,

32969,

33019,

33070,

33120,

33170,

······//中间省略

32466,

32517,

32567,

32617,

32667,

32718;//最后为分号

（2）在Vivado中设置ROM IP核

选择IP Catalog-Distributed Memory Generator

注：老师后来提醒我，Distributed Memory适合数据没那么多的时候，比如：16、32、64、128这样，1024、4096等数据规模较大的时候还是用Block Memory Generator比较好。

根据设计要求，ROM的深度为4096，数据位宽为16，即ROM的地址为12根，寻址范围0~4095，每一个地址上存储的数据为16位。

选择上一步中MATLAB生成的COE文件上传，Radix为10（与MATLAB生成时保持一致即可）。

（3）代码

<code>`timescale 1ns / 1ps

module rom(

input clk, //50M时钟

//input rst, //复位

output [15:0] data, //rom 输出的数据

output [11:0] addr

);

dist_mem_gen_0 Rom (

.a(addr), // input wire [11 : 0] a

.clk(clk), // input wire clk

.qspo(data) // output wire [15 : 0] qspo

);

endmodule

4.任务五：例化ROM产生不同频率的正弦信号

（1）思路

例化两个正弦信号ROM，通过控制地址的步进量的方式产生两个正弦信号输出，频率近似为3MHz和19MHz。即对ROM中的数据进行采样，采样频率不同，产生的信号频率也不同。其中，频率控制字的计算公式为：

M

=

f

o

u

t

∗

2

N

f

c

l

k

M=\frac{f_{out}*2^{N}}{f_{clk}}

M=fclk​fout​∗2N​（1）。其中，N为位宽，M为频率控制字，

f

o

u

t

{f_{out}}

fout​为输出频率，

f

c

l

k

{f_{clk}}

fclk​ 为时钟频率200MHz。产生3MHz的信号时，根据公式（1）可知，频率控制字M为12’d61，同理产生19MHz信号的频率控制字为12’d389。

（2）产生3MHz正弦波代码

产生3MHz正弦波代码如下：

`timescale 1ns / 1ps

module DDS_3M(

input clk,

input reset,

output [15:0] data_3M

);

parameter fre_word = 12'd61; //频率控制字 fre_word = f_out * 2^N / fclk N为累加器位宽

reg [11:0] addr_sin;

//相位累加器

always @(posedge clk or reset)

begin

if(!reset)

addr_sin <= 12'b0;

else

addr_sin <= addr_sin + fre_word;

end

wire [11:0]addra = addr_sin[11:0];

rom rom_3M (

.clk(clk), // input wire clk时钟

.addr(addra), // input wire [11 : 0] addra相位累加器输入给rom的地址

.data(data_3M) // output wire [15 : 0] douta从ROM返回的数据（3M正弦波的采样点）

);

endmodule

（3）产生19MHz正弦波代码

产生19MHz正弦波代码如下：

`timescale 1ns / 1ps

module DDS_19M(

input clk,

input reset,

output [15:0] data_19M

);

parameter fre_word = 12'd389; //频率控制字 fre_word = f_out * 2^N / fclk N为累加器位宽

reg [11:0] addr_sin;

//相位累加器

always @(posedge clk or reset)

begin

if(!reset)

addr_sin <= 12'b0;

else

addr_sin <= addr_sin + fre_word;

end

wire [11:0]addra = addr_sin[11:0];

rom rom_19M (

.clk(clk), // input wire clk时钟

.addr(addra), // input wire [11 : 0] addra相位累加器输入给rom的地址

.data(data_19M) // output wire [15 : 0] douta从ROM返回的数据（3M正弦波的采样点）

);

endmodule

（4）封装3MHz和19MHz正弦波产生代码

下面的代码让后续代码编写起来更简洁，即调用DDS时同时例化两路信号。

`timescale 1ns / 1ps

module DDS(

input clk,

input reset,

output [15:0] data_3M,

output [15:0] data_19M

);

DDS_3M DDS_3M_inst(

.clk(clk),

.reset(reset),

.data_3M(data_3M)

);

DDS_19M DDS_19M_inst(

.clk(clk),

.reset(reset),

.data_19M(data_19M)

);

endmodule

（5）结果

如图所示，wave_3M为频率为3MHz的正弦波，wave_19M为频率为19MHz的正弦波。

5.任务六：滤波

（1）testbench代码

在Vivado中，为了与MATLAB中保持一致，合路信号wave_sum[16:0] 需要除以4，即wave_sum右移两位，然后取[15:0]即可。

<code>`timescale 1ns / 1ps

module wave_simulation();

reg clk;

reg reset;

wire [15:0] wave_3M;

wire [15:0] wave_19M;

wire [15:0] wave_in;

wire s_tready;

wire m_tvalid;

wire [39:0] m_tdata;

reg [16:0] wave_sum;

wire[16:0] wave_inn;

always #2.5 clk = ~clk;

always @(posedge clk or reset)

if(!reset)

wave_sum = 'b0;

else begin

wave_sum = wave_3M + wave_19M;

end

assign wave_inn = wave_sum>>2;

assign wave_in = wave_inn[15:0];

initial begin

clk = 1'b0;

reset = 1'b0;

#1 reset = 1'b1;

end

DDS wave(

.clk(clk),

.reset(reset),

.data_3M(wave_3M),

.data_19M(wave_19M)

);

LowBandPassFilter LBF(

.clk(clk),

.s_tready(s_tready),

.wave_in(wave_in),

.wave_out(m_tdata),

.m_tvalid(m_tvalid)

);

endmodule

（2）结果分析

由下图可知，经过低通滤波器滤波后，输入的3M+19MHz信号（wave_in[15:0]）只剩下3MHz的单频信号。且当FIR滤波器输出有效时，m_tvalid信号置高。

三、结果分析

1.要求：

（1）误差分析

FFT的结果可以采用对数值分析，10log(x)（功率值）或20log(x) （幅度值），这样方便结合滤波器的滤波特性进行对比分析。

（2）时序分析

用250MHz、300MHZ、350MHZ和400MHz时钟进行分析，看设计的电路能够工作在哪个频率上。

（3）SNR分析

需要输入单音、多音冲击响应、多种随机信号进行分析。

2.误差分析

（1）思路

在Vivado中，两个频率的单音信号幅度都是一样的，合路信号经过滤波器后，两个频率的单音信号的幅度出现差异（我们期望的是3MHz的信号幅度没有衰减，而19MHz的信号幅度尽可能多地衰减），分析其真实差异值是否与MATLAB中滤波器的频率响应图上的两个频率的幅度差异一致。

（2）在MATLAB中，查看并绘制FIR滤波器的幅值响应

①在APP-滤波器设计工具（ 或者在命令行窗口直接输入filterDesigner）打开任务一中设计的滤波器，即可以看到滤波器的幅值响应。

②在图上单击鼠标右键，可以选择分析参数和采样频率。

③在分析参数中，可以选择归一化频率（因为MATLAB中绘制的滤波器的幅值响应横坐标也是归一化频率，为方便对比在此统一）。

**频率归一化的定义：**滤波器中各个频率值以截止频率作归一化后，频率都是截止频率的相对值，没有了量纲。信号处理工具箱中经常使用的是奈奎斯特频率（也叫做分析频率、截止频率），它被定义为采样频率的二分之一，在滤波器的阶数选择和设计中的截止频率均使用奈奎斯特频率进行归一化处理。

实际物理频率表示AD采集物理信号的频率，fs为采样频率，由奈奎斯特采样定理可以知道，fs必须大于等于信号最高频率的2倍才不会发生信号混叠，因此fs能采样到的信号最高频率为fs/2。角频率是物理频率的2π倍，这个也称模拟频率。归一化频率是将实际物理频率除以分析频率，分析频率为采样频率的二分之一。

在本设计中，采样频率为200MHz，归一化频率即实际物理频率/分析频率，即际物理频率/100MHz。如果将归一化频率转换为角频率，则将归一化频率乘以2*pi。如果将归一化频率转换为Hz，则将归一化频率乘以分析频率。

④我们设计的滤波器的幅值响应如下图所示，可以单击图上的点建立游标数据，方便后续我们对比。

⑤通过代码绘制滤波器的幅值响应

添加的代码是：

<code>%% 绘制滤波器的幅频响应

figure;

freqz(b);

更改后的低通滤波器的代码如下：

function Hd = lowpass_filter

% All frequency values are in MHz.

Fs = 200; % Sampling Frequency

N = 15; % Order

Fpass = 8; % Passband Frequency

Fstop = 15; % Stopband Frequency

Wpass = 1; % Passband Weight

Wstop = 1; % Stopband Weight

% Calculate the coefficients using the FIRLS function.

b = firls(N, [0 Fpass Fstop Fs/2]/(Fs/2), [1 1 0 0], [Wpass Wstop]);

Hd = dfilt.dffir(b);

%% 绘制滤波器的幅频响应

figure;

freqz(b);

% Set the arithmetic property.

set(Hd, 'Arithmetic', 'fixed', ...

'CoeffWordLength', 16, ...

'CoeffAutoScale', true, ...

'Signed', true, ...

'InputWordLength', 16, ...

'inputFracLength', 15, ...

'FilterInternals', 'FullPrecision');

denormalize(Hd);

得到的幅值响应如下图所示：

（3）从Vivado中导入滤波前后信号的采样值

①在testbench文件中写入滤波后信号（m_tdata[39:0]）的值。

需要添加的代码如下：

<code>

reg [11:0] cnt;

integer file_sig_before_fir;

integer file_sig_after_fir;

always @(posedge clk) begin

if(m_tvalid) begin

if(cnt <= 12'd4095) begin

cnt <= cnt + 1'b1;

$fdisplay(file_sig_before_fir,"%d\n",wave_in);

$fdisplay(file_sig_after_fir,"%d\n",m_tdata);

end

else begin

$fclose(file_sig_after_fir);

$fclose(file_sig_before_fir);

end

end

else

;

end

initial begin

file_sig_before_fir = $fopen("This is your path","w");

file_sig_after_fir = $fopen("This is your path","w");

cnt = 12'b0;

clk = 1'b0;

reset = 1'b0;

#1 reset = 1'b1;

end

这部分引入了一个12位的cnt计数，是因为如果以m_tvalid是否置高（当FIR滤波器输出有效时，m_tvalid一致信号置高）判断是否写入的话，文件将无法关闭，所以写入的文件一直是空白的，只有执行了$fclose语句后,文件才会被写入并保存。但是这里共写入了8192个数据（4096*2=8196），这里我还没有想明白。

Testbench完整代码如下：

`timescale 1ns / 1ps

module wave_simulation();

reg clk;

reg reset;

wire [15:0] wave_3M;

wire [15:0] wave_19M;

wire [15:0] wave_in;

wire s_tready;

wire m_tvalid;

wire [39:0] m_tdata;

reg [16:0] wave_sum;

wire[16:0] wave_inn;

reg [11:0] cnt;

integer file_sig_before_fir;

integer file_sig_after_fir;

always #2.5 clk = ~clk;

always @(posedge clk or reset)

if(!reset)

wave_sum = 'b0;

else begin

wave_sum = wave_3M + wave_19M;

end

always @(posedge clk) begin

if(m_tvalid) begin

if(cnt <= 12'd4095) begin

cnt <= cnt + 1'b1;

$fdisplay(file_sig_before_fir,"%d\n",wave_in);

$fdisplay(file_sig_after_fir,"%d\n",m_tdata);

end

else begin

$fclose(file_sig_after_fir);

$fclose(file_sig_before_fir);

end

end

else

;

end

assign wave_inn = wave_sum>>2;

assign wave_in = wave_inn[15:0];

initial begin

file_sig_before_fir = $fopen("This is your path","w");

file_sig_after_fir = $fopen("This is your path","w");

cnt = 12'b0;

clk = 1'b0;

reset = 1'b0;

#1 reset = 1'b1;

end

这部分并不困难，但是我当时做了很久，是因为文件写入一直是空白，弄了很久才发现是/写成\，我是现在本地建了文档，然后直接复制文档路径，系统中默认的路径是\，读者可以注意这个问题。

（4）在MATLAB中分析滤波前后信号的幅度变化

参考代码如下：

clc

clear

close all

Fs = 200e6;

t = 0:1/Fs:(4096-1)/Fs; % 时间向量

f1 = 3e6; % 第一个正弦信号频率为3MHz

f2 = 19e6; % 第二个正弦信号频率为19MHz

depth = 4096; %存储器的单元数4096

widths = 16; %数据宽度为16位

sig_3M = [];

sig_19M = [];

decay_new = [];

%% 读Vivado中滤波前的数据

fid=fopen('This is your path');

VivadoBeforeFir=fscanf(fid,'%f');

fclose(fid);

VivadoBeforeFir = VivadoBeforeFir/2^15; % 输入小数为15位

BeforeFir = VivadoBeforeFir(2:4097,:);

%% 读Vivado中滤波后的数据

fid=fopen('This is your path');

VivadoAfterFir=fscanf(fid,'%f');

fclose(fid);

VivadoAfterFir = VivadoAfterFir/2^33; % 输出小数为33位

AfterFir = VivadoAfterFir(2:4097,:);

%% 绘制滤波前后信号频谱

figure;

subplot(2,1,1);

[~,MagBeforeFir] = DrawFFT(BeforeFir,Fs);

title('滤波前合路信号频谱','Fontname','SimSun');

xlabel('归一化频率(Hz)','Fontname','SimSun');

ylabel('幅度(dB)','Fontname','SimSun');

subplot(2,1,2);

[~,MagAfterFir,f] = DrawFFT(AfterFir,Fs);

title('滤波后合路信号频谱','Fontname','SimSun');

xlabel('归一化频率(Hz)','Fontname','SimSun');

ylabel('幅度(dB)','Fontname','SimSun');

%% 绘制滤波前后信号波形

figure;

subplot(2,1,1);

plot(t,BeforeFir(1:4096,1));

title('滤波前合路信号波形','Fontname','SimSun');

xlabel('时间 (秒)','Fontname','SimSun');

ylabel('幅度','Fontname','SimSun');

subplot(2,1,2);

plot(t,AfterFir(1:4096,1));

title('滤波后信号波形','Fontname','SimSun');

xlabel('时间 (秒)','Fontname','SimSun');

ylabel('幅度','Fontname','SimSun');

%% 幅值响应

decay = 20*log10(MagAfterFir./MagBeforeFir);

%% 绘制滤波器的幅值响应

figure;

plot(f/(Fs/2), decay);

title('Vivado滤波幅度响应','Fontname','SimSun');

xlabel('归一化频率(Hz)','Fontname','SimSun');

ylabel('幅度(dB)','Fontname','SimSun');

（5）结果分析

①滤波前后信号频谱

由下图可知，3MHz信号滤波前后衰减为20log10（0.263612/0.249992）=0.46078dB，同理，19MHz信号滤波前后衰减为20log10（0.00262554/0.249993）= -39.57418dB。

②滤波前后信号波形

③Vivado滤波幅值响应

由游标数据知，3MHz信号的幅值响应为0.460778，19MHz信号的幅值响应为-39.5742，这与上文MATLAB中设计的滤波器的幅值响应一致。不过这里画出来的幅值响应图和MATLAB中FilterDesigner得到的幅值响应不太一样（仅仅0Hz、3MHz、19MHz）的幅值衰减一致，原因我还没有弄清楚。