FIR滤波器的MATLAB与FPGA的设计实现
引言
数字滤波器是数字信号系统里面最常用、最重要、最基本的元件之一。有限脉冲响应(Finite Impulse Response,FIR)滤波器可以设计成任意幅频特性，同时保证精确，严格的线性相位特性。因此在电子通信、图像处理、模式识别等领域，FIR数字滤波器被广泛应用。本文主要通过以低通滤波器的MATLAB与FPGA的设计实现来介绍FIR滤波器的原理介绍，架构研究，设计实现等。
原理介绍
FIR滤波器的方程表达式为： y(n)=∑x(k)h(n-k) k=0,1…N-1;0≦n≦N-1;

其方程分解式为：y(n)=h(0)x(n)+h(1)x(n-1)+h(2)x(n-2)+…+h(k)x(n-k);

根据FIR滤波器的幅频特性可知，FIR滤波器为奇对称和偶对称两种形式。

例如当滤波器阶数为偶数n=3时，根据线性相位跟对称法则可得：h(0)=h(3)，h(1)=h(2);

可得到方程式为：y(n)=h(0){x(n-1)+x(0)}+h(1){x(n-2)+x(1)}；

当滤波器阶数为奇数n=4时，根据线性相位跟对称法则可得：h(0)=h(4)，h(1)=h(3)，h(2)独立一个;

可得到方程式为：y(n)=h(0){x(n-1)+x(0)}+h(1){x(n-2)+x(n-3)}+h(2)x(n-2)；

故其系统函数为



从其系统函数可以看得出FIR系统是线性时不变系统，只在原点上存在极点，这使得FIR系统具有全局稳定性。从FIR系统函数也可以看得出FIR滤波器其实是由一个抽头延迟线加法器和乘法器的集合构成的。每个乘法器的操作系数为FIR滤波器系数。因此这种抽头延迟线结构的架构实现图可如下所示：


MATLAB设计实现
在matlab里面设计数字滤波器既可以通过窗函数自己编写程序来设计FIR滤波器，也可以通过调用工具箱来设计FIR滤波器。通过编写程序来设计FIR滤波器比较普遍，滤波器的参数等改动也比较方便，与FPGA做联合功能仿真验证也方便。而直接调用工具箱就直接了当，较节省设计时间。本文将对两种方法进行介绍，不过本文设计采用的是凯塞窗函数来设计FIR滤波器。
窗函数设计法
利用窗函数法来设计FIR滤波器是设计FIR滤波器最简单，最普遍的方法，在MATLAB里面窗函数有海明（Hamming）窗，汉宁(Hanning)窗，凯塞(Kaiser)窗等，本文设计采用的是凯塞窗函数来设计FIR滤波器。通过最优算法来计算出过渡带，纹波等参数供给凯塞窗来计数出β值，以及滤波器的阶数。

凯塞窗函数的语法形式为：[n,wn,beta,ftype]=kaiserord(fc,mag,dev,fs)

其中各参数的意义介绍如下所述：

（1）fc及fs：fc是指过渡带的起始点和终点，fs指的是采样评率；

（2）mag：f指定了过渡带，向量a用于指定这些频率段的理想幅度值,例如a=[1 0]可以代指低通滤波为1，高频截止为0；

（3）dev：用于指定通带或阻带内的容许误差；

（4）n,wn：n是返回凯塞窗函数满足设计的最小阶数，wn是返回滤波器的截止频率点；

（5）beta：为凯塞窗函数计算得到的β值。

（6）ftype：返回设计滤波器的类型,low,hige等;

通过凯塞窗函数最终可以得到滤波器的最佳设计阶数n为多少,然后根据n来设计FIR滤波器，本文采用的是最优滤波器设计法，利用firpm函数根据最大误差最小准则来设计滤波器。

最优滤波器firpm函数的语法形式为：b_firpm = firpm(n,fpm,magpm)

其中各参数的意义介绍如下所述：

（1）n：为滤波器的设计阶数；

（2）fpm，magpm：指定滤波器的设计频率段对应的幅度值magpm；

（3）b_firpm：返回设计滤波器的系数h(n)。

到此为止设计FIR滤波器的matlab设计实现如下图所示，本文设计的是一个低通滤波器。



现在滤波器的系数b_firpm是可以求出来，但是注意是小数，在FPGA中实现还是比较麻烦还占用大量资源，所以在这里还进行了量化处理，对matlab数字滤波器的b_firpm进行14bit的量化处理（我在做量化的时候做了10，12，14bit的设计，发现14bit的滤波器效果最好）。



运行函数会发现生成了28阶的FIR滤波器系数如下所示，当n=27,滤波器阶数为28的时候，发现滤波器系数为对称相等的h(14)=h(15)=8191;

可验证到与前面的原理介绍是一致的。这样子后面在FPGA设计滤波器的时候可以用线性相位结构来设计。



运行函数，设计的低通滤波器性能如下图所示：



接下来生成高斯白噪声通过MATLAB设计的低通滤波器，来验证我们设计的低通滤波器性能如何，同时将生成的高斯白噪声生成数据文件txt供给后面的FPGA设计的低通滤波器做modelsim的功能仿真测试。生成白噪声的程序如下图所示，q_noise为高斯白噪声，filter_noise为滤波后噪声。



将生成的高斯白噪声生成数据文件txt供给后面的FPGA设计的低通滤波器做modelsim的功能仿真测试。



调用自己MATLAB设计的低通滤波器对高斯白噪声进行滤波，程序设计如下图所示：



运行程序看设计的低通滤波器的性能如何，如下图所示，matlab设计的低通滤波器很完美的把高频率的数据进行截止滤波掉。在低频率的时候保持原高斯白噪声的输出状态。


tool工具箱设计法
在matlab的应用程序里面有一个filter design & analysis tool点进去即为快速设计滤波器的工具。操作界面如下图所示，较为简单不做详细介绍。



点analysis里面的filter coefficient即可出现设计滤波器的FIR滤波器系数，通过file-export可以把滤波器系数导出到workplace。


FPGA设计实现
在前面的原理介绍中可以得知，根据阶数的系数为对称相等的，于是抽头延迟线结构可以升级为线性相位优化结构来设计，具体我的架构设计图如下所示：



从上面的架构设计图可以看得出，FIR的设计大概可以分为三步分，数据流的移位延迟模块，系数和乘法器模块，累加和模块三大部分。

整体的端口设计介绍如下图所示：


数据流的移位延迟模块
data_mem为28个14位宽的寄存器，进行数据流的采集跟移位保存。其设计如下图所示：


累加和的设计模块
这里的数据流有符号正负，所以在运算时为了不溢出，扩展一位，并用signed声明进行有符号的运算设计。


系数乘法器的设计模块
乘法器的设计端口介绍如下图所示：



乘法器的设计原理介绍如下图所示，具体的乘法设计可参考之前的博客章文：基于FPGA的乘法器原理介绍及设计实现，具体博文链接为：

https://blog.csdn.net/weixin_39015789/article/details/101430566



最后再将所有的乘积和进行累加便是对数据流进行低通滤波后的结果了。
MATLAB与FPGA的连仿验证
到现在为止，低通滤波器的FPGA设计实现也已经完毕了，接下来就是在modelsim里面读取matlab产生的高斯白噪声数据流，供给FPGA设计的FIR滤波器进行滤波计算，看两者的结果是否一模一样。

FPGA的textbench设计验证开始
读取高斯白噪声数据流
从外部txt文件(noise_in.txt)读入matlab产生的高斯白噪声数据流作为测试激励输入数据流。具体的设计如下所示：


将FPGA滤波后的数据流读取到文件
将FPGA滤波后的数据流保存读取到文件，然后提供给matlab进行傅里叶变换，看FPGA设计的低通滤波器性能是如何的。


MATLAB与FPGA的连仿
读取FPGA滤波后的数据，进行傅里叶变换，求高斯白噪声及滤波后信号的幅频响应，对比滤波前后的信号，查看滤波器的性能如何。



通过MATLAB将FPGA的滤波数据进行傅里叶变换的频谱滤波效果图如下所示，可以看到FPGA设计的低通滤波器在低频域的时候完美将高斯白噪声进行输出，在高频率的时候，对高斯白噪声进行滤波处理。效果性能与前面MATLAB设计的滤波器效果性能很是一致，滤波器的效果性能都很好。

用matlab的fdatool设计滤波器，导出参数，然后在FPGA中设计滤波器的方法。
1、设计fir插值滤波
2、设计半带滤波器
3、CIC滤波器
4、频谱相关
111111111111111111111111111111111111111111111111111
1、设计fir插值滤波
MATLAB 中fir插值滤波
MATLB产生系数的方法：在命令行中输入fdatool







注意：creat a multirate filter 设置中，其实设计的是一个多相滤波器。不管怎么更改这里的插值倍数和采样率，其实滤波器的系数是不改变的。原因是，这里的滤波器参数只是多相滤波器中原型滤波器的系数，
所以，滤波器系数不变。但是，为什么滤波器系数不变，反而在调用改滤波器的时候还能实现插值呢？这是因为，你已经设计了该滤波器是多相插值滤波器，在数据处理中，实际上已经按多相插值来操作了，
多相插值的结构本来就可以实现插值/抽取。
至此，整个滤波器就设计完成了。需要生成滤波器的FIR文件和滤波系数。


ISE 中fir插值滤波
采用FIR滤波器中第一个版本的


 说明：这里输出精度的截尾其实是有技巧的。如果是全精度输出，则高几位经常都是符号位。但是怎么截取才能不导致有效数据丢失呢？以前的做法是，把数据放在MATLAB里找一下最大值，然后，通过最大值来设计位宽。后来，周弘告诉我，可以直接用全精度输出的位数 - 系数扩大的位数即可。上面的 例子，全精度输出36位，则截尾为36 - 18 = 18即可满足要求。
FIR滤波器输出结果，应该为有符号数。将该数据和matlab中的数据对比，高几位差别不是很大，说明正确。 需要注意的是，FIR滤波器的前几位数据是无效数据，不弄用于和matlab对比。通常，在N/2(N为滤波器的阶数)个时钟延迟之后的数据才是真正的FIR输出，可以用于和MATLBA 对比。

2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222
2、设计半带滤波器
MATLAB 中半带滤波器




1是由2来自动确定阶数，也可以通过kaiser窗的公式，根据要求来计算滤波器阶数。现在用的方法是，如上图所示，不停的更改Apass的值，看频谱图。
还有一种是确定阶数，设置采样频率，截止频率自动为采样频率的一半。
注意，这个界面设计完成之后，还要进行2倍插值的配置。否则，这只是一个简单的低通滤波器，而不是插值滤波器。
ISE 中半带滤波器

1、首先，你导入的COE文件要是半带滤波器的格式。
2、如果你导入的文件是半带滤波器的格式，那么第二个配置界面会有如下变换.注意，差值倍数设置为2




FPGA直接调用IP核就可以用这个滤波器了，那MATLAB怎么用呢？
1、用滤波器系数，滤波器的本质就是卷积。所以：Ipulse= conv(Izero,ht); %ht就是滤波器系数，和输入数据卷积即可。----这种方法只适用于没有插值/抽取的滤波器
2、调用matlab生成的M文件。如果该滤波器中有插值或者抽取，就不能用方法1.因为在配置插值的那个界面的时候，其滤波器系数是不变的。只有用下面的方法调用滤波器的时候才能实现插值或者滤波
%---------------------------这个比较麻烦，下面的简单-----------------
% N= 32; %%滤波器阶数

% da = Izero (: , N/2 +1 : end - N/2); %%是为了和xilinx中的数据对应。

% Hd = fir_cos; %%FIR滤波器，在fdatool中产生

% d2= filter ( Hd , [da , zeros(1, N/2) ]); %%将数据通过FIR滤波器。

%---------------------------下面的简单-----------------
Hd = fir_cos; %%FIR滤波器，在fdatool中产生
d2 = filter (Hd,da ); %da是需要滤波器的数据。

如果只知道滤波器的系数，怎么用系数来操作呢？
I_DATA = filter(fir_c,1,I_Symbols);

MATLAB产生的系数基本上都是小数，且位数特别宽。可以尝试进行归一化为整数之后再用。归一化的代码为（以归一化16位为例）：
ht_16 = round (ht/max(abs(ht))*(2^15-1));
在真正使用的时候发现，将系数按照上面方式归一化之后，FPGA结果和MATALB结果对不上。最后还是放弃了这种方法。还是按照上面将的那种系数带小数点的方法好。不过，输入输出数据要对应经过一定的放大或缩小之后 FPGA结果和MATALB结果 才能对的上。


3、CIC滤波器
在MATLAB中使用CIC滤波器的方法很简单。打开fdatool,直接点击这个:


至此,设计完毕.
用这种方法设计好滤波器之后，发现用psd不能画频谱。把滤波之后的数据 double 一下再用就可以了
在FPGA中,直接点击插值滤波器即可.

FIR滤波器的matlab设计

1.采样fir1函数设计

b= fir1(n,wn),

b= fir1(n,wn,’ftype’)

b= fir1(n,wn,’ftype’,window)

b= fir1(…,’noscale’)

函数参数说明

b:返回的FIR滤波器单位脉冲响应，脉冲响应为偶对称，长度为n+1

n:滤波器的阶数，需要注意的是，设计出的滤波器长度为n+1.

wn:滤波器的截至频率，；需要注意的是，wn的取值范围为0< wn<1, 1对应为信号采样频率的1/2；如果wn是个单个数值，且ftype参数为low,则表示设计截止频率为wn的低通滤波器，如ftype,参数为high,则表示设计截止频率为wn的高通滤波器，如果wn是由两个数组成的向量[wn1,wn2],ftype为stop,则表示设计带阻滤波器，ftype为bandpass，则表示设计带通滤波器；如果wn是由多个数组成的向量，则表示根据ftype的值设计多个带通或带阻的滤波器；ftype为DC-1，表示设计的第一个频带为通带，ftype为DC-0，表示设计的第一个频带为阻带。

window:指定使用的窗函数向量，默认为海明窗（Hamming）,最常用的窗函数有汉宁窗（Hanning）、海明（Hamming）窗、布拉克曼（Blackman）窗和凯塞（kaiser）窗，可以在MATLAB界面中输入help window命令查询各种窗函数名称。

noscale: 指定是否归一化滤波器的幅度。

2.函数使用方法

例：设计一个归一化截止频率为0.2，阶数为11、采用海明窗的低通滤波器

b= fir1(11,0.2);

plot(20log(abs(fft(b)))/log(10))

例： 分别设计长度为41的低通（截止频率为200HZ）、高通（截止频率为200HZ）、带通（通带为200400HZ）、带阻滤波器（阻带为200400HZ）FIR滤波器采样频率为2000HZ

N = 41;                                   %滤波器长度

fs =2000;                              %采样频率

%各种滤波器的特征频率

fc_lpf =200;

fc_hpf =200;

fp_bandpass= [200,400];

fc_stop =[200,400];

%以采样频率的一半，对频率进行归一化处理

wn_lpf =fc_lpf2/fs;

wn_hpf =fc_hpf2/fs;

wn_bandpass= fp_bandpass2/fs;

wn_stop =fc_stop2/fs;

%采样fir1函数设计FIR滤波器

b_lpf =fir1(N-1,wn_lpf);

b_hpf =fir1(N-1,wn_hpf,‘high’);

b_bandpass= fir1(N-1,wn_bandpass,‘bandpass’);

b_stop =fir1(N-1,wn_stop,‘stop’);

%求滤波器的幅频响应

m_lpf =20log(abs(fft(b_lpf)))/log(10);

m_hpf =20log(abs(fft(b_hpf)))/log(10);

m_bandpass=20log(abs(fft(b_bandpass)))/log(10);

m_stop =20log(abs(fft(b_stop)))/log(10);

%设置幅频响应的横坐标单位为HZ

x_f=[0:(fs/length(m_lpf)):fs/2];

%绘制单位脉冲响应

subplot(421);stem(b_lpf);xlabel(‘n’);ylabel(‘h(n)’);

subplot(423);stem(b_hpf);xlabel(‘n’);ylabel(‘h(n)’);

subplot(425);stem(b_bandpass);xlabel(‘n’);ylabel(‘h(n)’);

subplot(427);stem(b_stop);xlabel(‘n’);ylabel(‘h(n)’);

%绘制幅频响应曲线

subplot(422);plot(x_f,m_lpf(1:length(x_f)));xlabel(‘频率(HZ)’,‘fontsize’,8);ylabel(‘幅度(dB)’,‘fontsize’,8);

subplot(424);plot(x_f,m_hpf(1:length(x_f)));xlabel(‘频率(HZ)’,‘fontsize’,8);ylabel(‘幅度(dB)’,‘fontsize’,8);

subplot(426);plot(x_f,m_bandpass(1:length(x_f)));xlabel(‘频率(HZ)’,‘fontsize’,8);ylabel(‘幅度(dB)’,‘fontsize’,8);

subplot(428);plot(x_f,m_stop(1:length(x_f)));xlabel(‘频率(HZ)’,‘fontsize’,8);ylabel(‘幅度(dB)’,‘fontsize’,8);

例：采样各种窗函数，利用MATLAB软件分别设计截止频率为200HZ、采样频率为2000HZ的FIR低通滤波器，滤波器长度为81

N = 81;                         %滤波器长度

fs =2000;                    %采样频率

fc = 200;                       %低通滤波器的截止频率

%生成各种窗含数

w_rect =rectwin(N)’;

w_hann =hann(N)’;

w_hamm =hamming(N)’;

w_blac =blackman(N)’;

w_kais =kaiser(N,7.856)’;

%采用fir1函数设计FIR滤波器

b_rect =fir1(N -1 ,fc2/fs,w_rect);

b_hann =fir1(N -1 ,fc2/fs, w_hann);

b_hamm =fir1(N -1 ,fc2/fs, w_hamm);

b_blac =fir1(N -1 ,fc2/fs, w_blac);

b_kais =fir1(N -1 ,fc2/fs, w_kais);

%求滤波器的幅频响应

m_rect =20log(abs(fft(b_rect,512)))/log(10);

m_hann =20log(abs(fft(b_hann,512)))/log(10);

m_hamm =20log(abs(fft(b_hamm,512)))/log(10);

m_blac =20log(abs(fft(b_blac,512)))/log(10);

m_kais =20log(abs(fft(b_kais,512)))/log(10);

%设置幅频响应的坐标单位为HZ

x_f =[0:(fs/length(m_rect)):fs/2];

%只显示正频率部分的幅频响应

m1 =m_rect(1:length(x_f));

m2 =m_hann (1:length(x_f));

m3 =m_hamm (1:length(x_f));

m4 =m_blac (1:length(x_f));

m5 =m_kais (1:length(x_f));

%绘制幅频响应曲线

plot(x_f,m1,’.’,x_f,m2,’’,x_f,m3,‘x’,x_f,m4,’–’,x_f,m5,’-’);

xlabel(‘频率(Hz)’,‘fontsize’,8);ylabel(‘幅度(dB)’,‘fontsize’,8);

legend(‘矩形窗’,‘汉宁窗’,‘海明窗’,‘布拉克窗’,‘凯塞窗’);

grid;



二、采用kaiserord函数设计

语法形式

[n,wn,beta,filtype]=kaiserord(f,a,dev,fs)

参数说明

f及fs:如果f是一个向量，其中的元素为待设计滤波器的过渡带的起始点和结束点；

如果没有fs参数，f中元素的取值范围为0~ 1，即相对于采样频率一半的归一化频率；有fs参数，则fs为信号采样频率，f中元素即为实际的截至频率。例如需要设计滤波器的过滤带宽为1000~ 1200Hz、2000~2100Hz;信号采样频率为8000Hz；如没有设置fs参数，则f=[0.25 0.3 0.5 0.525];如设置fs为8000，则f=[1000 1200 2000 2100]。

a:a是一个向量，参数f确定了待设计滤波器的过渡带，向量a用于指定这些频率段的理想幅度值；a的第一个参数a1对应为f中的0f1频段，第二个参数a2对应f中的f2f3频段，后续对应关系一次类推。比如，对于上面讲述f及fs参数的例子，设置a=[1 0 1]，则表示需要设计的带阻滤波器。可以看出，由f及a可以表示滤波器的类型。

dev: dev是一个向量，用于指定通带或阻带内的容许误差。同样的是上述例子，需求通带容许误差为0.01，阻带容许误差为0.02，则dev=[0.01 0.02 0.01]。

n: 返回值n为kaiserord函数根据滤波器要求，得到满足设计的最小阶数。

wn:返回值wn是一个向量，kaiserord函数计算得到的滤波器截止频率点。

beta:返回值beta是根据滤波器要求，kaiserord函数计算得到的β值。

ftype:返回值ftype是根据设计要求获得的滤波器类型参数。

三、采用fir2函数设计

任意响应滤波器：滤波器的幅度频率响应在指定的频段范围内有不同的幅值，如在00.1的理想幅值为一，在0.20.4频段内的幅值为0.5，在0.6~0.7频段的幅值为1等。fir2H函数的算法是：根据要求的幅频响应的向量形式进行插值，而后进行傅里叶变换得到理想滤波器的单位脉冲响应，最后利用窗函数对理想滤波器的单位脉冲响应进行截短处理，fir2函数的6种语法形式如下。

b=fir2(n,f,m);

b=fir2(n,f,m,window);

b=fir2(n,f,m,npt);

b=fir2(n,f,m,npt,window);

b=fir2(n,f,m,npt,lap);

b=fir2(n,f,m, npt,lap ,window);

n及b：滤波器的阶数，与fir1函数类似，返回值b为滤波器系数，其长度为n+1。

同时，根据FIR滤波器的结构特点，当设计的滤波器在归一化频率为1处的幅度值不为0时，n不能为奇数。

f及m: f是一个向量，取值为01之间。对于为滤波器的归一化频率；m是长度与f相同的向量，用于设置对应频率范围内的理想幅值。例如，要求设计的滤波器在00.125范围内幅度为1,在0.1250.25范围内幅度为0.5，在0.250.5范围内幅度为0.25，在0.5~1范围内幅度为0.125，则f可以表示为[0 0.125 0.125 0.25 0.25 0.5 0.5 1],m可以表示为[1 1 0.5 0.5 0.25 0.25 0.125 0.125]。

window: window是一个向量，用于指定窗函数的种类，其长度为滤波器长度n+1,当没有指定窗函数时，默认为海明窗。

npt: npt 是一个正整数，用于指定在对幅度响应进行插值时的插值点个数，其默认值为512。

lap：lap是一个正整数，用于指定对幅度响应进行插值时，对于不连续点转变成连续时的点数，其默认值为25

例：利用fir2函数设计120阶的FIR滤波器，要求设计的滤波器在归一化频率0~ 0.125范围内幅度为1，在0.125~ 0.25范围内幅度为0.5，在0.25~ 0.5范围内幅度为0.25，在0.5~1范围内幅度为0.125.

N=120;                                                       %滤波器阶数

fc=[0 0.125 0.125 0.25 0.25 0.5 0.5 1];     %截至频率

mag =[1 1 0.5 0.5 0.25 0.25 0.125 0.125];%理想滤波器幅度

b=fir2(N,fc,mag);                                        %设计海明窗滤波器

freqz(b);                                                      %绘制频率响应曲线

三、采用firpm函数设计

firpm函数的语法形式

b=firpm(n,f,a);

b=firpm(n,f,a,w);

b=firpm(n,f,a,’ftype’);

b=firpm(n,f,a,w,’ftype’);

[b,delta]=firpm(…)

参数说明

n及b:滤波器的阶数,与fir1函数类似，返回值b为滤波器系数，其长度为n+1.

f及a:f是一个向量，取值为0~1之间，对应为滤波器的归一化频率；a是长度与f相同的向量，用于设置对应频段范围内的理想幅值。

w：w是长度为f的1/2的向量，表示设计滤波器时，实现对应频率幅度值得权值。

用下标表示向量的元素，则w0对应的是f0f1频段，w1对应的是f2f3频段，以此类推。权值越高，则实现时对应频段的幅值越接近理想状态。

ftype:ftype用于指定滤波器的结构类型，如没有设置参数，则表示设计偶对称脉冲响应的滤波器；如设置为“hilbert”,则表示设计奇对称结构的滤波器，即具有90°相移特性；如设置为“differentiator”,则表示设计奇对称结构的滤波器，且设计时针对非零幅度的频带进行了加权处理，使滤波器的频带越低，则幅值误差越小。

delta:delta为返回的滤波器的最大纹波值。

例：利用凯塞窗函数设计一个低通FIR滤波器，过渡带为1000~1500Hz,采样频率为8000Hz,通带纹波最大为0.01，阻带纹波最大为0.05，利用海明窗及firpm函数设计相同的低通滤波器，截至频率为1500Hz，滤波器阶数为凯塞窗函数求取的值。绘出三种方法设计的幅度频率的响应曲线

fs = 8000;                            %采样频率

fc=[1000 1500 ];                %过渡带

mag= [1 0];                         %窗函数的理想滤波器幅度

dev = [0.01 0.05];              %纹波

[n,wn,beta,ftype]= kaiserord(fc,mag,dev,fs)            %获取凯塞窗函数参数

fpm = [0 fc(1)2/fs fc(2)2/fs 1];                                %firpm函数的频段向量

magpm=[1 1 0 0];                                                    %firpm函数的幅值向量

%设计凯塞窗及海明窗滤波器

h_kaiser =fir1(n,wn,ftype,kaiser(n+1,beta));

h_hamm = fir1(n,fc(2)2/fs);

%设计最优滤波器

h_pm=firpm(n,fpm,magpm);

%求滤波器的幅频响应

m_kaiser =20log(abs(fft(h_kaiser,1024)))/log(10);

m_hamm =20log(abs(fft(h_hamm,1024)))/log(10);

m_pm =20log(abs(fft(h_pm,1024)))/log(10);

%设置幅频响应的横坐标单位为Hz

x_f =[0:(fs/length(m_kaiser)):fs/2];

%只显示正频率部分的幅频响应

m1 =m_kaiser(1:length(x_f));

m2 =m_hamm(1:length(x_f));

m3 =m_pm(1:length(x_f));

%绘制幅频响应曲线

plot(x_f,m1,’-’,x_f,m2,’-.’,x_f,m3,’–’);

xlabel(‘频率(Hz)’);ylabel(‘幅度(dB)’);

legend(‘凯塞窗’,‘海明窗’,‘最优滤波器’);grid；

从图中可以看出使用kaiserord函数获得滤波器的阶数为36，截至频率为0.3125，凯塞窗的β值为3.3953从幅频响应曲线可以看出，最优滤波器与凯塞窗滤波器相比，最优滤波器的第一瓣电平约低2.5dB,且阻带纹波相同，而凯塞窗的阻带纹波却逐渐减小。



四、采用FDATOOL工具设计

例：利用FDATOOL工具设计设计一个带通FIR滤波器，带通范围为1000~ 2000Hz，低频过渡带为700~ 1000Hz,高频过渡带为2000~2300Hz,采样频率为8000Hz的等阻带纹波滤波器要求阻带衰减大与60dB。

启动MATLAB软件后，依次单击主界面左下方的“start→Toolboxes→Filter Design→Filter Design&Analysis Tool”即可打开FDATOOL工具界面

1.首先单击FDATOOL界面左下方的滤波器设计（Design filter）按键图标，进入滤波器设计界面。

2. 单击FDATOOL界面上方工具栏上的滤波器参数设置（Filter Specifications）工具，进入滤波器参数设置界面。

3. 单击滤波器响应类型（Response Type）部分的带通滤波器（Bandpass）选项，指定设计带通滤波器。

4. 在设计方法（Design Method）部分的FIR下拉列表框中，选择等纹波（Equiripple）

设计方法。需要注意的是，由4.2节可知，除了凯塞窗可以通过调整β参数来改变滤波器阻带衰减外，其他窗函数的阻带衰减是无法调整的，增加滤波器阶数只能改变滤波器的过渡带性能。最优滤波器设计方法（等纹波滤波器）可通过增加滤波器阶数改善阻带衰减性能。

5.  根据设计要求，设置滤波器截至频率。

6. 设置滤波器阶数（Specify order）参数后，单击FDATOOL界面下方的“Design Filter”按键开始滤波器设计。

7. 观察FDATOOL中的频率响应，调整滤波器阶数，直到满足设计要求。

通过单击“Analysis→Filter→Coefficients”菜单来察看滤波器系数，或通过单击“Targets→XILINX Coefficient(.COE) File”菜单直接生成FPGA所需的滤波器系数配置文件。
ex:设计一个高通最优FIR滤波器，过渡带为1000~1500Hz,采样频率为8000Hz，通带纹波最大为0.01，阻带纹波最大为0.001绘制比较系数12位、14位量化以及无量化是的幅度响应曲线，将14位量化结果转换成十六进制补码格式，并将14位量化的十进制存放在文本文件中。

设计一个对纹波及截至频率有确定要求的滤波器，需要使用firpm函数设计最优滤波器。而firpm函数无法确定滤波器的阶数，因此需要使用kaiserord函数先求取满足设计需求的最小滤波器阶数。

function

hn=fir_1;

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

h_pm14 = round(h_pm/max(abs(h_pm))(2^13-1));

hn=h_pm14;

%转换进制

q14_hex_pm =dec2hex(h_pm14+2^14*(h_pm14<0));

%将生成的滤波器系数写入FPGA所需的文件中

fid=fopen(‘f:\lbq\fir_1.txt’,‘w’);

fprintf(fid,’%8d\r\n’,h_pm);

fprintf(fid,’;’);

fclose(fid);

%求滤波器的幅频响应

m_pm = 20log(abs(fft(h_pm,1024)))/log(10);

m_pm = m_pm-max(m_pm);

m_pm12 =20log(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;

matlab 运行结果如图

IIR滤波器的MATLAB与FPGA的设计实现
博主微信：flm13724054952，不懂的有疑惑的也可以加微信咨询，欢迎大家前来投稿，谢谢！


引言
数字滤波器是数字信号系统里面最常用、最重要、最基本的元件之一。无限脉冲响应(Infinite Impulse Response,IIR)滤波器具有很高的滤波效率，在相同的幅频响应条件下，所要求的IIR滤波器阶数会明显比FIR滤波器的低。也就是说在不介意滤波器的相位，实现同等滤波性能下的IIR滤波器会比FIR滤波器更加节省电路面积，节省电路功耗。在前面的博文中介绍了FIR滤波器的设计实现，本文将介绍无限脉冲响应IIR滤波器(Infinite Impulse Response,IIR)的原理介绍，架构研究，设计实现等。
原理介绍
IIR滤波器的结构式有直接型跟级联型结构，而级联型结构的IIR滤波器更不容易受到系统的极点零点的影响。



从系统函数等可以总结IIR滤波器的几个显著特性就是：

（1）IIR滤波器同时存在不为零的极点和零点。要保证IIR滤波器是稳定不振荡系统，要求系统的极点必须在单位圆内，因为系统的稳定性就是由系统函数的极点决定的。

（2）由于线性相位滤波器(FIR滤波器)所有的零点和极点都是关于单位圆对称的，所有只允许极点位于单位圆的原点。由于IIR滤波器存在不为零的极点，因此只能实现近似的线性相位特性，也就是说如果电路对滤波器的线性相位要求高的话，那IIR滤波器便不可使用了。

（3）在FPGA设计IIR滤波器的时候，由于存在反馈结构，也就是y(n)一直在累加自己前面的y(n-1)，这种无限脉冲响应与我们有限的寄存器长度相矛盾。所以在设计IIR滤波器，其有限字长效应是实际工程设计应该考虑的问题。

从IIR滤波器的差分方程式可以看得出，输出信号是由两部分组成的：第一部分∑x(n-i)b(i)表示对输入信号的延时组成M阶网络，其实相当于前面博文介绍的FIR滤波器的结构，实现了系统的零点。第二部分 ∑y(n-j)a(j)表示对输出信号的延迟，组成J阶网络，每阶经延时抽头后与常系数相乘，并将乘法结果相加，由于这部分是对输出延时，故为反馈网络，实现了系统的极点。

根据系统的差分方程式跟上面的分析后，例如一阶的IIR滤波器差分方程为:y(n)=b0x(n)+b1x(n-1)-a1y(n-1);其一阶系统的架构设计图可以定为以下直接型结构：


MATLAB设计实现
一般来说，IIR滤波器的设计方法分为三种：原型转换设计法，直接设计法以及直接调用MATLAB函数设计法。本文主要介绍的是通过matlab来实现对IIR滤波器的设计实现。在matlab里面设计数字滤波器既可以通过窗函数自己编写程序来设计IIR滤波器，也可以通过调用工具箱来设计IIR滤波器。通过编写程序来设计IIR滤波器比较普遍，滤波器的参数等改动也比较方便，与FPGA做联合功能仿真验证也方便。而直接调用工具箱就直接了当，较节省设计时间。本文将对两种方法进行介绍，一般matlab采用的几个IIR滤波器设计函数有巴特沃斯，切比雪夫，椭圆三种低通滤波器，不过本文设计采用的是切比雪夫函数来设计IIR滤波器。

一个数字滤波器的传输函数表达式为如下所示：



其中H(w)叫做幅频特性 ,Ψ(w)叫做相频特性。幅频特性表示信号通过该滤波器以后频率成分幅度衰减情况，而相频特性反映各频率分量通过滤波器后在时间上的延时情况。
切比雪夫低通滤波器
巴特沃斯低通滤波器的频率特性在通带和阻带内部都是随着频率单调变化的，也就是如果在通带的边缘能够满足指标，在通带的内部肯定超过设计的指标要求从而造成设计滤波器的阶数较高。如果将指标的精度要求均匀分布在整个通带内，或者在分布在阻带内，或者同时分布在通带跟阻带内，则可以设计出满足设计需要，阶数又低的滤波器了。而切比雪夫滤波器正是如此的。

MATLAB用切比雪夫函数设计滤波器的语法为：

[b,a] = cheby1(n,Rp,Wn,‘type’);

b,a   ：分别为滤波器的分子，分母系数因子；

n      ：为滤波器的设计阶数；

Rp   ：为滤波器在通带内的最大衰减db

Wn  ：为滤波器的截止频率，当Wn为[w1,w2]的向量则代表滤波器为带通滤波器；

type：指定设计的滤波器类型，为hige或者slow;

求滤波器的极点零点增益的函数语法为：

[z,p,k] = cheby1(n,Rp,Wn,‘type’);

z,p,k分别为所设计滤波器的零点，极点，增益表达式。

利用切比雪夫函数设计滤波器的matlab工程如下所示，hb,ha分别为量化后的滤波器分子，分母系数；



运行设计的matlab工程可以得到IIR滤波器的输出为如下图所示，及也可得到hb,ha滤波器分子分母系数；





因为滤波器的系数量化位数，运算位长，输出量化都会对滤波器的性能产生影响。所以设计的时候也应该考虑这些因素。具体这里就不做详细介绍了，因为在前面的博文FIR滤波器的设计里面已经做详细介绍了；这里只详细结果：IIR滤波器可以采用14bit的量化位数跟14bit的量化输出，这样子其滤波器性能才是最佳的。


tool工具箱设计法
在matlab的应用程序里面有一个filter design & analysis tool点进去即为快速设计滤波器的工具。操作界面如下图所示，较为简单不做详细介绍。



点analysis里面的filter coefficient即可出现设计滤波器的FIR滤波器系数，通过file-export可以把滤波器系数导出到workplace。


FPGA设计实现
由上面的matlab的设计可得到IIR滤波器的极点零点,如下图所示：



根据上面的极点零点，来设计本文的IIR滤波器的直接型结构，如下图所示：从IIR滤波器的差分方程式可以看得出，输出信号是由两部分组成的：上面部分∑x(n-i)b(i)表示对输入信号的延时组成4阶网络，其实相当于前面博文介绍的FIR滤波器的结构，实现了系统的零点。下面部分 ∑y(n-j)a(j)表示对输出信号的延迟，组成7阶网络，每阶经延时抽头后与常系数相乘，并将乘法结果相加，由于这部分是对输出延时，故为反馈网络，实现了系统的极点。


直接型IIR滤波器
直接型的IIR滤波器的RTL设计图如下所示，分为两部分实现pole的极点设计及zero的零点设计。


极点零点的设计
极点零点的设计其实都是加法器及乘法器的设计堆砌，涉及的原理基本一致，前面也在FIR滤波器一文里面详细描述了，这里不再做详细介绍，以极点的设计为例做简单描述；如下极点函数的RTL设计图所示，对零点函数模块处理后的数据，进行延迟后于相应的系数进行相乘，然后再累加，输出滤波器的数据。


系数乘法器的设计模块
乘法器的设计端口介绍如下图所示：



乘法器的设计原理介绍如下图所示，具体的乘法设计可参考之前的博客章文：基于FPGA的乘法器原理介绍及设计实现，具体博文链接为：

https://blog.csdn.net/weixin_39015789/article/details/101430566



最后再将所有的乘积和进行累加便是对数据流进行低通滤波后的结果了。
直接型IIR滤波器的功能仿真
下图为零点函数的功能仿真波形图，如下所示，data_in为输入的高斯白噪声，依次经过延迟存入data_mem,然后经过与系数的相乘得到data_reg，再进行累加得到零点函数的输出结果x_data。再传到极点函数模块进行处理。

前面经过零点函数处理得到的数据x_data存入到极点函数的data_mem,与相应的系数相乘累加得到mul_sum。再经过最后的截尾处理就可以得到滤波后的最后结果dout。



直接型IIR滤波器的整体功能仿真波形图如下所示：


直接型IIR滤波器的傅里叶变换
对直接型IIR滤波器的RTL设计输出的滤波后数据进行matalb的傅里叶变换，看其设计的滤波效果如下图所示：


级联型IIR滤波器
级联型的IIR滤波器的RTL设计图如下所示，其实每一个级联的滤波器相当于直接型IIR滤波器的设计，只是每一个滤波器的输入都是上一个滤波器的滤波后数据输出。


级联型IIR滤波器的傅里叶变换
如下图所示，直接型IIR滤波器与级联型IIR滤波器的滤波效果如下所示，直接型跟级联型的滤波效果其实相当接近，但是对于系统的稳定性来说，在实际的设计应用里面应该采用级联型滤波器的设计。